1.0 Objectives:

1. To write a Matlab script-file (program) to compute an energy density spectrum using the Matlab FFT function fft(x).

 
2. To test this program by generating synthetic data composed of a sinusoid plus Gaussian random noise (white noise).

2.0 Procedure:

1. Read the background information on the Fast Fourier Transform.

 
2. Before you write the matlab program, be sure you know what the spectrum of a sine wave plus noise would look like.  Use the sample spectral diagrams in your class notes to draw a sketch of the expected spectrum. Write down the formula for the variance of a sine wave plus noise.

 
3. Use a text editor to create your matlab script file.  Remember, your matlab program name must end in a .m.

Be sure to document your program.  The first few lines should include the program name, your name, the date, and a brief description of what this program does. Intersperse comments in appropriate places throughout the program.  Remember, comments are preceded with a %
 

4. Type in the following - using the values you have decided upon.


fs =     ;             % fs is the sampling frequency. Try 1,2, or 4

dt =     ;             % dt is the time between sampled points (inverse of fs)

N =     ;             % N is the number of sampled points. Try 128 or 256.

NS =     ;         % NS is the number of samples (for Lab 4, set NS=1)

T = N*dt;         % T is the sampling period

df = 1/T;         % df is the fundamental frequency

Nf= N/2 + 1;     % Nf is the index of the Nyquist frequency

fn = ??             % the Nyquist frequency - what is it?
 

5. Generate a time vector, t, which starts at dt and ends at T. Then create your synthetic data, ydata,  using the matlab functions cos(x) and randn(x). Your synthetic data will be made up of a sinusoid, y1, combined with some random white noise, y2. Your code should look something like this:


y1=A*cos(2*pi*f*t - phase) + off;             % generate a sinusoid vector

y2 = randn(size(y1));             % generate Gaussian noise vector, same size as y1
                                               % (note that it won't be perfectly Gaussian, since
                                               %  y1 contains only a finite number of points)

y2= y2*sd + xmean;             % Set standard deviation and offset of noise.

ydata=y1 + y2;                     % this is your synthetic data!

In creating y1, you'll also need to select values for off, the mean value or offset, A, the amplitude, f, the sinusoid's frequency, and phase, the angle in radians. Be careful what value you use for f so you don't get leakage or aliasing. Before generating y2, assign values to the constants sd and xmean. The values you choose will become the standard deviation and offset, respectively, of the random y2 data.

figure(1)                                % open 1st graphics window

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Stop here, run your m-file and check that the synthetic data is what you expect.
How would you do this?
Take advantage of matlab's interactive features by checking your code frequently.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

6. Using the reshape function, convert ydata into a matrix with N rows and NS columns.

 
7. Demean your data.

 
8. Do the FFT of your synthetic data and calculate the energy density spectrum. Then check the computed amplitude and variance with your expected amplitude and variance.  Here are some of the commands you'll need:


X = fft(ydata);                 % two-sided matlab amplitude spectrum (see backgrd info sheet)

X = X/N;                        % convert to our form of amplitude spectrum
                                              % (see backgrd info sheet - discrete equations)

X(Nf+1:N, :) = [];            % Create single-sided spectrum by folding (i.e., deleting freq>fn
Ax = 2*X(2:Nf, :);          % and doubling the magnitude of remaining components; leave out freq=0)

When matlab performs an FFT, it returns a double-sided spectrum of N amplitude values corresponding to frequencies f₁ (= 0*df = 0) through f_N (= (N-1)*df). Matlab's amplitude values are symmetric about fn, the Nyquist frequency, rather than zero frequency as you've seen in class. The above folding calculation will result in a single-sided spectrum running from df to fn.
 

9. Next, convert your amplitude spectrum into a mean energy density spectrum.

   Gx = abs(Ax) .^2/(2*df);         % the single-sided energy density spectrum from df to fn (Note `.^')
   Gxmean=mean(Gx,2);             % Gxmean becomes a column vector (average over the ensembles)
    

10. Plot Gxmean vs freq and the log of Gxmean vs freq on two panels on the same page using subplot.  (Use matlab's log10 function to calculate the log of Gxmean, or use semilogy.)  Label and title your plots using ylabel, xlabel, and title.

 
11. Using max and find, obtain the magnitude and index of the peak amplitude. Also find the peak frequency.  Using Parseval's Theorem (and matlab function sum), calculate the variance.   Don't forget to generate a vector called freq with the appropriate range of frequencies corresponding to the Gxmean values.

 
12. Print a copy of your spectral plots.

 
13. Save values of variance, peak-amplitude, df, and peak-frequency, as well as a few key values of freq vs Gxmean.  You'll also want to remember off, A, phase, sd, and xmean from your synthetic data so you can compute the theoretical variance and peak amplitude.

3.0 Report:

1. Turn in a copy of your Matlab program (Script file).

 
2. Turn in the page with plots of your Energy Density Spectrum.

 
3. Turn in a text file containing the computed variance, peak amplitude, peak frequency, important values of freq vs Gxmean (including at least the values at and around the peak, and the lowest/highest/Nyquist/average noise), and the input parameters for your synthetic data.

 
4. Modify the sketch you did of the expected spectrum so that the spectral peak is discrete, of finite height and width (i.e., more like your Matlab plot). Label important frequencies and amplitudes along with their estimated values, including (at least) peak frequency and amplitude, Nyquist frequency, and average amplitude of the noise. You'll need to actually calculate estimated values for the heights of the amplitude peak and the noise. Remember that Matlab must represent an infinite spike with zero width as a finite triangle with base = 2*df, because df > 0. Hand in this modified sketch with your report. How well do your estimates for the spectral peak and the average noise height compare to the values determined by Matlab? Why might your program's values be different from your estimates? Be specific.

 
5. Calculate the theoretical variance for your synthetic data (show work). Compare this to the value obtained from integrating the spectrum generated by your program, and try to explain any difference.