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:
- using just the phase shift from the complex value of the fourier transform computed using Matlab's spectrogram() function
- 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)