- use aryule() to obtain the AR coefficients
- return the error variance along with the AR coefficients from aryule
- use freqz() to estimate the frequency response
- use e*abs(h).^2 as the basis for PSD estimate where e is the prediction error from aryule() and h is the frequency response from freqz()
- scale appropriately with the sampling frequency (if in Hz) or (2*pi) if in radians/sample.
Reference level P0 in psd() logarithm?
1 view (last 30 days)
Show older comments
Hi,
I have written code to implement Yule-Walker algorithm myself and compared it against psd() result In my code, I directly take the log10 of of the squared DTFT, and found a constant difference around 100 that my value is above psd() value.
According to http://en.wikipedia.org/wiki/Decibel the definition of power decibel requires a reference level P0, and I therefore suspect it is this P0 that lowered psd() value. Could anyone tell me the value of this P0?
Bob
0 Comments
Accepted Answer
Wayne King
on 14 Jan 2012
Hi Bob, In PSD estimates, you can take this to be 1.
The use of the logarithm in PSD calculations is different than the use of the logarithm in measurements like sound pressure level where the reference level is critical.
The reason you typically take the log in PSD estimation is for variance stabilization. The variance in a PSD estimate like the periodogram is frequency-dependent. Taking the log() removes that frequency dependence and makes the PSD estimate much easier to interpret.
Have you looked at aryule() in the Signal Processing Toolbox?
Basically the Yule-Walker spectral estimators:
More Answers (1)
Wayne King
on 14 Jan 2012
You're welcome Bob, just to reiterate what I stated above, if you try the following:
x = randn(1000,1);
[a,e] = aryule(x,4);
[h,w] = freqz(1,a,512);
Pxx = e*abs(h).^2;
Pxx(2:end-1) = 2*Pxx(2:end-1);
Pxx = Pxx./(2*pi);
plot(w./pi,10*log10(Pxx),'k'); grid on;
set(gca,'xlim',[0 1]);
figure;
pyulear(x,4,512);
You'll see what I mean about the scaling. aryule() obtains the biased estimate of the autocorrelation and then solves the Yule-Walker equations using levinson().
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!