Рассчитать автокорреляцию с использованием БПФ в Matlab

Question

Я прочитал некоторые объяснения того, как автокорреляция может быть более эффективно рассчитана, используя fft сигнала, умножая действительную часть на комплексное сопряжение (область Фурье), затем используя обратное fft, но у меня возникают проблемы с реализацией это в Matlab, потому что на детальном уровне.

Answer 1

Как вы и заявили, возьмите fft и умножьте поточечно на его комплексное сопряжение, затем используйте обратную fft (или в случае взаимной корреляции двух сигналов: Corr(x,y) <=> FFT(x)FFT(y)*)

x = rand(100,1);
len = length(x);

%# autocorrelation
nfft = 2^nextpow2(2*len-1);
r = ifft( fft(x,nfft) .* conj(fft(x,nfft)) );

%# rearrange and keep values corresponding to lags: -(len-1):+(len-1)
r = [r(end-len+2:end) ; r(1:len)];

%# compare with MATLAB's XCORR output
all( (xcorr(x)-r) < 1e-10 )

На самом деле, если вы посмотрите на код xcorr.m, это именно то, что он делает (только он имеет дело со всеми случаями заполнения, нормализации, векторного / матричного ввода и т. Д ...)

Answer 2

По теореме Винера – Хинчина спектральная плотность мощности (PSD) функции является преобразованием Фурье автокорреляции. Для детерминированных сигналов PSD является просто квадратом амплитуды преобразования Фурье. См. Также теорема о свертке .

Когда речь идет о дискретных преобразованиях Фурье (т. Е. С использованием БПФ), вы фактически получаете циклическую автокорреляцию. Чтобы получить правильную (линейную) автокорреляцию, вы должны обнулить исходные данные до удвоенной исходной длины перед выполнением преобразования Фурье. Так что-то вроде:

x = [ ... ];
x_pad = [x zeros(size(x))];
X     = fft(x_pad);
X_psd = abs(X).^2;
r_xx = ifft(X_psd);