AA: Power Spectral Density and Filtering

Assignment Summary

If your last (family) name begins with A-N, please complete the "Appl. Assignment 6 Quiz" on WebCT after finishing work on this assignment. If your last (family) name begins with O-Z, please complete the "Appl. Assignment 6 Quiz" on WebCT PRIOR TO STARTING work on this assignment. The quiz is found in the "Assessments" tab of the WebCT page.

Last Name Begins With When To Take the Quiz
A-N After completing the assignment
O-Z BEFORE STARTING the assignment

In Application Assignment 5, we computed the autocorrelation function for the temporal fading signal. In this assignment, we will compute the power spectral density from the autocorrelation function. Then, we will filter the temporal fading signal, and compute the power spectral density of the output.

Motivation: In a secret key generation system, the high-frequency components of the temporal fading signal are not as correlated across the two directions of the link. They contain quite a bit of noise, and relatively little signal power. So filtering with a low pass filter the measured fading signal at both ends of the link is useful in the reduction of the noise. We will show that, using analysis of power spectral density and the effect of LTI filters, we can design a filter to reduce the high-frequency components of the filter in a predictable manner.

There are three parts to this assignment:

  1. Compute and plot the power spectral density of the temporal fading signal. Turn in: A plot of the PSD.
  2. Filter the temporal fading signal with a low-pass filter. Turn in: A plot of the frequency response of your filter.
  3. Compute and plot the power spectral density of the output of the filter. Turn in: A plot of the PSD.

You will need Matlab, a text file from an RSS measurement experiment, and some way to produce a document (eg., Word, OpenOffice). Create one PDF containing the outputs of your assignment, and the Matlab code you used to generate them. Turn it into the WebCT assignment drop box, by midnight on the due date.

Assignment Description

Load the data as you did in Application Assignment 2. This involves loading, separating, and interpolating the data so that you have data vectors "interpa" and "interpb" in memory.

1. Compute the PSD of the fading signal

Compute the PSD of either the interpa or the interpb fading signal. There are two ways to compute the PSD:

  1. Direct method: The PSD is effectively what we'd see if we set a spectrum analyzer on average. The definition of the PSD (see Def'n 11.4 on page 419 of Y&G) is a direct expression of this intuition. So, we can compute the PSD by computing the FFT on lots of subsequent "windows" of data and averaging them together.
  2. Compute the Fourier Transform of the autocorrelation function: The Wiener-Khintchine theorem says that the autocorrelation function $R_X[k]$ and PSD $S_X(f)$ are a Fourier Transform pair.

You should try both methods -- the process should convince you that the analysis actually works.
Turn in a plot of the PSD (using one of the two described methods).

1.1 Direct Method:

The PSD is defined (Def'n 11.4 on page 419 of Y&G) as
$$ \lim_{L\rightarrow \infty} \frac{1}{2L+1} E\{|X_L(f)|^2\} $$
where $X_L(f)$ is the $(2L+1)$-length Discrete-Time Fourier Transform of the random sequence $ X_n $.

To compute an estimate of the PSD in interpa from the direct method, I used this code:

len = length(interpa);
windowLen = 400; % How long of an FFT to take each time.
trials = len - windowLen - 1; % How many different windowed FFTs can be taken
totalPS = zeros(1,windowLen); % Total of power spectrum
for i = 1:trials,
% Compute the magnitude squared of the FFT of part of interpa
windowPS = abs(fft(interpa((1+i):(i+windowLen)))).^2;
totalPS = totalPS + windowPS; % Sum for later average.
PSD = (1 / windowLen) .* (totalPS./trials);

Note this last line is normalized by 1/windowLen. Our "windowLen" is equal to the $2L+1$ in the above equation.

1.2 Compute the Fourier Transform of the autocorrelation function:

This code assumes that you've computed the autocorrelation function as described in Application Assignment 5. In that assignment, you stored in autocorr computed values of the autocorrelation function. In particular autocorr(1) stored the value of $R_X(0)$, autocorr(2) stored the value of $R_X(1)$, and so on. The shift by one was required since Matlab won't accept a vector index of 0. For purposes of the Fourier transform, we actually need to take the Fourier Transform of the vector of $ [R_X(-L), \ldots, R_X(L)] $. We did not compute $R_X(\cdot)$ for negative time delays! But, it is an even function: $R_X(-L) = R_X(L)$ for all $L$. We fix this in one line of Matlab code by rearranging autocorr to repeat itself:

symmetricAutoCorr = [autocorr(1:200) fliplr(autocorr(2:199))];

Then, applying the Wiener-Khintchine theorem:

PSD = fft(symmetricAutoCorr);

This method requires only these two lines of code. If you've done it correctly, the PSD will be purely real; if not, there is something wrong with your symmetricAutoCorr.

1.3 Plotting the PSD:

In either case, plot the PSD (in dB) with the commands:

f = 0:2/length(PSD):2-2/length(PSD);
plot(f, 10.*log10(PSD));
set(gca,'xlim',[0 1])
xlabel('Normalized Frequency (\times\pi rad/sec)');
ylabel('Power Spectral Density (dB)');

2. Design a low-pass filter.

Design a low pass filter to attenuate power in frequencies higher than about 0.2 (normalized frequency). The exact design is up to you. If you want, for example, get help on Matlab's "butter" or "cheby1" to get information on some easy Matlab functions for standard filter designs.

Use Matlab's freqz command to plot the frequency response (in dB) of your designed filter. Note that the upper subplot of the freqz plot shows $20 \log_{10} |H(f)|$, that is, the power gain vs. frequency for the filter. Turn in this plot.

3. Filter the fading signal and plot the PSD of the filter output

Use Matlab's filter command to filter the fading signal (either interpa or interpb). Let's call the output filteredOut. Compute the PSD in the output signal. You may use either method to compute the PSD, or both, if you want to check your work. Turn in one plot of the PSD of the output signal.

We can check the output PSD by looking at the input signal PSD and the frequency response. We know that
$$ S_Y(f) = |H(f)|^2 S_X(f)$$
where $|H(f)|^2$ is the squared magnitude of the filter frequency response, $S_Y(f)$ is the PSD of the output of the filter. On all of these plots, we are looking on a (dB) scale, that is $10\log_{10}$ of the linear value. So the above expression becomes:
$$ 10\log_{10}S_Y(f) = 20\log_{10}|H(f)| + 10\log_{10} S_X(f)$$
You should be able to see that the shape of the filter power gain in dB has added in to the PSD of the original signal. However, note that it is difficult to come up with a PSD in (dB) in Matlab which is very very low. This is because the average is done in linear terms, and whatever "noise" which exists in each FFT is still positive; it is difficult to average this to close to zero when the noise always contributes a positive amount.


Submit your results as described above.

If your last (family) name begins with A-N, please complete the "Appl. Assignment 6 Quiz" on WebCT now.