Prev: variable object names
Next: how to clear the cache (rehash clear mex ???)
From: dbd on 12 Aug 2010 16:04
On Aug 12, 2:42 am, "Lucila " <lucpat...(a)gmail.com> wrote:
> Hello,
> I am generating a random process by passing white noise throught a second order system.
> a/(s^2+2*z*wn*s+wn^2)
> The autocorrelation function calculated using xcorr gives me what I expected, but the psd (calculated using fft) does not.
> I expected to see a decreasing line with slope -4 (loglog-Plot)

> If I increase the number of samples or the magnitude of the natural frequency (wn), I can see the correct slope.

You cannot see the correct slope in either of your plots.

> Why?

You have greatly abused poor old pwelch()

> Thank you and Regards,
> Lucy

1) Your plots default to different fft sizes with no averaging. You
are plotting different points in your plots. pwelch will average for
you. Give it a transform size and enough samples to average.

2) You are defaulting to a rectangular window. When you are looking
for steep rolloff of a psd, you need to use a window that has a
sidelobe rolloff that is faster than the slope you seek to observe.
Otherwise, 'spectral leakage' from your low frequency energy will
obscure your high frequency response as in your case. Hann rolls off
faster that rectangular, point-wise squaring of Hann will roll off
even faster. Note that for fft application you should use the
'periodic' flag for the built-in windows.

Example:

nData=1e4;
n=randn(1,nData); % Input signal
Fs = 200; %Hz
zeros_plano_S = [];
polos_plano_S = [-0.01+0.003i -0.01-0.003i]; %%%%%
[b,a] = zp2tf(zeros_plano_S,polos_plano_S,1); %
[numd,dend] = impinvar(b,a,Fs);
N2 = filter(numd,dend,n);%%%%%
[Pxx1,freq]=pwelch(N2,hann(10000,'periodic').^2,[],10000,Fs);

nData=3e4;
n=randn(1,nData); % Input signal
N2 = filter(numd,dend,n); %%%%%
[Pxx2,freq]=pwelch(N2,hann(10000,'periodic').^2,[],10000,Fs);

nData=1e6;
n=randn(1,nData); % Input signal
N2 = filter(numd,dend,n); %%%%%
[Pxx3,freq]=pwelch(N2,hann(10000,'periodic').^2,[],10000,Fs);

figure, loglog(freq,Pxx1,'r-',freq,Pxx2,'g-',freq,Pxx3,'b-');
axis([1e-2 1e2 1e-12 1e2]), grid on, zoom on
title(['Properly Windowed Plots'])
xlabel(['Unaveraged: Red Averaged by 3: Green Averaged by 100:
Blue'])
figure, loglog(freq,Pxx1,'r-',freq,Pxx2,'g-',freq,Pxx3,'b-')
axis([0.8*1e-0 4*1e1 0.5*1e-10 2*1e-6]), grid on, zoom on
title(['Properly Windowed Plots'])
xlabel(['Unaveraged: Red Averaged by 3: Green Averaged by 100:
Blue'])

Dale B. Dalrymple
From: Walter Roberson on 12 Aug 2010 16:26
sudipta pramanik wrote:
> sir, I have two psd which are certainly different.Give a parameter by
> which i can distinctly identify them

Goedelize(psd) is guaranteed to give different scalar numbers for different psd's.

You will need either a symbolic toolkit or the FEX contribution for indefinite
precision arithmetic in order to implement this.
First  |  Prev  | 
Pages: 1 2
Prev: variable object names
Next: how to clear the cache (rehash clear mex ???)