PSD from Wiener Khintchine and FFT

조회 수: 15 (최근 30일)
Jan-Niklas Krogmann
Jan-Niklas Krogmann 2020년 6월 17일
댓글: David Goodmanson 2020년 6월 21일
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

답변 (1개)

David Goodmanson
David Goodmanson 2020년 6월 19일
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.
  댓글 수: 2
Fego Etese
Fego Etese 2020년 6월 21일
Hello David Goodmanson, sorry to contact you like this, but its really an emergency, please I need your help with this question
I will be grateful if you can help me out, thanks
David Goodmanson
David Goodmanson 2020년 6월 21일
Hello Fego,
Sorry, I would provide information if I could, but I don't have enough knowledge in this area.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Correlation and Convolution에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by