Estimate phase using STFT in Matlab?

Question

My signal oscillates back and forth on the interval [-1, 1] at varying frequencies, with some added higher frequency noise components (example data shown below). I want to to estimate the phase of the slow oscillation at any given time. This means I would like to filter out some of the high frequency noise, and then calculate the phase at each time for the smoothed signal. However, my method isn't working so I'd really appreciate any ideas about how to do this properly!

I've attempted to do this with the Short Time Fourier Transform (FTST). The method I've tried is to calculate the Fourier transform for a sliding window across my signal, choose the frequency component with the highest magnitude, and then use the phase from that component.

However, I'm having trouble getting the result to match up with the slowly oscillating component that I can see in the signal. I've tried two methods and neither looks good:

1. using just the phase shift from the complex value of the fourier transform computed using Matlab's spectrogram() function
2. using the total phase based on the time and frequency return variables from the STFT plus the phase shift

Am I missing something? I'm new to signal processing, so it might be a conceptual misunderstanding or just an error in my implementation. I would really appreciate any help or suggestions!

Here is the Matlab code I'm using to generate the figure:

% 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)

Answer 1

Jumping between bins of STFTs will lead to the discontinuities seen in your phase estimates. Interpolating both frequency and phase from a sequence of STFTs is possible, but requires non-standard techniques, such as fftshifts and 2D interpolators.

If you know the frequency range of your desired slow oscillation, then a bandpass filter in conjunction with a PLL follower might provide a better looking phase result using more standard DSP techniques.

Answer 2

Filter the signal using a bandpass filter centred on the frequency of interest, then create the analytic signal. In MATLAB this is:

signal_analytic = hilbert(signal);

The amplitude and phase of the signal are then given by:

amplitude = abs(signal_analytic);
phase = angle(signal_analytic);