# Comparison of PSD calculating methods.

36 Ansichten (letzte 30 Tage)
Alex Bobrovnik am 20 Nov. 2011
Hello everyone. I have some problems with my Matlab code. I need to compare 3 methods of calculating PSD:
1. Welch-method;
2. FFT of autocorrelation function;
3. Averaging of frequency bands;
close all;
clear all;
clc;
fs=16000;
t=0.1;
b1=rand(1,fs*t);
1)
% PSD method 1 (welch)
figure(1)
h=spectrum.welch;
psd(h,b1)
% it works correctly, here i use some different windows and so on
figure(2)
wlen=256 ;
h=spectrum.welch('hamming',wlen,0.5)
psd(h,b1)
2)
%PSD method 2 (xcorr)
rxx=xcorr(b1); %getting ACF
PSD=abs(fft(rxx));
figure(3);
plot(PSD);
• i don't think that it's right, because of strange shape of PSD;
• how to make ACF symmetrically of "0"?
• if i want to "cut" ACF by the "triangular" window, what should i do?
3) ...
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Akzeptierte Antwort

Wayne King am 21 Nov. 2011
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+randn(size(t));
Rxx = xcorr(x,'biased');
Rxxdft = abs(fftshift(fft(Rxx)));
freq = -Fs/2:Fs/length(Rxx):Fs/2-(Fs/length(Rxx));
plot(freq,Rxxdft);
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Weitere Antworten (3)

Alex am 20 Nov. 2011
xcorr adds some padding numbers to either side. I forget how to manually truncate the series to take out the side bands (and I don't have Matlab at home).
So, here is one idea: calculate the auto correlation yourself. The function is pretty easy to implement with a few for loops. However, this method will take more processing time (Matlab's xcorr function using the fft & ifft to calculate it, which is much faster).
Second, do you know if your auto-corr should be biased or not? that will give you slightly different vairiations.
Third, the plot((abs(fft)) does not plot 0 at the center.
lets say x = -4:5, then plot(abs(fft(ifft(x))) will end up plotting (1,2,3,4,5,-4,-3,-2,-1).
to get the plot to come out right, use something like this:
rez = abs(fft(xcorr))
rez = circshift(rez',round(length(rez)/2))
plot(rez)
##### 5 Kommentare3 ältere Kommentare anzeigen3 ältere Kommentare ausblenden
Alex Bobrovnik am 21 Nov. 2011
i must use both of them)
Alex Bobrovnik am 21 Nov. 2011
and about windowing... "cutting" ACF gives us smoother estimation, isn't it? I want to show this effect. I mean if i'll use 64 points of ACF or 128 points it'll give different shape.

Melden Sie sich an, um zu kommentieren.

Wayne King am 21 Nov. 2011
Hi Alex, If the time series is zero-mean, the periodogram is equivalent to Fourier transform of the biased autocorrelation sequence.
##### 4 Kommentare2 ältere Kommentare anzeigen2 ältere Kommentare ausblenden
Alex am 21 Nov. 2011
According to http://en.wikipedia.org/wiki/Spectral_density, it is Ok to take calculate the PSD of a non-zero mean signal by using the FFT(auto corr) if the signal is Wide Sense Stationary (WSS).
Both of the, Uniform and Normal, distributions used in this project are WSS
Wayne King am 21 Nov. 2011
yes, it's fine to have the psd of a nonzero mean process, I never said it wasn't. All I said was, if you are comparing PSD methods, then you can easily make your process zero mean, that's not going to affect any comparison of PSD methods, and then you can use the periodogram which is the same as the Fourier transform of the biased autocorrelation sequence.

Melden Sie sich an, um zu kommentieren.

Alex Bobrovnik am 21 Nov. 2011
Today I tried to improve my code:
... % % PSD method 2 (xcorr) % Calculate the Power Spectral Density of Original Signal%
Rxx=xcorr(b1,b1,length(b1)-1);
Pxx=abs(fftshift(fft(Rxx)));
figure(2)
subplot(211)
plot(Rxx)
title('Autocorrelation function');
subplot(212)
plot(Pxx)
title('PSD');
...
but i'm not sure... i can't understand why ACF isn't symmetric on 0
##### 2 KommentareKeine anzeigenKeine ausblenden
Alex am 21 Nov. 2011
For the Future, please update your original post with the updated information.
Alex am 21 Nov. 2011
I just ran your code, and I'm getting what I expect, the spikes centered at 0.
Going back to what you mentioned in commenting on my post earlier, you can see the differences between using different # of points with:
Rxx64 = xcorr(b1,64,'biased');
Rxx128 = xcorr(b1,128,'biased');

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Correlation and Convolution finden Sie in Help Center und File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by