This shows you the differences between two versions of the page.
matlab [2015/05/14 10:48] schultz [Using the freqz command] |
matlab [2023/01/18 16:12] (current) scott Figure scale |
||
---|---|---|---|
Line 1: | Line 1: | ||
- | ====== Matlab Basic Functions ====== | + | ====== MATLAB Basic Functions ====== |
* len | * len | ||
* plot | * plot | ||
Line 11: | Line 11: | ||
You can type <code> help 'name'</code> to get a description of the particular function. | You can type <code> help 'name'</code> to get a description of the particular function. | ||
- | Let's do an example to use these basic functions. We are going to plot 4 different signals on 2 separate plots. | + | Let's do an example to use these basic functions. We are going to plot 4 different signals on 2 separate plots. |
* Plot 1 will contain a square wave and sine waves with amplitude of 2 and frequency of f=1 kHz. | * Plot 1 will contain a square wave and sine waves with amplitude of 2 and frequency of f=1 kHz. | ||
* Plot 2 will contain a square wave and sine waves with amplitude of 2 and frequency of f=10 kHz. | * Plot 2 will contain a square wave and sine waves with amplitude of 2 and frequency of f=10 kHz. | ||
- | Here is the Matlab code for the first plot. | + | Here is the MATLAB code for the first plot. |
<code> | <code> | ||
f1=1e3; | f1=1e3; | ||
Line 34: | Line 34: | ||
</code> | </code> | ||
- | Here is a description of the Matlab code. | + | Here is a description of the MATLAB code. |
- We use a variable for the signal frequency (//f1//) and the sampling frequency. | - We use a variable for the signal frequency (//f1//) and the sampling frequency. | ||
- | - We create a vector for the time (//t//). It has an increment of the sampling time (1/fs) and a total length of 10 periods (//10/f1//). | + | - We create a vector for the time (//t//). It has an increment of the sampling time (1/fs) and a total length of 10 periods (//10/f1//). |
- We create the vectors for the sine and square waves. | - We create the vectors for the sine and square waves. | ||
- We use the //subplot// function to plot the signals in the upper half of the plot window. | - We use the //subplot// function to plot the signals in the upper half of the plot window. | ||
- | - The //hold on// makes it so that the next plot command overlays onto the same plot. //hold off// does the opposite. | + | - The //hold on// makes it so that the next plot command overlays onto the same plot. //hold off// does the opposite. |
- Use //xlabel// and //ylabel// to add in axes. | - Use //xlabel// and //ylabel// to add in axes. | ||
- | The second part of the code is similar and just changes the frequency, the time, and the plot location. Here is the resulting Matlab code. | + | The second part of the code is similar and just changes the frequency, the time, and the plot location. Here is the resulting MATLAB code. |
<code> | <code> | ||
f2=10e3; | f2=10e3; | ||
Line 61: | Line 61: | ||
The resulting plot should look like the following. | The resulting plot should look like the following. | ||
- | {{ :example1.jpg?100 |}} | + | {{ :example1.jpg?600 |}} |
- | ====== Matlab Filtering Functions ====== | + | ====== MATLAB Filtering Functions ====== |
* butter: Butterworth filter design | * butter: Butterworth filter design | ||
* filter: Time domain filtering of a signal | * filter: Time domain filtering of a signal | ||
Line 79: | Line 79: | ||
</code> | </code> | ||
- | on the Matlab command line. Read through the description. Let’s look at parts of the description and what it means. | + | on the MATLAB command line. Read through the description. Let’s look at parts of the description and what it means. |
**Description**: “butter Butterworth digital and analog filter design.” | **Description**: “butter Butterworth digital and analog filter design.” | ||
Line 103: | Line 103: | ||
</code> | </code> | ||
- | on the Matlab command line. Read through the description. Let’s look at parts of the description and what it means. | + | on the MATLAB command line. Read through the description. Let’s look at parts of the description and what it means. |
**Description**: “One-dimensional digital filter.” | **Description**: “One-dimensional digital filter.” | ||
Line 132: | Line 132: | ||
You should get a plot that looks like | You should get a plot that looks like | ||
- | {{ :matlab_help_fig1.jpg?100 |}} | + | {{ :matlab_help_fig1.jpg?600 |}} |
- | If we create a low pass filter with a cutoff frequency of //fcut=1.5kHz// then the signal should pass through the filter. Let’s make some mistakes and see what happens and see if we can figure out what we did wrong. We are doing this to learn about debugging. | + | If we create a low pass filter with a cutoff frequency of //fcut=1.5kHz// then the signal should pass through the filter. Let’s make some mistakes and see what happens and see if we can figure out what we did wrong. We are doing this to learn about debugging. |
**Mistake**: Pass the cutoff frequency as the actual frequency. | **Mistake**: Pass the cutoff frequency as the actual frequency. | ||
Line 162: | Line 162: | ||
</code> | </code> | ||
The resulting plot should look something like. | The resulting plot should look something like. | ||
- | {{ :matlab_help_fig2.jpg?100 |}} | + | {{ :matlab_help_fig2.jpg?600 |}} |
- | The signal y has a lower amplitude. This mean that it is being filtered. We must have made an error in the cutoff frequency. If we make the following change | + | The signal y has a lower amplitude. This mean that it is being filtered. We must have made an error in the cutoff frequency. If we make the following change |
<code> | <code> | ||
Wn=Wn/(2*Fs); | Wn=Wn/(2*Fs); | ||
</code> | </code> | ||
The following plot shows that it is now worse because we were supposed to divide by Fs/2 not 2*Fs. | The following plot shows that it is now worse because we were supposed to divide by Fs/2 not 2*Fs. | ||
- | {{ :matlab_help_fig3.jpg?100 |}} | + | {{ :matlab_help_fig3.jpg?600 |}} |
The correct value is | The correct value is | ||
Line 177: | Line 177: | ||
Resulting in the correct plot. | Resulting in the correct plot. | ||
- | {{ :matlab_help_fig4.jpg?100 |}} | + | {{ :matlab_help_fig4.jpg?600 |}} |
- | To finish up let’s make sure that we can use the butter command to perform a bandpass filter. Create a signal that is x=sin(2*pi*1e3*t) + sin(2*pi*2e3*t) + sin(2*pi*3e3*t) and has the same sampling frequency of Fs=10kHz. Use a bandpass filter to pass the f=2kHz signal through and attenuate the other two. The resulting plots should look like the following. | + | To finish up let’s make sure that we can use the butter command to perform a bandpass filter. Create a signal that is x=sin(2*pi*1e3*t) + sin(2*pi*2e3*t) + sin(2*pi*3e3*t) and has the same sampling frequency of Fs=10kHz. Use a bandpass filter to pass the f=2kHz signal through and attenuate the other two. The resulting plots should look like the following. |
- | {{ :matlab_help_fig5.jpg?100 |}} | + | {{ :matlab_help_fig5.jpg?600 |}} |
===== Using the fft command ===== | ===== Using the fft command ===== | ||
- | MATLAB has a built-in function for performing numerical Fourier transforms called the fast Fourier transform (FFT). The FFT is actually a discrete Fourier transform (DFT). We can relate the DFT to the continuous time domain. This [[https://www.youtube.com/watch?v=z7X6jgFnB6Y| youtube video]] explains the conversion between the FFT and frequency. MATLAB also has some help on using the FFT. In the command window type | + | MATLAB has a built-in function for performing numerical Fourier transforms called the fast Fourier transform (FFT). The FFT is actually a discrete Fourier transform (DFT). We can relate the DFT to the continuous time domain. This [[https://www.youtube.com/watch?v=z7X6jgFnB6Y| youtube video]] explains the conversion between the FFT and frequency. MATLAB also has some help on using the FFT. In the command window type |
<code> help fft</code> | <code> help fft</code> | ||
Then click on the | Then click on the | ||
Line 190: | Line 190: | ||
doc fft link | doc fft link | ||
</code> | </code> | ||
- | Go through the fft MATLAB document and make sure that you can use the function properly. | + | Go through the FFT MATLAB document and make sure that you can use the function properly. |
Let's go through another example besides the one provided in the MATLAB help to make sure that you understand how to use the FFT. | Let's go through another example besides the one provided in the MATLAB help to make sure that you understand how to use the FFT. | ||
- | **You need to create an frequency domain plot of the sum of two sinusoidal signals. The first sinusoidal signal has frequency of f1=3 kHz and amplitude of 1. The second sinusoidal signal has frequency of f1=6 kHz and amplitude of 3. You need to use a sampling frequency of Fs=100 kHz and a total of 20 cycles of the sine wave.** | + | You need to create a frequency domain plot of a sinusoidal signal with frequency of f1=3 kHz and amplitude of 1. We are using a sinusoidal signal because we know that the Fourier transform of a sine wave is a delta function. So the FFT should produce a spike at a frequency of 3 kHz and have an amplitude of 1. |
- | Here is the code of how to create the time signal | + | Here is the MATLAB code to create the time signal |
<code> | <code> | ||
Fs=100e3; %sampling rate | Fs=100e3; %sampling rate | ||
- | f=1e3; %signal frequency | + | f1=3e3; %signal frequency |
- | A=1; %signal amplitude, This is peak amplitude not peak-to-peak | + | A1=1; %signal amplitude, This is peak amplitude not peak-to-peak |
- | t=0 : 1/Fs : 20/f; | + | t=0 : 1/Fs : 10/f1; |
- | V=A*sin(2*pi*f*t); | + | y=A1*sin(2*pi*f1*t); |
</code> | </code> | ||
- | Now you need to create the Fourier transform of the signal using the **fft** command. You need to create an x-axis variable for the frequency with the correct scaling and you need to scale the amplitude. You should only be using the first half of the data since we only need to plot the positive frequencies. Your plot should look like this. | + | Here is the plot of our signal. |
- | {{ :sin_fft.jpg?200 |}} | + | |
- | Make sure that the x axis is scaled correctly by zooming in on the spike and verifying that it is at 1 kHz. The maximum frequency is Fs/2=50 kHz. | + | {{ :fft_example1.jpg?600 |}} |
+ | |||
+ | Now you need to create the Fourier transform of the signal using the FFT command. | ||
+ | <code> | ||
+ | Y=fft(y) | ||
+ | plot(abs(Y)) | ||
+ | </code> | ||
+ | Notice that we plotted the absolute value of Y since it is complex. | ||
+ | |||
+ | {{ :fft_example2.jpg?600 |}} | ||
+ | |||
+ | Notice that there are two spikes and that the amplitude is not equal to our amplitude of 1. | ||
+ | |||
+ | Let's start with the frequency. We only want the positive frequencies so we are only going to use the first half of the vector. We know that maximum frequency of a sampled signal is Fs/2. So we create an x-axis vector that goes from 0 to Fs/2 and has a length of half of Y. Here is the resulting MATLAB code. | ||
+ | <code> | ||
+ | L=length(Y); %the total length of our fft | ||
+ | Y1=Y(1:L/2); %we take the first half of the vector | ||
+ | f=linspace(0,Fs/2,L/2); %we create the frequency vector | ||
+ | plot(f,abs(Y1)) | ||
+ | </code> | ||
+ | The linspace command creates a vector with equally spaced elements between a starting and an ending value and with a specified number of elements. Here is the resulting plot. | ||
+ | {{ :fft_example3.jpg?600 |}} | ||
+ | You can zoom in on the spike and verify that it has a frequency close to 3kHz. | ||
+ | |||
+ | Our signal does not have an amplitude greater than 100 so we obviously need to scale the amplitude. We simply need to divide the signal by the number of elements. | ||
+ | <code> | ||
+ | Y2=Y1/(L/2); | ||
+ | plot(f,abs(Y2)) | ||
+ | </code> | ||
+ | Here is the resulting plot. | ||
+ | {{ :fft_example4.jpg?600 |}} | ||
+ | |||
+ | Zoom in on the spike to make sure that it is centered close to 3kHz. | ||
+ | |||
+ | //NOTE: As you learned by watching the [[https://www.youtube.com/watch?v=z7X6jgFnB6Y| youtube video]], the FFT creates an array of frequency domain spectrum samples. The number of samples is the same as the length of your data (N) and the frequency range is based on the sampling frequency fmax=Fs/2. This means that the frequency bin is Fs/(2*N).// | ||
+ | |||
+ | In our example the number of samples is L/2=167 and the frequency band is Fs/2=50kHz. This means that our frequency resolution is 299.4 Hz. If we increase the number of points in our time domain signal without increasing the total time of the signal then we are actual changing the sampling frequency, which would increase our total frequency band but not change our frequency resolution. In order to change the resolution we need to increase the total time width of our signal without changing the sampling frequency. | ||
+ | |||
+ | This MATLAB code produces a sine wave with 10 periods and a corresponding frequency resolution of 299 Hz. | ||
+ | <code> | ||
+ | Fs=100e3; %sampling rate | ||
+ | f1=3e3; %signal frequency | ||
+ | A1=1; %signal amplitude, This is peak amplitude not peak-to-peak | ||
+ | t=0:1/Fs:10/f1; | ||
+ | y=A1*sin(2*pi*f1*t); | ||
+ | Y=fft(y); | ||
+ | L=length(Y); % the total length of our fft | ||
+ | Y1=Y(1:L/2); %we take the first half of the vector | ||
+ | Y2=Y1/(L/2); | ||
+ | f=linspace(0,Fs/2,L/2); | ||
+ | subplot(2,1,1) | ||
+ | plot(t,y) | ||
+ | subplot(2,1,2) | ||
+ | plot(f,abs(Y2)) | ||
+ | </code> | ||
+ | Here are the resulting time domain and frequency domain plots. | ||
+ | {{ :fft_example5.jpg?600 |}} | ||
+ | |||
+ | This MATLAB code produces a sine wave with 100 periods and a corresponding frequency resolution of 29.9 Hz. | ||
+ | <code> | ||
+ | Fs=100e3; %sampling rate | ||
+ | f1=3e3; %signal frequency | ||
+ | A1=1; %signal amplitude, This is peak amplitude not peak-to-peak | ||
+ | t=0:1/Fs:100/f1; | ||
+ | y=A1*sin(2*pi*f1*t); | ||
+ | Y=fft(y); | ||
+ | L=length(Y); % the total length of our fft | ||
+ | Y1=Y(1:L/2); %we take the first half of the vector | ||
+ | Y2=Y1/(L/2); | ||
+ | f=linspace(0,Fs/2,L/2); | ||
+ | subplot(2,1,1) | ||
+ | plot(t,y) | ||
+ | subplot(2,1,2) | ||
+ | plot(f,abs(Y2)) | ||
+ | </code> | ||
+ | Here are the resulting plots. | ||
+ | {{ :fft_example6.jpg?600 |}} | ||
- | If you are not able to create this plot then you need to go through the MATLAB documentation again. | ||
===== Using the freqz command ===== | ===== Using the freqz command ===== | ||
Line 217: | Line 292: | ||
**Description**: “freqz Frequency response of digital filter.” | **Description**: “freqz Frequency response of digital filter.” | ||
- | **Why we care**: We want the frequency response of a digital filter. So it looks like we found the correct MATLAB command. | + | **Why we care**: We want the frequency response of a digital filter. So it looks like we found the correct MATLAB command. |
**Description**: “[H,W] = freqz(B,A,N) returns the N-point complex frequency response vector H and the N-point frequency vector W in radians/sample of the filter:” | **Description**: “[H,W] = freqz(B,A,N) returns the N-point complex frequency response vector H and the N-point frequency vector W in radians/sample of the filter:” | ||
- | **Why we care**: H is the complex frequency response. N is the number of points that are calculated. | + | **Why we care**: H is the complex frequency response. N is the number of points that are calculated. |
- | We don't really want the plot in terms of radians/sample. We would prefer to create the plot in terms of frequency. So we are going to scan down the description until we find the form of freqz that is in terms of frequency. | + | We don't really want the plot in terms of radians/sample. We would prefer to create the plot in terms of frequency. So we are going to scan down the description until we find the form of freqz that is in terms of frequency. |
- | **Description**:"H = freqz(...,F,Fs) returns the complex frequency response at the frequencies designated in vector F (in Hz), where Fs is the sampling frequency (in Hz)." | + | **Description**: "H = freqz(...,F,Fs) returns the complex frequency response at the frequencies designated in vector F (in Hz), where Fs is the sampling frequency (in Hz)." |
- | **Why we care**: This tells us how to create the plot in terms of frequency rather than samples. | + | **Why we care**: This tells us how to create the plot in terms of frequency rather than samples. |
Here is an example to show how to use the freqz command. | Here is an example to show how to use the freqz command. | ||
Line 257: | Line 332: | ||
Here is the resulting plot. | Here is the resulting plot. | ||
- | {{ :freqz_1.jpg?200 |}} | + | {{ :freqz_1.jpg?600 |}} |
- | Notice how the filter is flat over the band between the two frequencies and then slopes down outside of this band. This is the shape of the frequency response of a bandpass filter. | + | Notice how the filter is flat over the band between the two frequencies and then slopes down outside of this band. This is the shape of the frequency response of a bandpass filter. |
- | Often the data is scaled to be in decibels so that it is easier to distinguish small changes. When you do this you often need to scale the y-axis of the plot to your desired range of interest. | + | Often the data is scaled to be in decibels so that it is easier to distinguish small changes. When you do this you often need to scale the y-axis of the plot to your desired range of interest. |
<code> | <code> | ||
H_db=20*log10(abs(H)); | H_db=20*log10(abs(H)); | ||
Line 267: | Line 342: | ||
axis([0 5000 -50 3]) | axis([0 5000 -50 3]) | ||
</code> | </code> | ||
- | Here is the plot in decibel scale. You can zoom in on the 1kHz and 2kHz frequencies and see that these are the 3dB points. | + | Here is the plot in decibel scale. You can zoom in on the 1kHz and 2kHz frequencies and see that these are the 3dB points. |
- | {{ :freqz_3.jpg?200 |}} | + | {{ :freqz_3.jpg?600 |}} |
You can change the order of the filter and see that the filter becomes steeper and change the frequency corners and see the change in the width of the filter. | You can change the order of the filter and see that the filter becomes steeper and change the frequency corners and see the change in the width of the filter. | ||