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 a0 = 1, a1 = -1, a2 = 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.

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.

th = ar(y,2,'burg');

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.

clear;

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  Ts2 = 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.

K = 4; M = N/K;

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.

SAV = mean(S');

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 b0 = 0.5 and b1 = 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')