In this project you will explore and compare various parametric and non-parametric techniques for power density spectrum (PDS) estimation. The goal is to approximate the PDS with non-parametric direct methods, such as the periodogram, periodogram averaging, and the non-parametric indirect Blackman-Tukey method. You will also investigate parametric all-pole modeling for PDS estimation. In the process, you will have a chance to sort through many issues involving finite length signals and the DFT.
Suggested Reading: Sections 10.6, 10.7, and 11.0 to 11.5.
You are given a signal that is the result of passing white Gaussian noise with variance one ( ) through a system .
Figure 1: Generating
The goal is to estimate the PDS of , . Recall that:
Thus estimating is equivalent to estimating . Your goal is to construct estimates of from the samples of the colored noise .
The data for this project is contained in the file pj2data.mat.
Copy this file into the local directory where you are working. This file contains two vectors. The vector y is a 512-point vector representing the discrete time signal . The vector Hejw2 is a 512-point DFT representing samples of the magnitude-squared response . It is the desired response that you are trying to estimate from . Hejw2 is given to you as a baseline for error calculation.
You can load the data into your MATLAB environment using the load command. Just cd to the directory where you copied pj2data.mat and at the prompt type:
Your MATLAB environment will now have two vectors y and Hejw2 defined in it. Warning! If you had other variables named y and Hejw2 they will be overwritten by the load command.
Figure 2: Plots of the data in pj2data.mat
As a first exercise you may want to plot the data y and Hejw2. They are plotted above for your convenience. Note that the k-indices on Hejw2 represent sampled values in the interval (i.e. ).
In this project you are using the data to construct estimates of , where .
Error Criterion: The error criterion is defined as the discrete-squared-error in the frequency domain. An expression for the estimation error is given in (1):
Autocorrelation Estimate: The autocorrelation function of a wide-sense stationary signal $y[n]$ is defined as:
In this project, the autocorrelation estimate provided by the observed data of length is:
Functions you may find useful in this project include fft(), fftshift(), fliplr(), conv(), downsample(), xcorr(), toeplitz(), levinson(), freqz(), sum(), abs(), transpose(), ones(), zeros(), triang(). Look at the appropriate help files in MATLAB to see how to use them.
A. Autocorrelation and MATLAB
The function xcorr() in MATLAB calculates the finite autocorrelation between two sequences. Calculate the following sequence
= xcorr(y, y, ‘biased’)
This is similar to the convolution of two finite sequences and . For the first part of this exercise, show this fact using the following example:
1. Take the first 32 points of (we define this signal as ).
2. Plot the autocorrelation of using xcorr().
3. Then plot the autocorrelation using the conv() function.
The two plots you get should be the same, except for a scaling factor.
What does xcorr calculate if you replace ‘biased’ with ‘unbiased’?
(a) The deterministic autocorrelation of is
where is the length of . Explain why the Fourier transform of this autocorrelation function is a positive, real function.
(b) Take the 64-point DFT of using fft command in MATLAB. Plot the absolute value and the phase of the DFT of . What is the relationship between the 64 point DFT of and the 64 point DFT of ?
(c) Use and generate a signal whose 64 point DFT, which is calculated by fft, is the same as the 64 point DFT of .
The PDS is the Fourier transform of the autocorrelation function defined in (2). The Fourier transform of is an estimate of the PDS, which is the 64-point periodogram of using a rectangular window. The periodogram can be calculated directly from without first convolving it with itself: Take a 64-point DFT of and then find its magnitude squared. This will give you a 64-point DFT representing .
For the last part of this exercise, plot three figures:
(a) The absolute value of the 64-point DFT of , the deterministic autocorrelation sequence.
(b) The magnitude squared of the 64-point DFT of .
(c) The magnitude squared of the 64-point DFT of the first 64 points of .
What is the relationship between the three signals in A.3.a, A.3.b and A.3.c ?
B. Nonparametric PDS Estimation
Estimate the PDS of using only the first 32-points of the data. Do this by taking the 64-point periodogram. Plot the 64 DFT points of the desired frequency response, , and your estimate of the PDS at on the same plot. Calculate the estimation error as given in (1).
Estimate the PDS using all 512 points. Do this by taking the 1024-point periodogram of . Similar to Problem B.1, plot the 64 DFT points of the desired frequency response, , and your estimate of the PDS at those frequencies on the same plot. Calculate the estimation error as given in (1).
Estimate the PDS using periodogram averaging. Write a script to find the 64-point periodogram of every 32 non-overlapping samples and average the results (i.e. there will be 512/32 = 16 periodograms to average). The expression for the periodogram averaging is:
Similar to Problem B.1, plot the 64 DFT points of the desired frequency response and your estimate on the same plot. Comment on any differences between this estimate and the previous two. Calculate the estimation error as given in (1).
In this problem you will use the indirect Blackman-Tukey method to estimate the PDS of . The Blackman-Tukey method requires you to first estimate the autocorrelation of . Use the following three steps to generate your Blackman-Tukey estimate:
(1) Estimate the autocorrelation of using (3) and .
(2) Truncate your autocorrelation estimate by using only for (i.e. Do not use the estimate of autocorrelation for .)
(3) Take the 64-point DFT of the truncated autocorrelation function.
Similar to Problem B.1, plot 64 DFT points of the desired frequency response and your estimate on the same plot. Calculate the estimation error as given in (1).
Format the errors you found in Problems B.1-B.4 in a table and discuss your results.
(a) Which method performed the best?
(b) Explain why the Blackman-Tukey method does so well without any averaging.
(c) Multiply your estimated autocorrelation obtained in part B.4.(2) with , a triangular window of size (nonzero for ) with unit center tap. Take a -point DFT of this new estimate. Calculate the estimation error as given in (1). Discuss why the performance is different than your result for problem (B.4).
C. Parametric PDS Estimation
Estimate the PDS of using all-pole modeling. Here the stable impulse response in Figure 1 is modeled with a class of stable filters with the following structure
The goal is to fit an all-pole model to the data and use that model to estimate .
Formulate the Yule-Walker equations for 2nd through 7th order all-pole models, and solve for the coefficients to estimate . Use the values of as estimated in Problem B.4. In a table, list the coefficients and the prediction error you found for . For which order is the prediction error minimum? (hint: check out the levinson command in MATLAB)
Draw the lattice implementation for .
In addition to finding the denominator coefficients, you must also find the gain factor . In the Yule-Walker estimation method, the minimum mean-squared value for in the th-order predictor is given by
Plot the estimation error as defined in (1) as a function of . You may find the command freqz() helpful. Comment on any interesting features of this error curve. What does this imply about the system function ? What is your best estimate of the number of poles in ? Compare the prediction error to the estimation error.
(a) Another method to find the gain factor is to use the Yule-Walker equations. In this case, the estimate of is given by:
Show how this equation is obtained. Use and find the estimation error in (1). Compare s with the estimates s in Problem C.2.
(b) Another estimator is defined based on the prediction error signal, . The estimator is
Explain how this estimator is obtained. Repeat C.3.a for this estimator.
MATLAB contains functions called
View the help files for these functions and explore their use in the context of Problems B.1, B.2, B.3 and C. Show how these commands can provide the same results which you obtained for these relevant problems . Do not use these functions to find answers in the previous problems in the project.