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:

- Compute and plot the power spectral density of the temporal fading signal.
*Turn in*: A plot of the PSD. - Filter the temporal fading signal with a low-pass filter.
*Turn in*: A plot of the frequency response of your filter. - 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.

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.

Compute the PSD of either the `interpa`

or the `interpb`

fading signal. There are two ways to compute the PSD:

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

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.

end

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.

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`

.

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,'FontSize',20);

set(gca,'xlim',[0 1])

xlabel('Normalized Frequency (\times\pi rad/sec)');

ylabel('Power Spectral Density (dB)');

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.

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.