Remember that the basic operation of the system involves the following operations:
One of the essential parts of this basic operation is the bank of bandpass filters. You already designed and tested this bank of filters on sampled data that was provided to you during Lab #5 in ECEN 380. In this task, you are going to make sure that your filter bank is working by testing it with simulated signals at the different player frequencies.
This milestone requires the design, analysis, and testing of the bank of 10 bandpass filters in MATLAB. You will look at both frequency and time domain analysis of the involved signals and filters.
The 10 bandpass filters should be IIR filters with the following specifications.
As you know, the laser tag system supports 10 players each with an assigned frequency in the range of 1.4 - 4.2 kHz. If the receiver is being illuminated or 'hit' by a laser tag gun then the sampled signal will consist of both a square wave at a particular frequency and noise. The system needs to analyze the received signal to determine whether the sampled signal contains a 'hit' by a laser tag gun. As the shooting gun gets farther away from the receiver, the amplitude of the square wave signal from the shooting gun gets smaller. Since the noise amplitude stays the same, this corresponds to a lower SNR, and it gets harder to distinguish the signal from the noise. This figure shows a square wave added to a random noise signal:
It will be easy for you to see the square wave in the close range signal, when the SNR is high. However, for longer-range signals, it will just look like noise even though there is a square wave embedded in the signal. The beauty of your signal processing system will be in its ability to detect a hit in the noisy received signal even at relatively large ranges.
The 'hit' signal has a frequency spectrum concentrated at the modulating frequency, while the noise is spread out over a wider frequency band (as we'll see in Task 2). A simple method to extract the signal from the noise would be to compute a Fourier transform of the received signal, and then see if there is a spike in frequency at one of the player frequencies.
Unfortunately, computing a Fourier transform of the incoming signal requires too much computation time to enable it to be performed in real-time by our system. Therefore, we are effectively going to make a 'coarse' Fourier transform using our bank of 10 bandpass filters, and see if there is a spike in the signal energy coming out of any of the filters. This is computationally much more efficient, and can be implemented in real time (provided we keep the order of our filter low enough).
You should start by reviewing the lab from ECEN 380 in which you designed digital filters. For convenience, here are links to both Lab #5 in ECEN 380 and the supplementary filter design handout:
Lab #5 from ECEN 380 Supplementary filter design handout
We are going to be using an Infinite Impulse Response (IIR) digital filter. Wikipedia has a description of this type of filter (http://en.wikipedia.org/wiki/Infinite_impulse_response). You should read through this description. The IIR digital filter is characterized by two sets of coefficients called a and b.
There are lots of different filter designs that have already been developed. A few of the common types of filters are Butterworth, Chebyshev, and Elliptic. MATLAB has functions that produce the coefficients (a and b) for these filters. Let’s use the Butterworth filter for this example. Review the description of the Butterworth filter on Wikipedia (https://en.wikipedia.org/wiki/Butterworth_filter). At the bottom of the Wikipedia page there is a comparison to the other common types of digital filters.
The MATLAB command is 'butter'. Typically when you are learning how to use a new command in MATLAB you type
help butter
on the MATLAB command line to get some help. In this case here is a more detailed help on using the butter and filter MATLAB commands.
IMPORTANT: The 'butter' command in MATLAB takes frequencies in units of 'half-cycles/sample'. You can think of this as a normalized frequency in the range of 0 to 1, where one corresponds to half your sampling frequency. It simply tells you how many half-cycles your signal goes through per sample. The easiest way to get these normalized frequencies in half-cycles/sample is to divide your frequency in Hz by half the sample rate. I suggest keeping all of your frequencies in your MATLAB code as human readable Hz, and then simply divide your frequencies in Hz by half your sample rate when you pass them to the butter command.
The design parameter (what you can change) are (1) the filter type (i.e., Butterworth, Chebyshev, etc.), (2) the filter order, and (3) the cutoff frequencies.
To make it easier for the TAs to help we are going to limit the number of design parameters. We want you all to use the Butterworth filter design. To limit computational complexity so your filters can operate in real time on the ZYBO board, we are also going to limit the length of each IIR filter to 11 (i.e., the b and a vectors each have a length of 11). Note that this corresponds to a Butterworth filter order of 10, since a polynomial with 11 coefficients has 10 roots (the numbers of poles or zeros). This means that the only design parameters we are varying are the cutoff frequencies of the filters.
Each Butterworth bandpass filter has two frequency corners. It should pass signals with frequency between these two frequencies and attenuate all other frequency components. Review the code provided in ECEN 380 Lab #5 Task #3 that shows you how to set the corner frequencies for player 1, create the bandpass Butterworth filter, and generate the a and b coefficients for the filter.
Once you have created the filter coefficients for each of your filters, please save the coefficients off into double precision ascii files. This will allow you to easily load the filters into the ZYBO board in the future. If the MATLAB variable a1 holds your a coefficients for the player 1 frequency filter, then the following code would save those coefficients to an ascii file called a1.txt:
save a1.txt a1 -ascii -double
In a future milestone, you will need to convert your MATLAB filter coefficients into C code for use on the ZYBO board. A script will be provided to ease this conversion. The script expects the following convention:
Note: If we were making a true Coarse Frequency Domain calculation the filter would need to cover the entire band between the different player frequencies. However, in this system we have a discrete set of allowable frequencies. The filter just needs to pass one of these allowable frequencies to see if the receiver is 'hit'. Therefore, the bandwidth of the filter just needs to cover the specific allowable frequency that it is assigned. This means that the narrower the bandwidth the more noise that is rejected and the better the performance of the system. However, this assumes that the frequency of the transmitter does not change at all. Therefore, we need to create some bandwidth to allow wobble in the player frequency. How much bandwidth is actually needed would need to be determined experimentally. In this example we created a bandwidth of 50 Hz. This can be optimized after the system is built if you want to optimize the range of your system.
To determine whether the Butterworth filters are designed properly we need to analyze the filters. The filter analysis will include both a frequency domain and time domain analysis.
In the frequency domain each bandpass filter should have a flat transfer function around the frequency that the filter is supposed to pass and then decrease as the frequency gets farther away from the frequency corners. You can plot the frequency domain transfer function using the MATLAB freqz function.
NOTE: I suggest using the following form of freqz: freqz(b, a, f, f_s) where b and a are the filter coefficients, f is a frequency vector containing exactly the frequencies (in Hz) that you want to evaluate the filter at, and f_s is the sampling frequency in Hz. Then you can simply plot the output of freqz against your frequency axis f.
You should create a frequency domain plot for each of the bandpass filters and overlay them on the same plot. The result should look similar to the following.
The horizontal axis should be in frequency and cover the entire allowable frequency band (0<f<5kHz).
The bandpass filters will eventually be implemented (in a later milestone) in the time domain on the ZYBO board, since frequency-domain processing would be too computationally expensive (as we've previously mentioned). We therefore will need to analyze the output of each filter in the time domain, and determine whether a hit was registered. We will now walk you through the basic process of testing your filter bank by (1) simulating a time domain square wave signal at one of the player frequencies, (2) filtering the square wave signal with your filter bank, (3) computing the energy (in the time domain) of each filtered signal, and (4) determining which frequency the hit was at based on the signal energies.
The basic steps of the time domain analysis are the following.
Let's walk through this process for the first player.
Step 1: Create a square wave
Use the MATLAB function 'square'. The square command is similar to the sin command except that it produces a square wave rather than a sine wave. Here is the MATLAB code to create the square wave:
Fs = 10e3; % The sampling frequency t = linspace(0, 0.2, 2000); % The time vector with a length of 0.2 seconds and 2000 total samples freq = 1471; % The frequency of player 1 x = square(2*pi*freq*t); % Create the time domain square wave
Step 2: Filter the square wave
As you know, the bandpass filters have two arrays of coefficients (a and b) that define the IIR digital filters. In a future milestone, you will implement these filters in the time domain using the queue code you developed in Milestone 1. However, we are going to keep things simple during our testing and simply use the built-in MATLAB function 'filter' to perform the actual filtering in this milestone. We need to filter the square wave signal that we created using each of our 10 bandpass filters (we called the square wave signal x in our previous code snippet). You already reviewed how to use the filter command when you learned how to use the butter MATLAB command. If you need another review go here.
Plot the output of the filter command (operating on your square wave signal) for each of the 10 bandpass filters. The following plot shows the output of the first 4 bandpass filters when the square wave input has a frequency of 1471 Hz (the frequency for player 1):
You can see from this plot that filters 2, 3, and 4 attenuate the signal a lot compared to filter 1 (look at the amplitudes of each filtered signal). The output of filter 1 has a much higher amplitude. If you zoom in on the output of filter 1 you should be able to see that the square wave was changed into a sine wave by the filter. (Do you understand why this happens? Remember that the square wave signal contains a sinusoid at the fundamental frequency, as well as higher frequency harmonics that give it a square shape. Our bandpass filter effectively takes out the higher harmonic frequencies that make it a square wave, leaving only the sinusoid at the fundamental frequency!)
Step 3: Compute the energy in each filtered signal
Looking at the outputs of the 10 filters, it is evident that when the square wave has a frequency that matches the frequency of a filter, the amplitude of the filtered signal will be much higher. You need a way of robustly determining which filtered signal has the highest average amplitude. However, if you simply take the largest value in each filtered signal, the system will be very susceptible to noise. A better way is to determine the total energy in each signal.
The total energy in a signal is found by squaring each sample and then summing all of the squared samples. In our implementation, we are looking at 200ms snapshots of our signals, because that is the duration of our shots. These 200ms snapshots are sampled at 10 ksamples/s, yielding 2000 total samples. Thus, the summation for our signal energy calculation is only going to be over 2000 samples.
Now that we have computed the energy of each filtered signal, we can create our coarse Fourier transform. We take the summation of the 2000 squared samples through each of the 10 bandpass filters and create an array. We then plot a bar graph for each of the array values. MATLAB has a command to create a bar graph called 'stairs'.
Your resulting bar graph produced using the square wave with the frequency for player 1 should look like this:
You should produce 10 of these plots, one for each of the 10 player frequencies.
The following items need to be shown to the TAs for pass off: