Though we discussed in class how the FFT (an algorithm for finding the discrete Fourier transform) is carried out, you need not write your own code to implement this algorithm. This is because it is already implemented in Matlab. For an n-vector y, fft(y) is its n-vector DFT. The setting for what follows is that we will think of y as containing n sampled values (y-coordinates measured at evenly-space t-coordinates) from some continuous function of t. The spacing between these t-values is dictated by the sampling rate (expressed in Hertz, or cycles per second).
Along with this common setting, there will be a common list of Matlab commands that we use in some way, shape or form. Here they are with explanation, assuming that we already have a y vector to work with.
Set the sampling frequency for y: >> sf = value; Play y as a sound: >> sound(y, sf) The FFT of y (entries each associated with a specific frequency): >> Y = fft(y); Get the strength of frequencies in first half of DFT: >> pow = abs(y(1:n/2+1)); The frequencies with which the strengths in pow are associated: >> f = (0:n/2)*sf/n; Plot strengths of frequencies vs. frequency: >> plot(f, pow)
An example from class
Let's sample the sine function whose frequency is 440 Hz. We choose
a sampling frequency that is more than double this, like 1024
(= 210) Hz, and sample for 5 seconds.
>> sf = 1024;
>> n = 5*sf;
>> t = (0:n-1)/sf;
>> y = sin(440*2*pi*t);
>> sound(y, sf)
>> Y = fft(y);
>> pow = abs(y(1:n/2+1));
>> f = (0:n/2)*sf/n;
>> plot(f, pow)
Aliasing
Why do we need to sample at a rate at least twice that of the
highest frequency in the signal? Consider the sine curve
y = sin(10*pi*t), whose frequency is 5 Hz. Let's first sample
it for one second at a frequency of 5 Hz (i.e., 5 times each
second).
>> sf = 5;These commands show you what has been sampled, not its DFT. You can see that the vector of sampled values is (except for roundoff errors) full of zeros. If you overlay the actual sine curve you have sampled, it is clear why this is so:
>> n = 1*sf;
>> t = (0:n-1)/sf;
>> y = sin(5*2*pi*t);
>> plot(t, y, 'o')
>> axis([0 1 -1 1])
>> moret = (0:1000)/1000;If you stop and think for a second, you'll realize that there are lower-frequency sine curves that pass through these very same points. For instance, the sine curve whose frequency is 2.5 Hz:
>> plot(t, y, 'o', moret, sin(5*2*pi*moret))
>> axis([0 1 -1 1])
>> plot(t, y, 'o', moret, sin(5*2*pi*moret), moret, sin(2.5*2*pi*moret))Simply raising the sampling frequency a little will not fix the problem. Let's try a sampling frequency of 8 Hz:
>> sf = 8;This plot shows the 8 sampled values (for 1 second), the 5-Hertz function being sampled, and another function, the sum of 4 and 3-Hertz sine curves, passing through the same points. Now look at the DFT of y:
>> n = 1*sf;
>> t = (0:n-1)/sf;
>> y = sin(5*2*pi*t);
>> plot(t, y, 'o', moret, sin(5*2*pi*moret), moret, ...
sin(4*2*pi*moret)-sin(3*2*pi*moret))
>> Y = fft(y)(no semicolons this time). The 4th and 5th coefficients (the ones corresponding to 3 and 4 Hz) are the only nonzero ones. When one frequency is mistaken for a lower frequency (or for a linear combination of lower frequencies), this is called aliasing. (You can be sure that the sound corresponding to the red curve is quite different than that for the green one.) One can only avoid aliasing by sampling at a high enough frequency—more than twice the highest frequency present in the sampled function.
Sunspots and el Nino
Matlab comes with a data file called sunspot.dat.
Cleve B. Moler, one of the founders of The Mathworks, the company
that makes Matlab, has a book the chapters of which are all
available for free download at
http://www.mathworks.com/moler/chapters.html.
Download
Chapter 8 and work through Section 4, which carries you
through an analysis of this data. Then download the file
elnino.dat (Place it in the “scofield” subdirectory)
and do Exercise 8.7, found at the end of Moler's Fourier Analysis
chapter.
This page maintained by:
Thomas L. Scofield
Department of Mathematics and Statistics,
Calvin College