I'm trying to calculate an autocorrelation on a platform where the only accelerated primitive I have available is the (I)FFT. I'm having a problem though.
I prototyped it in MATLAB. I am, however, slightly confused. I assumed that it works simply as follows (this is from memory so apologies if I've got it slightly wrong).
autocorr = ifft( complex( abs( fft( inputData ) ), 0 ) )
However I get a different result than I get from using the xcorr
function. Now I'm fully expecting not to get the left hand side of the auto correlation (as it's a reflection of the right hand side and thus not needed anyway). However, the problem is my right hand side appears to be, itself, reflected around the halfway point. Which effectively means that I get about half the amount of data that I'm expecting.
So I'm sure I must be doing something very simple wrong, but I just can't figure out what.
Answer
pichenettes is right, of course. The FFT implements a circular convolution while the xcorr() is based on a linear convolution. In addition you need to square the absolute value in the frequency domain as well. Here is a code snippet that handles all the zero padding, shifting & truncating.
%% Cross correlation through a FFT
n = 1024;
x = randn(n,1);
% cross correlation reference
xref = xcorr(x,x);
%FFT method based on zero padding
fx = fft([x; zeros(n,1)]); % zero pad and FFT
x2 = ifft(fx.*conj(fx)); % abs()^2 and IFFT
% circulate to get the peak in the middle and drop one
% excess zero to get to 2*n-1 samples
x2 = [x2(n+2:end); x2(1:n)];
% calculate the error
d = x2-xref; % difference, this is actually zero
fprintf('Max error = %6.2f\n',max(abs(d)));
No comments:
Post a Comment