信息处理 - 在 Matlab 中使用 STFT 估计相位？ - 吾爱随笔录

我的信号在区间 [-1, 1] 上以不同的频率来回振荡，并添加了一些更高频率的噪声分量（示例数据如下所示）。我想估计任何给定时间慢振荡的相位。这意味着我想过滤掉一些高频噪声，然后计算平滑信号的每次相位。但是，我的方法不起作用，所以我非常感谢有关如何正确执行此操作的任何想法！

我试图用短时傅里叶变换（FTST）来做到这一点。我尝试过的方法是计算信号中滑动窗口的傅里叶变换，选择幅度最高的频率分量，然后使用该分量的相位。

但是，我无法让结果与我在信号中看到的缓慢振荡的分量相匹配。我尝试了两种方法，但都不好看：

仅使用从使用 Matlab 的 spectrogram() 函数计算的傅立叶变换的复值的相移
使用基于来自 STFT 的时间和频率返回变量加上相移的总相位

我错过了什么吗？我是信号处理的新手，所以这可能是一个概念上的误解，或者只是我的实现中的一个错误。我真的很感激任何帮助或建议！

我的数据示例和我从低频振荡中提取相位的失败尝试

这是我用来生成图形的 Matlab 代码：

% define STFT options
window = 50;
num_overlap = window - 1; % because we need phase at every position
signal_length = length(signal);

% pad signal with zeros on both sides
pad_left = floor(num_overlap/2);
pad_right = num_overlap - pad_left;
padded_signal = [zeros(1,pad_left) signal zeros(1,pad_right)];

% calculate STFT
[S,F,T] = spectrogram(padded_signal,window,num_overlap);
component_magnitude = abs(S); 
[m max_index] = max(component_magnitude); % get the frequency with maximal amplitude

% calculate the phase when the frequency components are maximized
for ti = 1:size(S,2)
   phase_shift(ti) = phase(S(max_index(ti), ti));
   signal_phase_total(ti) =  T(max_index(ti))*F(max_index(ti))+phase_shift(ti);
end  

num_plots = 5;
% plot the initial signal
subplot(num_plots,1,1)
plot(signal, 'LineWidth', 2);
title('signal', 'FontSize', 18)
axis tight

% plot the STFT
subplot(num_plots,1,2)
imagesc(component_magnitude); 
hold on
plot(1:length(max_index), max_index, 'y', 'LineWidth', 2) 
hold off
colormap(bone); 
title('STFT magnitudes with max component in yellow', 'FontSize', 18)

% plot the phase info
subplot(num_plots,1,3)
plot(phase_shift, 'k');
title('phase shift of maximal component', 'FontSize', 18)
axis tight

% try to superimpose the phase shift onto original signal
subplot(num_plots,1,4)           
hs = plot(signal, 'LineWidth', 2);
hold on
% plot the phase of the max component
hp = plot(real(cos(phase_shift)), 'r');
lh = legend([hs hp], 'Signal', 'Phase');
set(lh, 'FontSize', 12);
axis tight
hold off
title('phase shift of max component superimposed onto signal', 'FontSize', 18)        

% try to superimpose the total phase onto original signal
subplot(num_plots,1,5)           
hs = plot(signal, 'LineWidth', 2);
hold on
% plot the total phase
hp = plot(real(cos(signal_phase_total)), 'r');
lh = legend([hs hp], 'Signal', 'Phase');
set(lh, 'FontSize', 12);
axis tight
hold off
title('total phase superimposed onto signal', 'FontSize', 18)