Description of this paper

Loading

Hi, Im wondering if anyone knows how to do the lab exercise worth-(Answered)

Description

Instant Solution ? Click "Buy button" to Download the solution File


Question

Hi, Im wondering if anyone knows how to do the lab exercise worth 50 points


Lab 10: An Introduction to Using the FFT

 

Function in MATLAB

 


 

Figure 1 ? 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.

 


 

1. Pre-Lab (10 pts.)

 

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. It is most

 

important that you watch the lecture videos that are provided in section 1.3. At the end of this prelab there is a small exercise, where you will rewrite given signals as a sum of sinusoids (basically

 

a review of the material at the beginning of the course).

 


 

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

 

complexity is less than that of the DFT. In this lab, we won?t care about the mechanics of the

 

1

 


 

algorithm (hurrah!). We are simply concerned with how to use the FFT function in MATLAB to

 

determine frequency domain information of signals.

 


 

1.2 Background: How does the FFT ?Read? Frequency Information?

 

Before starting the lab, it may be useful to gain some insight as to how the discrete Fourier

 

transform extracts frequency information from a discrete input signal. The FFT uses a process

 

called correlation, which utilizes a series of sinusoidal basis signals of varying frequency. The

 

FFT determines how closely the input signal ?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 Useful Lecture Series:

 

For more information on the DFT and MATLAB FFT function, refer to the following series of

 

lecture videos by David Dorran. They are very useful for gaining an intuitive understanding of the

 

DFT and FFT, and will prove useful in completing this lab.

 


 

[1] Using MATLAB?s FFT Function (17:40)

 

https://www.youtube.com/watch?v=dM1y6ZfQkDU

 

[2] Using MATLAB?s FFT Function (2) ? Zero Padding and Windowing (23:06)

 

https://www.youtube.com/watch?v=V_dxWuWw8yM

 

[3] How the Discrete Fourier Transform (DFT) Works - An Overview (4:23)

 

https://www.youtube.com/watch?v=h6QJLx22zrE

 


 

2

 


 

1.4 Pre-Lab Exercise: Decomposing Signals into Sinusoids

 

Rewrite the following signals as a sum of sinusoids (that is, expand them as much as possible).

 

Identify the amplitude, frequency, and phase of each sinusoidal component. Your answers to this

 

part should help you in the Warmup section. In your write-up, provide your answers with

 

explanations or scans of your working.

 


 

2. Warmup (40 pts.)

 

In this section, we will use MATLAB?s FFT function to analyze basic signals ? ones that can be

 

expressed mathematically. In the third component of this lab, we will 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.

 


 

2.1 FFT of a Single Sinusoid with Integer Frequency

 

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:

 


 

We can see that the frequency is 50 Hz and the phase is /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 one peak, since the signal is composed of a single

 

sinusoid.

 

1)

 

2)

 

3)

 

4)

 


 

Define the variable fs for the sampling frequency and assign it 1000 samples/second.

 

Construct a time vector that is 1.5 seconds long (exactly 1500 samples).

 

Define (), the sinusoid described above, as the vector xx.

 

Run the FFT function. The function fft() only takes one argument: the input signal

 

xx. We will use a semicolon to suppress the output (otherwise the command window

 

will fill with complex numbers). By convention, frequency domain variables are

 

denoted with capital letters, so we will use XX to store the outputs.

 

>> XX = fft(xx);

 


 

3

 


 

5) You should notice that the output XX is defined as a ?1x1500 complex double.? This is

 

again because the outputs of the FFT are complex numbers. There are 1500 frequency

 

bins because there are 1500 samples in our signal. To make sense of the outputs, we

 

need to take the magnitude of this vector.

 

a. Create a new vector that stores the values of abs(XX) and plot the result. The

 

resulting plot is called the magnitude spectrum, where the x-values are

 

frequency bins and the y-values are their corresponding magnitudes. Note: By

 

convention, MATLAB?s indexing starts at 1 instead of 0. In the plot of the

 

magnitude spectrum, the bin values will start at 1. To make them start at zero,

 

use >> plot(0:length(XX)-1,abs(XX))

 

b. 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.

 

What is the frequency bin (index) corresponding to the left peak? You might

 

want to zoom in on it first. Again, be sure that the plot starts at zero. What is

 

the value of the magnitude spectrum everywhere else?

 

c. To find the frequency in Hz of the sinusoidal component corresponding to a

 

peak in the magnitude spectrum, we can use the following formula:

 

f = bin*fs/length(XX)

 

In this formula, bin is the frequency bin corresponding to the peak, fs is the

 

sampling frequency, and length(XX) is the total number of bins, equivalent

 

to the number of samples in the signal. Calculate the frequency. Does this match

 

the frequency of x(t)?

 

6) The FFT can also be used to extract phase information from the signal. Recall that XX

 

is a vector of complex numbers, which have both magnitude and phase. Hence, we can

 

also use the FFT to plot the phase spectrum of our signal.

 

a. Create a new vector that will store the values of phase(XX). Now run the

 

command >> phase(bin+1), where bin is the index you identified in the

 

previous section. Then divide the output by . How does this value compare

 

with the phase of ()?

 

b. Plot the phase spectrum of the signal. What do you notice about the values in

 

the plot? Are they bounded? If so, between what values?

 


 

4

 


 

7) In parts (5) and (6), we identified the frequency and phase of (). It turns out that we

 

can also use the FFT to determine the exact magnitude of (). To do this, we divide

 

the values of the magnitude spectrum by half the length of (). This can be achieved

 

by running the following command:

 

A = 2*abs(XX(bin+1))/length(XX);

 

Here, A is the amplitude of the sinusoidal component and bin is the index where the

 

magnitude spectrum peaks.

 

In this case, what is the value of the amplitude A? Does this value match the

 

amplitude of ()?

 


 

2.2 FFT of Multiple Sinusoids with Integer Frequency

 

Now we will use the FFT function to analyze more complicated signals. In the pre-lab, we were

 

given the functions:

 


 

These are clearly not sinusoids, but can be decomposed into a sum of finitely many sinusoids.

 

MATLAB?s FFT function will be useful in determining whether your work from the pre-lab was

 

correct. Complete the following steps:

 

1) Using the same sampling frequency and duration from 2.1, create the vectors x1 and

 

x2, storing 1500 samples of the signals 1 () and 2 ().

 

2) Plot the magnitude spectrum of x1 and x2. Unlike in 2.1, we should observe multiple

 

peaks, each corresponding to a single sinusoidal component. Comment on the

 

appearance of each spectrum. Do you notice any differences in the number of peaks for

 

these signals?

 

*Hint: In the magnitude spectrum of x1, zoom in on the zeroth bin. What is the

 

magnitude? Why might this be the case? (Think about the result from the pre-lab).

 

3) Locate the indices corresponding to the peaks in the magnitude spectra for each signal.

 

Use the formulas from the previous section to identify frequency components of each

 

signal. Do these agree with your results from the pre-lab?

 

4) Using the angle() function, identify the phases of the sinusoidal components of

 

1 () and 2 (). Do these agree with your results from the pre-lab?

 

5) Lastly, identify the exact amplitudes of the sinusoidal components of 1 () and 2 ().

 

Do these agree with those amplitudes found in the pre-lab?

 

5

 


 

2.3 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 section 2.1, 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 xx to represent the given sinusoid.

 

2) Use fft() and assign the output to XX. 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. Obtain the bin value where this maximum occurs. Use the

 

command f = bin*fs/length(XX) to find the frequency at bin (the bin value

 

you identified). Does this exactly match the frequency 50.5 Hz?

 


 

2.4 Zero Padding

 

You may have realized in section 2.3 that the resolution of the magnitude spectrum was not great.

 

At the bin value you identified, the frequency may not have perfectly matched the frequency 50.5

 

Hz. To correct this, one might think to increase the sampling frequency, but this does not correct

 

the error. To more accurately measure the frequency of our signal using the FFT, we can use a

 

process called zero padding.

 

Zero padding works by concatenating a large sample of zeros (using the zeros(1,N) function,

 

where N is the number of zeros appended) at the end of the signal being analyzed. For example,

 

say our signal is () = 3 cos(2(2.5) + /4). The following script uses vector notation to add

 

15,000 zeros to the end of xx. Here the sampling frequency is 1000 samples/second and t is a

 

time vector of 1500 samples. The results are plotted below:

 


 

xx = 3*cos(2*pi*2.5*t+ pi/4);

 

xx = [xx zeros(1,15000)];

 


 

6

 


 

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. For more information on the mechanics of zero

 

padding, refer to the following video:

 

[4] How DFT Zero Padding Works

 

https://www.youtube.com/watch?v=7yYJZfc5q-I

 


 

Complete the following steps:

 

1) For our signal in section 2.3, append it with 15,000 zeros, and run the FFT again.

 

Observe the effect on the appearance of the magnitude spectrum and calculate the

 

frequency for the bin value where the magnitude is maximized. Is this closer to the

 

actual value?

 

2) The appearance of the magnitude spectrum may be somewhat odd compared to what

 

we have seen before. There should appear to be a main ?lobe? and several smaller ?side

 

lobes? on the sides. How does the width of these side lobes compare to the main lobe?

 

3) Try zero padding with 25,000 zeros, 50,000 zeros, and 100,000 zeros. Using the subplot

 

function, compare the appearance of the magnitude spectrum. What is the effect of

 

increasing the number of zeros in the accuracy of the frequency calculation? Does the

 

7

 


 

error increase or decrease? In the case of our signal, how many zeros would need to be

 

padded for the formula f = bin*fs/length(XX) to give the exact value of the

 

frequency of our input signal?

 


 

2.5 Windowing

 

When analyzing signals with sinusoidal components of non-integer frequencies, we will generally

 

see magnitude spectra with ?side lobes? surrounding the main peaks where frequency components

 

occur. You have seen in the previous example that regardless of how many zeros we append at the

 

end of our signal, we will still observe side lobes. When analyzing the FFT of ECG heart rate data

 

in Lab Exercise 3.1, these side lobes may actually interfere with our calculations of the unknown

 

heart rate. One method for reducing the size of side lobes is called windowing.

 


 

8

 


 

Suppose we have the signal () = 3 cos(2(2.5) + /4) which we padded with zeros in the

 

previous section. All we are doing is multiplying our input signal with a bell-shaped function (the

 

window) to reduce the amplitude of the endpoints of the function, resulting in 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 in the image on the previous page. If you

 

read the documentation of the hanning() function in MATLAB, you will notice that it takes

 

only one argument (the length of the input vector) and outputs a column vector. Hence, in order to

 

multiply the values of the Hann Window by the values of our input signal, they need to be of the

 

same dimension. Hence, we need to transpose the Hann vector (i.e. hanning(length(xx))?).

 

For a more in-depth explanation on windowing, refer to the following:

 

[5] DFT Windowing Explanation and Demo

 

https://www.youtube.com/watch?v=PYd-pzaZpYE

 


 

Complete the following steps:

 

1) Returning to our signal () = cos(2(50.5) + 0.25) from section 2.3 and 2.4, we

 

want to reduce the size of the ?side lobes? in the plot of the magnitude spectrum. Using

 

the script above, implement the Hann window on (). 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

 

() after windowing and zero padding. Make a side-by-side plot of the magnitude

 

spectra, zoomed in at the peak value, for XX with windowing and XX without

 

windowing. What can we conclude about the effect windowing has on the FFT? How

 

does the magnitude of the maximum peak change with windowing? Is it higher or lower

 

than that of the maximum peak without windowing? Has the width changed at all?

 


 

9

 


 

3. Lab Exercise (50 pts)

 

Determining Average Heart Rate Using FFT

 


 

A plot showing a sample of heart beat waveforms from one of the ECG measurements in the zip file.

 


 

In this exercise, you will employ the fft() function in MATLAB to automatically extract the

 

heart rate from electrocardiogram (ECG) measurements. The zip file on Canvas contains four ECG

 

signals (obtained from http://physionet.org/cgi-bin/atm/ATM). Decompress the zip file to obtain

 

the ECG signals in the .mat files. Load each .mat file, and the ECG signal samples are contained

 

in the val variable. The sampling rate for all data files is 125 Hz. Your goal will be to produce

 

four plots, showing how the average heart rate changes over each signal?s duration.

 


 

10

 


 

Key Steps to a Solution

 

Note: You are allowed and encouraged to use anything you have learned to help the FFT

 

identify the heart rate. This may include zero padding, windowing, filters, etc. In addition,

 

you can choose whatever window length or overlap length you feel will give the best result.

 

1. Load each signal in MATLAB. Because there are four .mat files (101am, 105am, 107am,

 

110cm), it might be a good idea to import each file, one at a time, and store val into a new

 

vector (one for each signal).

 

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 use a window

 

about 2 minutes long.

 

3. In order to plot the average heart rate over time, you will need to iterate through the

 

windows and take the FFT for each window. It is suggested that there be an overlap

 

between each window as the function iterates. This overlap is suggested to be about 10

 

seconds.

 

4. Take the absolute value of the FFT of the sample (abs(XX)) 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. Record the index/bin of the peak you identified. Convert this bin number to frequency in

 

Hz, and then convert this frequency to beats/minute.

 

6. Your function/script should store these average heart rates into a vector, and then produce

 

a plot showing the change in the average heart rate over time.

 

7. You will do this for each .mat file.

 


 

11

 


 

 

Paper#9210654 | Written in 27-Jul-2016

Price : $19
SiteLock