**POWER SPECTRAL
DENSITY FUNCTION**

**DEFINITION**

As
with deterministic signals, the frequency content of random signals is also
very important. The definition is different and is a function of how the
variation in a signal is caused by different frequency components. Thus it is
strictly a *variance density function*
that describes how the signal energy or power is distributed across frequency.
It is more traditionally called the *Power
Spectral Density (PSD) Function*. Refer to Chapter 6 for the formal
definition.

**PSD BASED ON THE AR MODEL**

**Derivation**

One of
the primary methods to derive the PSD is to develop an AR model of the signal
and then to calculate it through the *power
transfer function*. Let's go through this procedure by deriving the PSD of a
signal. Consider an AR(2) signal defined below with T = 10 sec and a_{0} = 1, a_{1}
= -1, a_{2} = 0.5. The input white noise signal will have a variance of
1. Notice that the function *freqz* is
useful for calculating points on the transfer function.

clear; close all; format compact; whitebg('w');

T = 10; %
sampling interval, seconds

%
generate some points of the signal model

N
= 200;

w
= randn(250,1); y = zeros(N,1);

t
= T*[0:N-1]';

y(1)
= w(1); y(2) = w(2);

for
n = 3:N

y(n) = y(n-1) - 0.5*y(n-2) + w(n);

end

%
now calculate the PSD of the signal model

a
= [1 -1 0.5];

[h
rf] = freqz(1,a,101);

psdy
= (1*T)*abs(h) .^2; fm = rf/(2*pi*T);

figure(1)

subplot(2,1,1);plot(t,y);
title('SIMULATED SIGNAL')

xlabel('TIME,
seconds'); ylabel('AMPLITUDE')

subplot(2,1,2);plot(fm,psdy);
title('THEORETICAL PSD')

xlabel('FREQUENCY, Hz'); ylabel('MAGNITUDE')

NOTE: here is the direct way to calculate the PSD if the 'freqz' command is not available

fm = (0:100)'/(100*2*T);

h = ones(101,1) ./(1 -1*exp(-i*2*pi*fm*T) +0.5*exp(-i*2*pi*fm*2*T));

psdy = (1*T)*abs(h) .^2;

**Estimation**

The
estimation of the PSD of a random signal is very straight-forward if one is
using a modeling approach. One simply estimates the model as was learned before
and then calculate the PSD based on the power transfer function concept. Being
consistent with the modeling approach*,
the variance of the input signal is equal to the variance of the modeling error*.
Now let's do this for the signal simulated above.

ai
= [1 th(3,1:2)]; %estimates of AR model parameters

disp('ai
= '); disp(ai)

[h
rf] = freqz(1,ai,101);

varerr
= th(1,1);

psdye
= (varerr*T)*abs(h) .^2; fm = rf/(2*pi*T);

figure(2)

plot(fm,psdy,fm,psdye);
title(' PSD & ESTIMATE')

xlabel('FREQUENCY,
Hz'); ylabel('MAGNITUDE');

legend('PSD','ESTIMATE')

__How close are the estimated and theoretical PSDs?
Simulate 2 more signals and estimate their PSDs and plot them. How close are
they to the theoretical one? How close are the AR coefficients to the
theoretical ones? Model the signal with a tenth order model and estimate the
PSD. What has happened and what is it called?__

Now
let's estimate the PSD of a signal with which we are familiar, the grinding
wheel data. Below is the estimated PSD using a fifth order model. Notice that
upper and lower confidence bounds are placed on the plot. __What do they mean?__
Please refer to Section 8.4.2. Some of you thought that a second order model
was good for the signal. __Estimate the PSD with a second order model, plot
it, and compare the two spectra.__

gr
= detrend(grnwheel); ord = 5; N = length(gr);

th
= ar(gr,ord,'burg');

ae
= [1 th(3,:)];

[h
rf] = freqz(1,ae,101);

varerr
= th(1,1);

psdgr
= (varerr*T)*abs(h) .^2; fm = rf/(2*pi*T);

psdgrup
= psdgr*(1+1.96*sqrt(2*ord/N));

psdgrlow
= psdgr*(1-1.96*sqrt(2*ord/N));

figure(3)

plot(fm,psdgr,fm,psdgrup,fm,psdgrlow);
title(' PSD & BOUNDS')

xlabel('CYCLES
PER INCH'); ylabel('MAGNITUDE');

legend('PSD ESTIMATE','UPPER BOUND','LOWER BOUND')

**PSD - NONPARAMETRIC
APPROACH**

**Definition**

There is another definition of the PSD that is based
directly on the DFT. This is the classical definition that is derived as the
Fourier Transform of the autocorrelation function. This fundamental definition
was presented in section 6.2. An alternative mathematical form of this definition
of the PSD is called the *periodogram*.
It and its estimate are found in equations 6.8 and 7.10 respectively. The
periodogram is now used extensively to estimate the PSD because it is simply
proportional to the square value of the DFT; refer to equation 7.19.

**Estimation**

The 'downside' of this estimator is that it is *not consistent*; that is, directly used
it is not accurate and that an additional operation must be applied to it. We
shall examine these concepts now by starting with a signal whose PSD is well
defined, white noise.

N = 1000;

x = randn(N,1); T=1; % Gaussian white noise

t = T*(1:N);

var = cov(x); ave = mean(x);

disp('average and variance'); disp([ave var])

%

% periodogram estimate

xd = detrend(x);

h = hanning(N);

xdw = xd .* h;

PL = 0.36; % process loss factor

dft = T*fft(xdw);

I = (abs(dft(1:N/2)) .^ 2)/(N*T* PL);

f = (0:(N/2)-1)/(N*T);

figure(4)

subplot(2,1,1), plot(t(1:100),x(1:100))

xlabel('TIME, secs'); ylabel('AMPLITUDE'); title('TIME SIGNAL')

subplot(2,1,2), plot(f, I)

xlabel('FREQUENCY, Hz'); title('PERIODOGRAM ESTIMATE')

disp('area of PSD'); disp(2*sum(I)/(N*T))

Figure 4 shows a sample of the signal and the periodogram
estimate. Notice that although the area is close to 1 the shape of the PSD does
not resemble at all the constant value of
Ts^{2}
= 1. If one uses imagination, one could say that the average value is close to
1. __Calculate the average value of I(m). Is it close to 1? Estimate the PSD
for N = 2,000 and 5,000. Are they any better?__ As you can see they are not;
this is the definition of an *inconsistent*
estimator; as N increases the estimate does not get any better. Let's try an
additional procedure. First segment the original signal into 4 sample signals
and calculate the individual PSDs as shown below.

xs = reshape(x,M,K); %K averages, M points per signal

for i=1:K; xsd(:,i) = detrend(xs(:,i)); end

%

h = hanning(M);

for i=1:K; xsdw(:,i) = xsd(:,i) .* h; end

%

for i=1:K

dft = T*fft(xsdw(:,i)); S(:,i) = (abs(dft(1:M/2)) .^ 2)/(M * T*
PL);

end

f = (0:(M/2)-1)/(M*T);

figure(5)

subplot(2,2,1), plot(f, S(:,1))

xlabel('FREQUENCY, Hz'); title('1st PSD ESTIMATE')

subplot(2,2,2), plot(f, S(:,2))

xlabel('FREQUENCY, Hz'); title('2nd PSD ESTIMATE')

subplot(2,2,3), plot(f, S(:,3))

xlabel('FREQUENCY, Hz'); title('3rd PSD ESTIMATE')

subplot(2,2,4), plot(f, S(:,4))

xlabel('FREQUENCY, Hz'); title('4th PSD ESTIMATE')

The results show the random nature of the periodogram
estimate. The four estimates not only do not equal 1 but they do not resemble
each other. We know from experience that averaging reduces the variance of
numbers, so let's do that. Create an *ensemble
average* of the four estimates, that is, average the four magnitudes at the
same frequency values.

figure(6)

plot(f, SAV)

xlabel('FREQUENCY, Hz'); title('AVERAGE PERIODOGRAM ESTIMATE, K =
4')

disp('area of AVE PSD'); disp(2*sum(SAV)/(M*T))

__Does this estimate more resemble the theoretical spectrum
than that in Figure 4? Why?__

__Now repeat the ensemble averaging with K = 10. Is there
an improvement in estimation? Use some measure such as the range or variance of
magnitudes in the ensemble average to show that as K increases the estimate of
I(m) improves. What happens to the frequency spacing as the number of segments
increases?__

**Application**

__Simulate a first order moving average signal with b _{0}
= 0.5 and b_{1} = 0.5, N = 1,000 and T = 0.2. Derive the theoretical
spectrum. Estimate the periodogram with an ensemble average of 4 segments. Plot
the estimate and the theoretical PSDs together. Is the estimate close. Repeat
the estimation with 8 segments. Is the latter estimate any better? For the last
estimate plot the 95% confidence bounds. Does the theoretical PSD lie between
these bounds?__ The code below should help. 'SAV' represents the ensemble
average.

% plot bounds

df = 2*8; lchi = 28.85; uchi = 6.91;

SUP = SAV * df/uchi; SLW = SAV * df/lchi;

figure

plot(fi,SAV,fi,SUP,'--',fi,SLW,'--',f,S,':'); title('ENSEMBLE AVERAGE AND BOUNDS')

xlabel('FREQUENCY, Hz'); ylabel('PSD MAG')