当前位置：首页 > news >正文

MATLAB中基于两步法的时频脊线提取与瞬时频率估计新算法

news 2026/7/27 2:13:21

MATLAB环境下一种基于两步法的瞬时频率（IF）估计算法算法运行环境为MATLAB r2018a，为一种更优秀的瞬时频率（IF）估计方法，可用于信号的时频脊线提取算法可迁移至金融时间序列，地震信号，语音信号，声信号，生理信号（ECG,EEG,EMG）等一维时间序列信号

最近捣鼓了个MATLAB下的瞬时频率（IF）估计算法，叫“两步法”，感觉挺有意思的。为啥突然研究这个？主要是之前处理非平稳信号时，传统方法要么精度不行，要么计算太慢，这个两步法据说能兼顾效率和准确性，还能套用到金融、地震这些一维信号上，今天就给大家唠唠怎么玩的。

啥是瞬时频率？为啥这么重要？

简单说，瞬时频率就是信号在某一时刻的“瞬时振动频率”。比如你听音乐时，基频随时间变化，这个变化的频率就是IF；地震信号里，我们想知道震源能量释放的频率变化；金融数据里，股价波动频率能反映市场情绪；语音、脑电这些生理信号，更是靠IF分析不同状态。

但传统IF估计有坑啊！比如直接对信号求相位导数，遇到噪声大的信号，结果就像喝了假酒——飘得没边儿。或者用短时傅里叶变换（STFT），时频分辨率互相“打架”：窗口小了时间准但频率糊，窗口大了频率准但时间模糊。这个两步法，就是想解决这些问题。

两步法咋玩？先粗估再精修！

核心思路特简单：第一步快速粗估，抓住信号整体趋势；第二步精细优化，把毛刺去掉。以最常见的“带噪声线性调频信号（chirp）”为例，它的瞬时频率是线性变化的，适合验证算法。

1. 造个测试信号

先模拟一个带噪声的chirp信号（线性调频，瞬时频率从100Hz线性增加到105Hz，1秒时长，采样率1000Hz）：

Fs = 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)); % 加高斯白噪声（信噪比10dB）

2. 第一步：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); % 滑动平均→精估IF

4. 对比结果

画图看看效果：

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');

运行后：