Reference level P0 in psd() logarithm?

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

Answer 1

[Wayne King]

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:

1. use aryule() to obtain the AR coefficients
2. return the error variance along with the AR coefficients from aryule
3. use freqz() to estimate the frequency response
4. 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()
5. scale appropriately with the sampling frequency (if in Hz) or (2*pi) if in radians/sample.

[Bob Li]

Wayne,

Thanks for the answer. I agree that reference level is not critical here. I will look back to this later.

Bob

Answer 2

[Wayne King]

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().

[Bob Li]

Wayne,

Thanks for this example. I have noted it down and will study it.

Bob