Description

This image depicts the overlay of hundreds of FFT plots of heart-rate samples, taken 1 second at a time. The spikes indicate the frequency components of the heart rates signal, where the first harmonic (at around the 1000th index) stores information about the average heart rate in beats/minute. Our ultimate goal in this lab will be to apply MATLAB’s FFT function to analyze the average heart rate of provided ECG signals.
Part 1
The goal of this lab is to understand the use of MATLAB’s implementation of the Fast Fourier Transform (FFT) and to apply it in order to find the heart rate in some real-world electrocardiogram (ECG) signals.
In this pre-lab, you will be introduced to the concept of the FFT. A basic definition is provided, as well as some general background and terminology that will be used throughout this lab.
1.1 Definition: What is the FFT?
FFT is an acronym for the Fast Fourier Transform, an algorithm that computes the discrete Fourier transform (DFT) of a discrete-time signal. In other words, the FFT takes a discrete-time signal and extracts important frequency domain information from that signal. It is considered “fast” because its computational complexity is far less than that of the DFT, the Discrete Fourier Transform (see https://en.wikipedia.org/wiki/Fast_Fourier_transform). In this lab, we won’t care about the mechanics of the algorithm (hurrah!); we are simply concerned with how to use the FFT function in MATLAB to determine frequency domain information of signals. IMPORTANT: The FFT function performs the DFT of a signal; the FFT is specific to an input signal with length 2^n, but the FFT function in MATLAB can be used for input signals that do not meet this criteria: this just performs a DFT (see help fft for more information). The two terms are used interchangeably throughout the lab.
1.2 Background: How does the FFT “Read” Frequency Information?
“matches” these basis signals, and outputs a vector of complex numbers – one for each iteration of the algorithm. In this lab, we will call these complex measurements frequency bins, or simply bins.
Because the outputs of the FFT are complex, we need an easier way to visualize them. Simply plotting the vector of frequency bins will result in a
complicated jumble of vectors in the complex plane. Therefore, we “read” the results of the FFT by plotting the magnitude spectrum of the frequency bins. The higher the magnitude of the frequency bin, the higher the correlation between the input signal and a sinusoid of that frequency. Therefore, the magnitude spectrum should spike at regions where the frequency of the sinusoidal basis signals match the input signal. Recall that any signal can be expressed as the sum of sinusoids of different frequencies. Consequently, the FFT is an extremely powerful tool because it can help us decompose a complicated discrete-time signal into its sinusoidal components.

1.3 Exercise: Decomposing Signals into Sinusoids
Rewrite the following signals, by hand, as a sum of sinusoids (that is, expand them as much as possible). Identify the amplitude, frequency, and phase of each sinusoidal component. In your write-up, provide your answers with explanations or scans of your work. (Note: this should be a trivial task that serves essentially as a review of the beginning of the course).

1.4 FFT of a Single Sinusoid with Integer Frequency
In this section, we will use MATLAB’s FFT function to analyze basic signals – ones that can be expressed mathematically. We will then apply what we know from this warmup to apply the FFT to a real-world example, where the frequency components of the signals (ECG heart rate data) are unknown.
To start, we’ll use the FFT function for a discrete-time signal that is a single sinusoid. This way it will be clearer as to what the FFT is really doing. Let’s consider the sinusoid defined below:

Plot the continuous-time frequency spectrum of this signal by hand.
We can see that the frequency is 50 Hz and the phase is pi/4 radians. When we run the FFT algorithm, we should expect to see the plot of the magnitude spike at the index corresponding to the frequency 50 Hz. There should only be two peaks, since the signal is composed of a single sinusoid (remember Euler’s formula for decomposing a cosine).
1. Define the variable fs for the sampling frequency, and assign it 1000 samples/second.
2. Construct a time vector using this sampling frequency that is exactly 1500 samples long, starting at zero.
3.
Define x(t) , the sinusoid described above, as the vector x .
4.
Run the fft() function on x (see help fft ).
5.
You should notice that the output is “1×1500 complex double.”. This is again because the outputs of the FFT are complex numbers. The output, in the frequency domain, is divided into a number of complex amplitudes in frequency “bins”; the number of frequency bins is equal to the number of samples. In this case, there are 1500 complex amplitudes, in 1500 frequency bins, because there are 1500 samples in our signal. To make sense of the outputs, we need to take the magnitudes and phases of the complex amplitudes in this complex vector. a
Create a new vector that stores the values of abs(X) and plot it. This plot is called the magnitude spectrum, where the xvalues are numbered frequency bins and the y-values are their corresponding magnitudes.
b Compare this plot to your continuous-time spectrum. Question 1: What’s different about it?
c
In the plot of the magnitude spectrum, you should notice two congruent peaks on either side. For the signals we will be dealing with in this lab (real signals), the magnitude spectrum will be two sided. That is, any peak will be a mirror image of the peak on the other side. Because of this redundancy, we can simply evaluate the peak on the left.
Question 2: What is the index corresponding to the left peak?
Hint: use the find() function (see help find ).
Question 3: What is the value of the magnitude spectrum everywhere else?
d
The result of the Fourier transform of a sampled signal goes into frequency bins that correspond to the normalized radian frequency, . Our convention thus far has been to draw magnitude spectra from -pi to pi. The FFT function, however, returns frequency bins ranging from 0 to 2*pi.
Plot FFT output against normalized radian frequency. Here’s the command for the normalized radian frequency output: ww = 0:(2*pi/length(XX)):(2*pi-1/length(XX)); plot(ww,abs(XX));
e
The fftshift() function (see help fftshift ) is a way to “fix” the output of the FFT function in order to make it range from −pi to pi.
Plot the shifted FFT output against normalized radian frequency: ww = -pi:(2*pi/length(XX)):(pi-1/length(XX)); plot(ww,abs(fftshift(XX))); f
The normalized radian frequency is related to the frequency in Hz by the sampling frequency and a factor of 2*pi.
Plot the FFT of X against the Hertz frequencies of the bins.
Question 4: What are the frequencies matching the peaks?
6.
The FFT can also be used to extract phase information from the signal. From the complex amplitude of the FFT function, you can use the angle() function (see help angle ) to find the phase in each bin. Plot the phase spectrum of the signal against Hertz frequency. Question 5: What are the phases of the peaks?
Question 6: Does this match what you would expect, and why?
7.
In parts [1.4.5] and [1.4.6], we identified the frequency and phase of x(t). It turns out that we can also use the FFT to determine the exact magnitude of x(t), by finding the value of the magnitude spectrum in that frequency bin.
Determine the exact magnitude of the peak. Note that the two peaks will have the same magnitude, being complex conjugates of each other. Recall that the magnitude of the peak is directly related to the length of the DFT.

1.5 FFT of Multiple Sinusoids with Integer Frequency
Now we will use the FFT function to analyze more complicated signals. In Part 1.3, we were given the functions:

These are clearly not sinusoids, but can be decomposed into a sum of finitely many sinusoids. The FFT function will be useful in determining whether your work from section 1.3 was correct. Complete the following steps:
1.
Using the same sampling frequency (1000 samples/s) and duration (1500 samples) from part 1.5, create the vectors x1 and x2 , storing

1500 samples of the signals x1(t) and x2(t).
.
2.
Plot the magnitude spectra of x1 and x2 against Hertzian frequency. Unlike in 1.5, we should observe multiple peaks, each corresponding to a single sinusoidal component. Comment on the appearance of each spectrum.
Question 7: Do you notice any differences in the number of peaks for these signals?
3.
Using MATLAB, find the frequencies, amplitudes, and phases of the sinusoidal components of x1 and x2. Compare these results to what you found by hand.

Part 2
2.1 FFT of a Sinusoid with Non-Integer Frequency
In the previous sections, we dealt with the FFT of sinusoids and signals composed of sinusoids with integer frequencies. In other words, these sinusoids exhibit an integer number of cycles per period. In this section, we will be analyzing the magnitude spectrum of the sinusoid:

Notice that this sinusoid is the same as that of an earlier section, only the frequency in Hz has changed from 50 Hz to 50.5 Hz. We want to see how the appearance of the magnitude spectrum is influenced by taking the FFT of a sinusoid with non-integer frequency.
1. Using 1000 samples/second for fs and a time vector of 1.5 seconds, construct the vector x to represent the given sinusoid.
2.
Use fft() and assign the output to X . Plot the magnitude spectrum. Zoom on the left-most peak. You should notice that the spectrum is nonzero at multiple values surrounding this peak. By introducing a non-integer frequency, we see a spread of energy across the spectrum. Find the frequency bins where this maximum occurs.
Question 1: Does this exactly match the frequency 50.5 Hz?

2.2 Zero Padding
Zero padding works by concatenating a large sample of zeros (using zeros() ) at the end of the signal being analyzed. For example, suppose our signal is

. The following script uses vector notation to add 1,500 zeros to the end of x . Here the sampling frequency is 1000 samples/second and t is a time vector of 1500 samples. The results are plotted below:

Recall that the DFT works by correlating the input signal with sinusoidal basis functions (i.e. sines and cosines with an integer numbers of cycles over the duration of the input signal). Zero padding increases the likelihood of the DFT identifying a sinusoidal basis function that perfectly matches a sinusoidal component of the input signal.
1.
Append 15,000 zeros to our signal from section 1.7 (the one with frequency 50.5 Hz), and run the FFT again.
Question 2: What is the effect of the zero padding on the magnitude spectrum of the signal?
Question 3: What is the frequency of the bin for which the magnitude is maximized?
2.
“lobe” and several smaller “side lobes” on the sides.
Question 4: How does the width of these side lobes compare to the main lobe?
3.
Try zero padding the signal with 25,000 zeros, 50,000 zeros, and 100,000 zeros. Using the subplot() function, compare the appearances of the magnitude spectra.
Question 5: What is the effect of increasing the number of zeros in the accuracy of the frequency calculation?
Question 6: Does the error increase or decrease?
Question 7: In the case of our signal, how many zeros would need to be added in order to obtain an exact value for the frequency of our input signal?

2.3 Windowing
xx = 3*cos(2*pi*2.5*t+ pi/4); xx = [xx zeros(1,15000)];

Suppose we have the signal
which we padded with zeros in the previous section. If we multiply our signal with a bell-shaped function (the window) to reduce the amplitude of the endpoints of the function, we obtain a signal that fades in and out. In this lab we will use the Hann Window, similar to the Hamming window from previous labs, only differing by a constant. The following code implements the Hann Window:
xx = 3*cos(2*pi*2.5*t+ pi/4); xx = xx.*hanning(length(xx))’;
The results of implementing this code snippet are shown above.

1.
Returning to our signal
from the past two sections, we want to reduce the size of the “side lobes” in the plot of the magnitude spectrum. Using the script above, window x(t) with a Hanning window.
Note: this needs to be done before zero padding! As before, use a sampling frequency of 1000 samples/second and a time vector of 1.5 seconds.

2.
Padding with an appropriately large number of zeros, plot the magnitude spectrum of x(t) after windowing and zero padding. Make a sideby-side plot of the magnitude spectra, zoomed in at the peak value, for X with windowing and X without windowing.
Question 8: What effect does windowing have on the FFT?
Question 9: How does the magnitude of the maximum peak change with windowing?
Question 10: Is it higher or lower than that of the maximum peak without windowing? Question 11: Has the width changed at all?

2.4 Determining Average Heart Rate Using FFT

In this exercise, you will employ the FFT function to automatically extract the heart rate from electrocardiogram (ECG) measurements. The data file on Canvas, heartratelabdata, contains four ECG signals (obtained from http://physionet.org/cgi- bin/atm/ATM). The sampling rate for all data files is 125 Hz. Your goal will be to use the Fourier transform to produce four plots, showing the average heart rate changing over each signal’s duration.
1.
Load the ECG data. There will be four arrays – am101 , am105 , am107 , and cm110 – each containing an hour of ECG data.
2.
In this step, you will program a function or write a MATLAB script to read each signal and determine the average heart rate. Since each signal is one hour long, you will need to divide the signal into subsamples (or ‘windows’) of equal width. It is suggested you start by using a window about 2 minutes long. DO THE MATH in order to condense it to 2 minutes with even windows.
3.
4.
Take the absolute value of the FFT of the sample and identify the peak closest to bin zero, but not at zero. You are looking for the first harmonic. Your algorithm should automatically identify the first harmonic each time the function iterates, because it is this peak that stores the frequency information of the heart rate.
5.
6.
Store these average heart rates into a vector, and then produce a plot showing the change in the average heart rate over time.
7.
Write the above into a function that takes a data vector ( am101 , am105 , am107 or cm110 ) as input, calculates the heart rate for each two-minute window, and plots heart rate over time.
8.
Use that function to produce the four plots, one for each signal.