PSD from Wiener Khintchine and FFT

Hey guys,

i am pretty new to using matlab and i need a little bit of help right now.

is the PSD of a statistical process x(t)

ist the autocorrelation function and E{ } the expected value

are a fourier pair

Now the Wiener Khintchine theorem states, that:

Using the average as the expectation:

Then the cross spectrum of two random processes gives :

I was trying to verify the above theory in matlab but didnt manage to do so, maybe you can help me detect what went wrong.

c = randn(1,1e6); %created random sequences

x = 5*randn(1,1e6)+c;

y = 5*randn(1,1e6)+c;

L=1e6;

[xy_cor,lags] = xcorr(x,y); %execute the fft on the cross-correlation on x,y

P2xy = abs(fftshift(fft(xy_cor))/1e6);

P1xy = P2xy(1:L/2+1);

P1xy(2:end-1) = 2*P1xy(2:end-1);

figure(1); %plotted the PSD ,length of showed plot shouldnt really matter

plot(f(1:100),10*log10(P1xy(1:100)))

title('Single-Sided Amplitude Spectrum of xy(t) selfmade')

xlabel('f (Hz)')

ylabel('|P1(f)|')

xf = fft(x);

yf = fft(y);

xf_wien = xf.*yf;

figure(2);

plot(1:100,10*log10(xf_wien(1:100)))

So why do figure 1 and 2 differ? As far as i understood the theory they should be the same???

Thanks in advance

Answer 1

[David Goodmanson]

Hi Jan-Niklas

It's easier done with convolution instead of correlatation, so convolution is first here.

fft does circular convolution or correlation. Neither conv nor xcorr do that. Consequently, for the fft approach, the x and y arrays (length n) are padded with n zeros so that the array can't 'come around the other end' and give a nonzero contribution.

x = [1 2 1 3 5 4];
y = [4 2 3 3 1 9];
conv12 = conv(x,y)
% by fft
n = length(x);
nzer = zeros(1,n);
xpd = [x nzer];                          % padded arrays
ypd = [y nzer];
conv12byfft  = ifft(fft(xpd).*fft(ypd))  % agrees with conv12
% in frequency domain
fft([conv12 0])                          % these two complex arrays agree
fft(xpd).*fft(ypd)

The conv12 array has 2n-1 entries and the conv12byfft array has 2n entries, with an extra zero at the end. To compare results in the frequency domain as you are doing, then you have to add a zero at the end of conv12 as shown, before doing the fft.

---> Note the nice symmetry between x and y, where fft applies to both the padded x and y arrays.

For xcorr, things do not work out quite as cleanly. In that case one of the arrays requires an fft, one an ifft. Also there is a circular shift by one place.

[The ifft(ypd) function is multiplied by 2n because ifft automatically divides by the length of the array and that factor has to be canceled out].

x = [1 2 1 3 5 4];
y = [4 2 3 3 1 9];
corr12 = xcorr(x,y)
% by fft
n = length(x);
nzer = zeros(1,n);
xpd = [x nzer];                                % padded arrays          
ypd = circshift([nzer y],1);
corr12byfft = ifft(fft(xpd).*(2*n*ifft(ypd)))  % agrees with corr12
% in frequency domain
fft([corr12 0])                                % these two complex arrays agree
fft(xpd).*(2*n*ifft(ypd))

I have not tried it for n = 1e6 but I'm going on the principle (there must be a name for it) that, time and memory problems aside, if a code works for n= 5 or 6 it will work for n = one billion.

[Fego Etese]

Hello David Goodmanson, sorry to contact you like this, but its really an emergency, please I need your help with this question

https://www.mathworks.com/matlabcentral/answers/552142-how-to-use-wavelet-transforms-for-template-matching-and-location

I will be grateful if you can help me out, thanks

[David Goodmanson]

Hello Fego,

Sorry, I would provide information if I could, but I don't have enough knowledge in this area.