MATLAB环境下一种基于两步法的瞬时频率IF估计 算法 算法运行环境为MATLAB r2018a为一种更优秀的瞬时频率IF估计方法可用于信号的时频脊线提取 算法可迁移至金融时间序列地震信号语音信号声信号生理信号ECG,EEG,EMG等一维时间序列信号最近捣鼓了个MATLAB下的瞬时频率IF估计算法叫“两步法”感觉挺有意思的。为啥突然研究这个主要是之前处理非平稳信号时传统方法要么精度不行要么计算太慢这个两步法据说能兼顾效率和准确性还能套用到金融、地震这些一维信号上今天就给大家唠唠怎么玩的。啥是瞬时频率为啥这么重要简单说瞬时频率就是信号在某一时刻的“瞬时振动频率”。比如你听音乐时基频随时间变化这个变化的频率就是IF地震信号里我们想知道震源能量释放的频率变化金融数据里股价波动频率能反映市场情绪语音、脑电这些生理信号更是靠IF分析不同状态。MATLAB环境下一种基于两步法的瞬时频率IF估计 算法 算法运行环境为MATLAB r2018a为一种更优秀的瞬时频率IF估计方法可用于信号的时频脊线提取 算法可迁移至金融时间序列地震信号语音信号声信号生理信号ECG,EEG,EMG等一维时间序列信号但传统IF估计有坑啊比如直接对信号求相位导数遇到噪声大的信号结果就像喝了假酒——飘得没边儿。或者用短时傅里叶变换STFT时频分辨率互相“打架”窗口小了时间准但频率糊窗口大了频率准但时间模糊。这个两步法就是想解决这些问题。两步法咋玩先粗估再精修核心思路特简单第一步快速粗估抓住信号整体趋势第二步精细优化把毛刺去掉。以最常见的“带噪声线性调频信号chirp”为例它的瞬时频率是线性变化的适合验证算法。1. 造个测试信号先模拟一个带噪声的chirp信号线性调频瞬时频率从100Hz线性增加到105Hz1秒时长采样率1000HzFs 1000; t (0:Fs-1)/Fs; % 时间轴0到1秒 f0 100; % 初始频率 beta 5; % 频率斜率每秒增加5Hz true_IF f0 beta*t; % 真实瞬时频率 signal chirp(t, f0, 1, f0 beta*1, linear); % 无噪chirp signal signal 0.5*randn(size(t)); % 加高斯白噪声信噪比10dB2. 第一步Hilbert初估快速抓趋势用Hilbert变换把实信号转为解析信号通过相位导数得到IF的“粗轮廓”。这一步快但噪声大hilbert_signal hilbert(signal); % 解析信号实部虚部 phase angle(hilbert_signal); % 瞬时相位-π到π phase unwrap(phase); % **相位展开**避免相位跳变导致IF错误 phase_diff diff(phase); % 相位差 time_diff diff(t); % 时间差 if_estimate1 phase_diff / time_diff; % 直接相位导数→粗估IF为啥要unwrap比如相位从π突然跳到-π差值会是-2π而实际应该是0。unwrap能把相位拉成连续的避免“假跳变”。3. 第二步滑动平均精修去毛刺第一步的IF有大量高频噪声第二步用滑动平均平滑毛刺保留趋势window_size 5; % 窗口大小可根据噪声调整比如噪声大就设大一点 if_estimate2 movmean(if_estimate1, window_size); % 滑动平均→精估IF4. 对比结果画图看看效果figure; subplot(2,1,1); plot(t, signal); title(带噪声的chirp信号); subplot(2,1,2); plot(t(1:end-1), if_estimate1, b, LineWidth,1); hold on; plot(t(1:end-1), if_estimate2, r, LineWidth,2); plot(t, true_IF, k--, LineWidth,2); legend(初估IF,精估IF,真实IF);运行后蓝色初估IF毛刺极多但能看出“100→105Hz”的趋势红色精估IF滑动平均后毛刺被平滑几乎和真实的黑色虚线重合两步法的优势快Hilbert相位导数计算量小适合实时处理抗噪滑动平均简单高效去掉高频噪声普适只要是一维时间序列金融、地震、语音、生理信号都能套。小Tips窗口大小噪声大→窗口大比如10噪声小→窗口小比如3适用场景线性/非线性变化的信号都能试比如非线性调频、带谐波的复杂信号。如果你的项目刚好需要分析一维信号的时频脊线这个两步法可以直接改改参数用起来代码简单效果比传统方法稳不少~ 有兴趣可以跑一下测试信号对比下结果PS下期试试非平稳语音信号看看能不能用这个方法提取基频变化~