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

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

Ответ 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 XCORR output
all( (xcorr(x)-r) < 1e-10 )

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

Ответ 2

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

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

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