Sequence Analysis II: Optimal Filtering and Spectral Analysis

File last modified 9 October 1996

7.3 Power Spectral Analysis

7.3.1 Types of Spectra

You need to be careful when you look at power spectra. There are a wide variety of ways that have been used to compute the total power. Three different schema that are widely used are (1) "sum squared amplitude", given by

(2) the mean squared amplitude, given by

and (3) the time-integral squared amplitude, given by

where "c_j" are the coefficients, N is the length of the series, and is the sampling interval. However you define the total power, something called Parseval's Theorem must apply. It states that the total power computed in the time domain must equal the total power computed in the frequency domain:

More confused yet is the Power Spectral Density (PSD) which is what you actually plot when you plot a spectrum (you plot the PSD vs. frequency). For any given definition of the total power, you could define the PSD to range over all positive and negative frequencies, between -f_c and f_c where "f_c" is the critical or nyquist frequency (we'll learn more about that in section 7.4), or over only the positive range (either all frequencies or up to f_c). We'll tell you what MATLAB does in a moment.

What does the power spectrum look like for different functions? For example, the spectrum of a pure sine wave (that goes on forever) of a fixed frequency would be a spike at that frequency. This would be referred to as "monochromatic". The spectrum of random (gaussian) noise would be a flat line (with some bumps and wiggles, but none significant). This is referred to as "white". A spectrum dominated by a lot of energy in the low frequency range is referred to as "red", and a spectrum dominated by high frequencies is called "blue".

7.3.2 Doing it with FFT

An important step in doing a spectrum, which is sometimes ignored, is to "pre-whiten" the spectrum. This is usually simply the task of removing any long term trends from the data, which will tend to redden the spectrum. Another possibility is that the signal of interest is fundamentally bandwidth limited (for example, if you are looking for seasonal signals in a daily record of temperature. You are not interested in the high frequency (or blue) end of the spectrum so you low pass filter or smooth the data. Having done this, you proceed with the next few steps.

We show you this procedure to give you a rough idea of what book-keeping is involved in the process, and so you will understand "what is going on under the hood". We'll use a simple, artificial example so that you can see what is going on. Feel free to play along with us at home…

First, let's create our little data set. We'll create a set sets of time "t" and create a time-series "y" which is the sum of two separate frequency sine waves of different amplitudes, plus some random noise. In MATLAB, this is done with

deltat = .001 ;    % size of time step

t = 0:deltat:1 ;   % NB has 1001 elements

f1 = 50; f2 = 120;  % our two frequencies

y = sin(2*pi*f1*t) + 2*sin(2*pi*f2*t) + 0.5*randn(size(t)); %
the time series

We now do a n=1024 point FFT (remember that FFTs like radix 2 numbers, and 1024 is the closest radix 2 number to 1001. All MATLAB does is to pad the time series out with zeros). When MATLAB does the FFT, it will create values of the transform "Y" which range from -f_c to +f_c with n+1 distinct values. The manner in which they are stored is that the series in "Y" starts at f=0, goes to f_c at the (n/2+1)^th element, and puts the negative frequencies in the latter half of the matirix. Hence for the 1024 point FFT, the f_c coefficient is stored in Y(513).

n = 1024;   % FFT length (radix 2 number!)

m = n/2;   % number of distinct frequency bins

Y = fft(y,n);   % FFT of length 1024

Pyy = Y.*conj(Y) / n;  % PSD = |Y|^2

fc = 1/(2*deltat);  % the critical frequency

f = fc * [0:m]/m ;   % the frequency bins (513).

Note that in the fourth line above, we calculated the PSD by multiplying the Y elements times their complex conjugates (note the little dot too!). This is because the fourier transform values are in general complex quantities - that is, they contain both real and imaginary components. The square-amplitude of a complex number is always obtained by multiplying itself by its complex conjugate. Further, we normalize the PSD by dividing by the series length. Finally, we define the critical frequency (more about that later) "f_c", which is the maximum observable frequency for the series, and we create a vector containing the frequencies of the PSDs.

We now do a bit of book-keeping. First, the sequence is always symmetric about the zero frequency component. We throw away the top half of the sequence, and double the remaining frequencies, except for the zero frequency component Pyy(1).

Pyy(m+2:n) = [ ] ;  % sneaky MATLAB trick for shortening a vector

Pyy(2:m+1) = 2*Pyy(2:m+1);  % compensate for missing negative
frequencies

semilogy(f,Pyy);

Note that the two frequencies come through loud and clear. You can see them at about 50 and 120 (hard to tell on this plot, but you can easily devise a way of finding the frequencies of the maxima with MATLAB. We'll leave that one up to you. Also difficult to tell from this plot, but something you can check with MATLAB is the size of the maxima. Oh heck, we may as well tell you:

[Peaks IFreqs] = sort(-Pyy) ;

The trick here is that the MATLAB sort routine sorts a vector in ascending order, so we negate the vector to sort it in the opposite direction. If you specify two outputs for sort, it also returns the indices of the sorted matrix, which also tells you the frequencies. We now look at the 5 largest values, and their frequencies with :

abs(Peaks(1:7))

f(IFreqs(1:7))

We'll pretty up the output here (doesn't MATLAB type some ugly numbers sometimes?), and don't forget that because we used a random number generator in this example, your mileage may vary…

abs(Peaks(1:7)) = [1898.6    425.1    57.1    32.    15.8   12.7
  11.4]

f(IFreqs(1:7)) = [120.1   49.8   119.1   50.8   118.1   121.1
  48.8 ]

Notice a few things. The biggest numbers are near the peaks, and spot on the frequencies we put into the data set. Note also that the next biggest peaks are nearest those frequencies. This is very interesting, because it says that the sine waves don't produce perfect spikes. There is a reason for this, and we will be discussing this in section 4. Note that the difference in PSD between the two peaks is not exactly 4 (since the amplitudes differ by exactly 2, the PSDs should differ by 4). This is related to the problem talked about above (the areas under the peaks are more important), and also to the presence of noise in the data. More on this in section 4.

7.3.3 Frequency Binning, and Using Specplot

Well, now you've done it the hard (and hopefully instructive) way, we'll show you the easy way. But first, we need to ask a general question about how well you know the PSD. Remember our basic philosophy: it is not enough to determine a quantity, you must also determine the uncertainty with which you know that quantity. Specifically, what is the variance in the PSD estimate as N (thte length of the time series) goes to infinity? The answer seems terribly unfair: the variance in PSD remains the same, independent of the length of the time series: it is 100% of the value. The reason for this is actually quite straight forward. Using the FFT technique described above, if we double the length of the time series, we halve the width of the corresponding frequency bins, and hence we double the number of frequency bins. All our efforts go into refining the frequency resolution.

We don't (and shouldn't) have to settle for 100% uncertainty in our PSD estimates. What we can do is to average frequency bins to gain precision at the expense of resolution. This averaging has the effect of reducing the variance by the square root of the number of bins involved in the averaging, so that averaging 100 bins reduces the error to 10%. An equivalent procedure (and one generally preferred for logistical reasons on longer data sets) is to break the time series up into "n" non-overlapping segments, computing the spectra for each of those segments, and averaging them together. This is exactly the same thing as doing one big spectrum and averaging the bins in groups of "n".

To do all of the above in MATLAB, we take our little series, and do the following:

P=spectrum(y,128);

specplot(P,1000);

The first statement computes the spectrum by breaking the original series into lengths of 128 elements each (because it uses FFT, the lengths of the sequences must be radix-2). This produces several estimates at each frequency, thus improving the uncertainty in each estimate. Note that specplot is smart enough to know that even though you specified the sampling frequency (i.e., 1/sample_interval), that the ultimate resolution is the nyquist frequency which is half the sampling frequency. Specplot gives you the following:

The yellow line is the estimated PSD, while the lines above and below delimit the 95% confidence interval for the PSD estimates. Note that the trade-off between precision and resolution is evident in this plot, because the peaks, although smoother, are now wider than the straight FFT based estimate we did earlier. We have an exercise for the student: try generating the time series sequence as we have done above, and try doing the specplot using a variety of segment lengths: say 32,64,128,256 and 512. Observe what happens to both the width of the spectral peaks and the confidence intervals. Be warned that the y-axis will change from plot to plot. You can look at the numbers themselves, since the "P" that the MATLAB "spectrum" routine returns is both the PSD (as column 1) and the confidence range in the PSD (as column 2) in P.

GoTo Next Section
GoTo Lecture TOC
GoTo 12.747 TOC