An FFT of a signal returns a complex number, with a phase angle related to the starting location of each frequency within the signal. If I have break my signal into blocks, and FFT them individually, the phase angle is changing in time. Is there a smart way to 'average' these phase angles - I know a direct average over time will give you a 0, if you have enough points.
I can offset the phase by the amount of my window shift, but I'm wondering if there is anything more robust, that can handle noisy signals. Also, if the FFT bin width is not sufficiently accurate, the offset will be slightly off and hence won't correct properly.
Here is a block of code that shows what I've done so far:
fs = 1000; % sampling frequency fp = 5; % peak frequency t = [1/fs : 1/fs : 100]'; % time signal x = sin(t * fp * 2 * pi + rand(length(t),1)); % data signal % FFT block sizes: N = 2^11; Ns = 2^4; % Split data into blocks: for ii = 1:ceil((length(x)-N+1)/Ns) % loop for each window xx(:,ii) = x(((ii-1)*Ns+1):((ii-1)*Ns+N)); tt(:,ii) = t(((ii-1)*Ns+1):((ii-1)*Ns+N)); end % Plot first five blocks for jj = 1:5 subplot(5,1,jj); plot(tt(:,jj),xx(:,jj)); title(sprintf('Block %g',jj)); end % FFT each block window_fn = hanning(N) * ones(1,size(xx,2)); % window function yy=window_fn.*(xx-ones(N,1)*(mean(xx))); %remove mean and window each block % Plot first five blocks for jj = 1:5 subplot(5,1,jj); hold on; plot(tt(:,jj),yy(:,jj)); ylim([-1 1]); xlim([0 t(N+Ns*4)]); end nfft = N * 4; % zero pad fftx = fft(yy,nfft); % Find peak frequency from mean of power of each window mx = abs(fftx).^2 / N/fs; NumUniquePts = ceil((nfft+1)/2); f = (0:NumUniquePts-1)*fs/nfft; fftx = fftx(1:NumUniquePts,:); [~,peakf] = max(mean(mx,2)); % Get phase information from each block; phi = fftx(peakf,:); figure; plot(angle(phi),'s-'); xlim([0 50]); %unwrap is helpful too % rotate the complex number by the phase offset (from Ns) shift = (0:(length(phi)-1)) * (Ns/fs) * f(peakf); newphi = phi .* exp(1i * shift*-2*pi); hold on; plot(angle(newphi),'.-'); xlabel('First 50 blocks'); ylabel('$\phi$');