0
$\begingroup$

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:

enter image description here enter image description here

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$');
$\endgroup$
  • $\begingroup$ Phase relative to what? For what relationship do you need phase information? Also, FFT bin phase is not that of all frequencies, but only those equal to basis vectors. $\endgroup$ – hotpaw2 May 2 '17 at 5:52
  • $\begingroup$ The phase relative to its (arbitrary maybe? it looks like phi=0 is the cosine) reference from the FFT. What do you mean the bin phase is not that of all frequencies? $\endgroup$ – James May 3 '17 at 7:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.