ECEN 390 - Laser Tag Project

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Matlab Basic Functions ======
  * len
  * plot
  * subplot
  * hold
  * sin
  * square
  * xlim
  * ylim

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

Here is the Matlab code for the first plot.
<code>
f1=1e3;          
fs=100e3;       
t=0:1/fs:10/f1;
y1=sin(2*pi*f1*t);
y2=square(2*pi*f1*t);

subplot(2,1,1)
hold off
plot(t*1e3,y1,'b')
hold on
plot(t*1e3,y2,'r')
hold off
ylabel('amplitude')
xlabel('time (ms)')
ylim([-1.2 1.2])
</code>

Here is a description of the Matlab code.
  - 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 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.
  - 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.

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>
f2=10e3;               
t=0:1/fs:10/f2;
y3=sin(2*pi*f2*t);
y4=square(2*pi*f2*t);

subplot(2,1,2)
hold off
plot(t*1e3,y3,'b')
hold on
plot(t*1e3,y4,'r')
hold off
ylabel('amplitude')
xlabel('time (ms)')
ylim([-1.2 1.2])
</code>

The resulting plot should look like the following.
{{ :example1.jpg?100 |}}


====== Matlab Filtering Functions ======
  * butter: Butterworth filter design
  * filter: Time domain filtering of a signal
  * fft: Fast Fourier Transform
  * freqz: complex frequency response of a digital filter


===== Using the butter and filter commands=====

==== butter ====
Start by typing
<code>
help butter
</code>

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

**Why we care**: We are designing a digital filter so this should help.

**Description**: “[B,A] = butter(N,Wn) designs an Nth order lowpass digital Butterworth filter and returns the filter coefficients in length N+1 vectors B (numerator) and A (denominator).”

**Why we care**: The B (numerator) and A (denominator) coefficients means that it calculates the coefficients of an IIR filter (see http://en.wikipedia.org/wiki/Infinite_impulse_response).

**Description**: “The cutoff frequency Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to half the sample rate.”

**Why we care**: When we design the cutoff frequencies we divide our desired cutoff frequency by half the sample rate to get it into the correct form for the butter function.

**Description**: “If Wn is a two-element vector, Wn = [W1 W2], butter returns an order 2N bandpass filter with passband W1 < W < W2.”

**Why we care**: Rather than simply designing a low pass filter we can calculate the IIR filter coefficients for a bandpass filter, which is what we want.

==== filter ====
Type
<code>
help filter
</code>

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

**Why we care**: We want to do a one dimensional filter so it looks like we have the correct function.

**Description**: “Y = filter(B,A,X) filters the data in vector X with the filter described by vectors A and B to create the filtered data Y.”

**Why we care**: Our sampled data is the X vector and the resulting sampled data is the Y vector.  

**Description**: “If a(1) is not equal to 1, filter normalizes the filter coefficients by a(1).”

**Why we care**: We need to make sure that the first element of our a vector is equal to 1. 

==== Using the butter and filter functions ====
Let’s do a simple example to make sure that we understand how to use the butter and filter MATLAB commands.  

We are going to start with a low pass filter to make sure that we have the cutoff frequency set properly.  

Plot a sinusoidal signal with a frequency of f1=1kHz, sampling frequency of Fs=10e3, and plot 10 periods of the signal.
<code>
Fs=10e3;
f1=1e3;
t=0:1/Fs:10/f1;
x1=sin(2*pi*f1*t);
plot(t,x1)
</code>

You should get a plot that looks like 
{{ :matlab_help_fig1.jpg?100 |}}

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.  
<code>
Wn=1.5e3;
[B,A]=butter(5,Wn);
</code>

We get an error message that say “…cutoff frequencies must be within the interval of (0,1).”
We would type help butter and find out that “Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to half the sample rate.”

**Mistake**: Divide the cutoff frequency by Fs rather than Fs/2.
<code>
Wn=1.5e3;
Wn=Wn/(Fs);
[B,A]=butter(5,Wn);

y=filter(B,A,x1);

subplot(211)
plot(t,x1)
ylabel('Original')

subplot(212)
plot(t,y)
ylabel('filtered')
</code>
The resulting plot should look something like.
{{ :matlab_help_fig2.jpg?100 |}}

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>
Wn=Wn/(2*Fs);
</code>
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 |}}

The correct value is
<code>
Wn=Wn/(Fs/2);
</code>

Resulting in the correct plot.
{{ :matlab_help_fig4.jpg?100 |}}

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 |}}

===== 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 
<code> help fft</code>
Then click on the 
<code>
doc fft link 
</code>
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.  

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

Here is the MATLAB code to create the time signal
<code>
Fs=100e3;  %sampling rate
f=1e3;     %signal frequency
A=1;       %signal amplitude, This is peak amplitude not peak-to-peak
t=0 : 1/Fs : 20/f;
V=A*sin(2*pi*f*t);
</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.
{{ :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.

If you are not able to create this plot then you need to go through the MATLAB documentation again.

===== Using the freqz command =====
Let's go through the MATLAB help and find the portions that are the most pertinent. Start by typing
<code>
help freqz
</code>
**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.

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

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)."

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

(1) We need to create the filter.
<code>
Fs=10e3;       % sampling frequency of 10kHz
Wn=[1e3 2e3];  % Bandpass filter with 1kHz < f < 2kHz
Wn=Wn/(Fs/2);  %divide by half of the sampling frequency

[B,A]=butter(4,Wn); %create the filter
</code>

(2) Set up the frequency vector.
<code>
F=0:Fs/(1000):Fs/2;     %create the array of frequencies
</code>
In this example we are creating an array of frequencies between 0 and the maximum frequency, which is half of the sampling frequency.

(3) Use the freqz command to create the complex frequency response.
<code>
H = freqz(B,A,F,Fs);    % create the complex frequency response   
</code>

(4) Plot the magnitude of the complex frequency response.
<code>
plot(F,abs(H))
</code>

Here is the resulting plot.
{{ :freqz_1.jpg?200 |}}

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.  
<code>
H_db=20*log10(abs(H));
plot(F,H_db)
axis([0 5000 -50 3]) 
</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.

{{ :freqz_3.jpg?200 |}}

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.
ECEN 390 - Laser Tag Project

User Tools

Site Tools

Page Tools
