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.
====== Milestone 2, Task 2: Design Decimating System and Analyze Noise Aliasing ======

===== Overview =====

In Task 2, you will design the system that takes your originally sampled signal (sampled at 100 ksamples/s on the ZYBO board) and down samples it to 10 ksamples/s so that it can be fed into your bank of 10 IIR filters to determine if there is a hit.

You will also be performing a noise analysis in this task. You will look at the spectrum of optical background noise sampled in a lab room at frequencies up to 50 kHz (the range that you can see using the data sampled at the full 100 ksamples/s), and then see how this noise aliases into the range from 1 to 5 kHz during the down sampling operation (when we down sample to 10 ksamples/s).

You will first try straight down sampling of the noise signal from 100 kHz to 10 kHz (just taking every 10th sample).  You will see that you get severe aliasing using this process. You will then design a digital FIR anti-aliasing filter as part of your down sampling system to minimize high frequency noise aliasing into the low frequency bands (< 5 kHz) that we care about. The process of digitally filtering the signal first and then down sampling (by taking every 10th sample of the filtered signal in this case) is called "decimation".

===== General Requirements =====
  * You must review aliasing
  * You must demonstrate aliasing of actual optical noise during the down sampling process when no digital anti-aliasing filter is used
  * You must design a digital FIR low pass anti-aliasing filter
  * You must analyze the performance of the FIR low pass filter and see how it reduces noise aliasing in the down sampled signal

===== General Notes =====
For the noise analysis, you will use a previously measured time-domain signal of optical noise in a lab room. Like you did in Milestone 1 Task 1, you will analyze a signal that is 500ms (half second) in duration, sampling at a rate of 100 ksamples/s, yielding a total of 50,000 samples. In this case, we simply want an optical noise signal without any player frequency signal on top of it. While you could use one of your 'noise only' signals acquired from last lab, I suggest you use the provided optical noise signal {{:light.zip|light.zip}}, since it is known to have sufficient noise to be detectable in this lab.

===== Specifications =====
Anti-aliasing FIR filter design specifications (we will discuss these more below):
  * Maximum filter length: 81
  * Maximum player variation: 1dB
  * Out of band rejection: 40 dB for frequencies between 10 and 50 kHz (also pay attention to attenuation in 5-10 kHz range, but doesn't need to meet the 40 dB rejection there)

Useful system specifications:
  * Stopband: 10kHz < f < 50 kHz  (with a corner frequency between 4.5kHz and 5.5kHz)
  * Player frequencies: 1471, 1724, 2000, 2273, 2632, 2941, 3333, 3571, 3846, 4167 (all in Hz)

===== Resources =====

==== Optical Noise Sample ====
  * Time-domain signal of optical noise: {{:light.zip|light.zip}}

After uncompressing the archive, use the MATLAB ''load'' command to read in the data.

<code>
load light % read optical noise data
</code>

==== Aliasing Background ====
To start off you should review aliasing with sampled signals. You covered this in ECEN 380. You can also read about it online on [[http://en.wikipedia.org/wiki/Aliasing|Wikipedia.]]

To improve your ability to use MATLAB you are going to start by looking at the problems associated with aliasing. The description of the aliasing problem, MATLAB plots, and sections of MATLAB code should be copied into your lab report.

As you know, it is hard to see the aliasing problem by looking at a time domain plot. However, in ECEN 380 you learned about looking at signals in the frequency domain. We are again going to use the built-in MATLAB function 'fft' for performing the fast Fourier transform (FFT) on your signal data.

[[matlab#using_the_fft_command |Click on this link to review how to use the MATLAB fft command.]]

Let's explore aliasing using some synthesized data. In MATLAB, you are going to generate two sine waves, one with frequency of f1 = 3 kHz and amplitude of A1 = 1 V, and one with frequency of f2 = 6 kHz and amplitude of A2 = 2 V. The MATLAB code to produce these signals (using a 100 kHz sampling rate) is: 

<code>
Fs = 100e3;  %sampling rate
f1 = 3e3;     %signal frequency
f2 = 6e3;
A1 = 1;       %signal amplitude, This is peak amplitude not peak-to-peak
A2 = 2;
% We're going to use a time axis from -0.25 to 0.25 seconds, which is 0.5s wide
% This yields Fs*0.5s total samples on our time axis
t = linspace(-0.25, 0.25, Fs*0.5);
x = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t);
</code>

Use the MATLAB fft command to produce the Fourier transform of your signal. Your signal should look like the image shown below (if you only plot the positive frequency axis after doing the FFT). Zoom in on the spikes and make sure that they are centered at 3 kHz and 6 kHz. If your signal does not look like this plot then you might want to go back and review the [[MATLAB |MATLAB fft command]].
{{ :sin_fft_2.jpg?300 |}}

Now you are going to down-sample the data to a sampling frequency of Fs = 10 kHz. (This is the sampling frequency that you designed your 10 bandpass filters for.) We will accomplish this by taking every tenth sample of the original signal. Here is MATLAB code to perform the down-sampling: 

<code>
x_ds = x(1:10:length(x));
t_ds = t(1:10:length(t));
Fs_ds = Fs/10;
</code>

Since the new sampling frequency is 10 kHz, the maximum frequency we can resolve is now reduced to 5 kHz. The 3 kHz sine wave is below 5 kHz, so we would expect that frequency to remain unchanged in the downsampled signal.  However, we would expect the 6 kHz signal to produce aliasing at 4 kHz with the new sampling rate of 10 kHz.  (Do you understand why?)

Use the MATLAB fft command to calculate the frequency domain plot of your new downsampled signal. You should produce a plot that looks something like this:

{{::alias_fig2_v2.png?200|}}

Again, you can see the formerly 6 kHz signal aliasing in to the new frequency range as a 4 kHz signal.

What would happen if the second sine wave originally had a frequency of 7 kHz instead of 6 kHz? (Answer: It would alias down to 3 kHz, on top of the original 3 kHz signal!)
==== Aliasing of Noise ====
The performance of the laser tag system can be improved by reducing the noise that ultimately ends up across the band of frequencies we care about (~1 - 5 kHz). We are not going to be investigating all of the noise sources in our system, but we are going to look into one of the easiest noise sources to reduce: the optical noise produced by the fluorescent lights in our lab room.

We measured the signal produced by the fluorescent lights and saved it as a MATLAB file. {{:light.zip|You can download the data here.}} This data is all noise. After unzipping this file you can load it directly into MATLAB. It contains two variables, t and y. (Use the MATLAB load command.)

You should get a plot similar to the following.

{{ :alias_fig3.jpg?200 |}}

As mentioned, this signal is all noise. It is hard to imagine the ability to detect a signal with an amplitude of ~10mV if the noise has an amplitude of ~1000 mV. That is the beauty of signal processing. To be able to understand this noise signal we need to plot it in the frequency domain.

As we've now done many times before, you can plot the signal in the frequency domain using the MATLAB function 'fft'. [[matlab#using_the_fft_command |Click here to learn how to use the MATLAB fft.]] As you know, in order to get the frequency axis right after performing the FFT you need to know the sampling frequency of the data.  If you didn't happen to know the sampling frequency of this data, you could extract it from the time vector by looking at the difference between subsequent time data points (which would be the sampling period, or 1/Fs).

<code>
Fs=1/max(diff(t))
</code>

Produce a frequency domain plot of the noise you loaded from the file provided. It should look something like the following:

{{ :alias_fig4.jpg?300 |}}

This plot shows that the noise in clumped around three main locations (1) low frequencies, (2) near 21 kHz, and near 42 kHz.

For the rest of this lab, you are going to use your own acquired optical noise data instead of the data that we provided you. Just like in Milestone 1, export 500ms of data sampled at 100 kHz (yielding 50,000 samples), but this time point your sensor upwards towards the fluorescent lights in the lab room. Plot the FFT of your new acquired optical noise data, with an appropriate frequency axis and paying attention to the amplitude scaling of your FFT data.

Zoom in on the plot and calculate the amplitude of the frequency components in the band where our player frequencies are located (1kHz < f < 5 kHz).

Now, as we discussed in lab lecture, this data will be acquired at a sampling frequency of 100 kHz by the ZYBO board. However, the 10 bandpass filters that we designed in Milestone 2 Task 1 are designed to process signals sampled at 10 kHz (since we only care about signals < 5 kHz and we want to keep our computations simple so our filtering can be done in real time on the ZYBO board). That means we need to down-sample the data to a 10 kHz sampling rate.

Let's first try getting to the new lower sampling rate by simply using every tenth sample like we did before when we reviewed aliasing. Down-sample your measured signal (which is 50,000 samples long, sampled at 100 kHz) to a 5,000 sample signal by taking every tenth sample.

Plot the down-sampled signal in both the time domain and frequency domain. It should look similar to the following (although you may have some more or less pronounced frequency peaks in your own data depending on how well your receiver board filtered some of the higher frequencies):

{{ :alias_fig5.jpg?300 |}}

Notice that the time domain signal looks very similar. However, the frequency domain signal now has a maximum frequency of 5 kHz (i.e., half of your new sampling frequency of 10 kHz). The noise in our desired frequency band (1kHz < f < 5kHz) is now drastically higher than it was prior to down-sampling. This increase in the noise is caused by aliasing of some of the higher frequency peaks in the original signal (in our noise data, there is a peak at f = 21 kHz and f = 42 kHz that alias into the band we are interested in).

So how can we avoid getting all of this higher-frequency optical noise aliasing into our 1 - 5 kHz band when we down-sample to 10 kHz?

The solution is to use a digital anti-aliasing (low-pass) filter before we down-sample the data. If we want to down-sample the data to a new sampling frequency of 10 kHz, our digital anti-aliasing filter would ideally have a cut-off frequency of 5 kHz. It would then filter out the higher frequency noise (like the 21 kHz and 42 kHz peaks we saw in our data) **before** we take every tenth sample, avoiding the aliasing problem. This process of first digitally filtering our signal (sampled at a high sampling rate) and then down-sampling the filtered signal is called 'decimation'.

So to summarize, what our system is going to do is (1) sample the incoming data at a sampling rate of 100 kHz, (2) filter our higher frequencies from this data using a digital anti-aliasing (low-pass) filter, (3) down-sample the filtered data to a sampling rate of 10 kHz, and (4) send the resulting signal through our bank of 10 filters to determine if there is a player hit.

==== Anti-Aliasing Filter Design====
Your next task is to design the digital anti-aliasing filter discussed above.

There is a very clever trick that we can use to reduce the computational complexity of this decimation step (digital filtering followed by down-sampling). If we use an FIR filter (refer back to the Lab #5 material from ECEN 380 for a review of FIR filter design), then the digital filtering step can be performed with a discrete finite convolution in the time domain. The convolution involves performing a sum **for each** of the samples in the output signal. However, after the digital filtering operation, we are going to throw away 9 out of every 10 samples in the resulting signal anyway. So why bother computing those 9 out of every 10 values that we are going to throw away?  If we only compute the 1 out of 10 values that we need, we will cut the operations needed to perform our filtering down by a factor of 10. See [[http://en.wikipedia.org/wiki/Decimation_(signal_processing)|decimating]] for more information.

When we actually implement decimation on our ZYBO board for real-time processing, that's exactly what we'll do. For now, however, let's design our FIR filter, determine the filter coefficients, and then simply use the 'filter' command in MATLAB, followed by down-sampling. This will be perfectly adequate for testing our signal processing algorithms before we implement them on the ZYBO board.

The specifications of the filter are provided at the top of the lab. You will be using the FIR filter design process that you used in ECEN 380 Lab #5 Task #2. For your convenience, here are links to the ECEN 380 Lab #5 and Dr. Rice's discrete filter design handout:

{{:lab_5_fall2016.pdf|ECEN 380 Lab #5}}
{{:filter_design_notes.pdf|Handout on discrete time filter design.}}

To refresh you memory, the basic process is:

  - Consider a **rect** function in the frequency domain with the appropriate width
  - Determine the resulting **sinc** function in the time domain
  - Multiply the time-domain **sinc** function times a windowing function (to limit the length of the impulse response, and to reduce the height of the side lobes in the frequency response of the filter)

The main design parameters of the filter are the following:
  - Filter length
  - Corner frequency
  - Filter windowing

==== Filter Characterization ====
Two MATLAB functions that are good for characterizing the filter are the following:
  - [[matlab#using_the_freqz_command|freqz: Frequency response of the digital filter]]
  - [[matlab#filter|filter: Function to filter a signal using filter coefficients]]

To understand how this works we are going to look at a filter without a windowing function. (In your design you will be adding in a windowing function.) Let's look at a low pass filter with a corner frequency of **f = 5 kHz**. In the frequency domain this would be a **rect** function. If you take the inverse Fourier transform of a **rect** function you get a **sinc** function. The length of the filter determines how much the filter looks like the ideal **sinc** function.

Let's start by looking at the filter in the time domain. (This means that we are looking at the impulse response of the filter.) Since we are working with data sampled at a frequency of 100 kHz, the maximum frequency we can resolve is Fs/2 = 50 kHz.

When you review your ECEN 380 code for designing digital filters, you may be confused by seeing frequencies described in a variety of different units. In this lab, we are dealing with time domain signals, and a natural way to think about frequency is in units of cycles/s, or Hertz. However, you may see frequencies in your FIR filter design code from ECEN 380 in units of cycles/sample (particularly if you took ECEN 380 prior to Fall 2016). These units are often used in digital filter design. To convert from frequency in Hz to frequency in cycles/sample, do the following:

  (Frequency in cycles/sample) = (Frequency in Hz) / (Sampling frequency)
  
You will further sometimes see frequencies in units of radians/sample. To convert from frequency in cycles/sample to frequency in radians/sample, do the following:

  (Frequency in radians/sample) = 2*pi*(Frequency in cycles/sample)
  
Finally, as you saw with the MATLAB 'butter' command, occasionally certain MATLAB commands will use other units, like half-cycles/sample (where a frequency of 1 corresponds to one half the sampling frequency).

The important thing to take away from all of this is to **be consistent in what units of frequency you are using**. Also, make sure to read the help for each MATLAB function that requires frequencies to ensure that you are passing in frequencies in the units that the MATLAB function is expecting. And remember that certain MATLAB functions (like 'freqz') have several different forms depending on what parameters you pass in, and that those different forms might require frequencies in different units. For example, here is the MATLAB help descriptions from two of the different forms of freqz:

<code>
H = freqz(B,A,W) returns the frequency response at frequencies
    designated in vector W, in radians/sample (normally between 0 and pi).
    W must be a vector with at least two elements.
    
H = freqz(B,A,F,Fs) returns the complex frequency response at the
    frequencies designated in vector F (in Hz), where Fs is the sampling
    frequency (in Hz).
</code>

The first form of freqz above takes the input frequencies in radians/sample, whereas the second form takes the input frequencies in Hz (and then requires a sampling frequency so that it can convert for you.)

With this in mind, a low pass filter with a corner frequency of 5 kHz that operates on data sampled at 100 kHz would be a **rect** function that goes from -0.05 to 0.05 in cycles/sample. (Can you see why? We want the low-pass filter to go from -5 kHz to 5 kHz, and the data was sampled at Fs = 100 kHz. So, in units of cycles/sample, our rect goes from -5kHz/100kHz = -0.05 to 5kHz/100kHz = 0.05 cycles/sample.)

IMPORTANT: If you are using the FIR filter design code from the Dr. Rice's filter design handout, note that the frequencies are done in units of cycles/sample. Also note that the code given produces a bandpass filter (so you'll need to set your center frequency to 0 for a lowpass filter), and the B parameter in that code is really the full width of the filter, not the half width. If you took ECEN 380 Fall 2016 or later, your code from Lab #5 Task #2 is likely written with frequencies in Hz.

Plot the resulting sinc function for both a long (N = 2001) and short (N = 21) filter length. Your resulting plots should look like the following:

{{ :filter_time.jpg?200 |}}
  
Notice how the longer filter looks a lot more like the sinc function than the shorter one.

Now we want to look at the resulting filters in the frequency domain. We are going to do this by using the **freqz** MATLAB function. The **freqz** MATLAB function is used for both FIR and IIR digital filters. Since we just designed an FIR filter, our impulse response (time domain) directly gives the **b** coefficients of our filter. The **a** coefficients are a0 = 1 and a1, a2, ... = 0. (If this isn't sounding familiar, review the FIR filter design section of the ECEN 380 Lab #5 handout.)
  
Use the **freqz** function to plot the frequency response of your N = 2001 and N = 21 length filters. Your plots should look like the following:

{{ :filter_freq.jpg?200 |}}

//Notice that the long filter looks a lot more like your desired low pass filter.//

Now let's look at how we specify the performance of a filter. This picture illustrates the various regions of a low pass filter:

{{ :lowpassfilterresponse.gif |}}

The ideal low pass filter has a flat frequency response in the passband region, complete attenuation in the stopband region, and an infinitely thin transition band.  

We have 1 specification related to the passband region and 1 related to the stop band region.  

**Specification 1 (We call this Player Variation.)**: This specification is similar to what is called Passband Ripple. Passband Ripple is the total variation in the frequency response for frequencies within the passband. However, our specification is a little bit looser than Passband Ripple. We just want to make sure that the filter has the same response for each player. So our specification is the total variation in the response for the 10 player frequencies.

**Specification 2 (This one is called Out of Band Rejection)**: This specification is the amount of attenuation experienced by any signal with a frequency component within the stopband. Since we are designing a low pass filter, the stopband is any frequency larger than the transition band.
=====Filter Characterization: (1) Player Variation (in band uniformity) =====

This specification is related to changes to the player signals. If one player signal is being attenuated more than another by the filter, then that player will have an advantage because the player shooting will need to be closer. Therefore, a system with a lot of player variation would be a poorly designed game.  

To determine Player Variation, we plug each of the player frequencies into the transfer function of our low-pass FIR filter and see what the gain is at each of the 10 player frequencies. We then take the ratio of the smallest value to the largest value (typically reported in decibels, so 20*log10(ratio of smallest player gain to largest player gain)).

[[MATLAB|In the investigation of the freqz command]] we learned that we can pass specific frequencies to the freqz command. This is the easiest way to find out what the gain is at each of the 10 player frequencies. You just pass freqz the 10 player frequencies, and it will return the gain at each frequency.

  * For the long (N = 2001) filter
    * The values are: 1.0005, 0.9999, 0.9976, 1.0003, 1.0012, 1.0026, 1.0018, 1.0010, 1.0048, 1.0031
    * This corresponds to a variation of 0.9976/1.0048=0.9928. As mentioned, we'll typically report this in decibels instead. So this variation becomes 20*log10(.9928)=-0.06 dB, which meets the specification because -0.06 > -1.
  * For the short (N = 21) filter
    * The values are: 1.0938, 1.0638, 1.0265, 0.9854, 0.9256, 0.8698, 0.7945, 0.7469, 0.6906, 0.6241.  
    * This corresponds to a variation of 0.6241/1.0938=0.62. So this variation becomes 20*log10(.62)=-4.87, which does not meet the specification because -4.3 < -1.

=====Filter Characterization: (2) Out of Band Rejection =====

The objective of the anti-aliasing filter is to attenuate any signals with frequency components above about 5 kHz (i.e., half the final sampling rate). The distance between the highest player frequency (4167 Hz) and the maximum frequency (Fs/2 = 5 kHz) is 833 Hz. This would make the transition band 4167Hz < f < 5833Hz.

This transition band is fairly small, which would make the design of the low pass filter difficult. Therefore, we are going to widen the transition band to make the design easier. We are going to use a transition band of 4.5kHz < f < 10kHz.

This wider transition band means that some noise within the transition band can alias into our passband. However, we fortunately don't have much optical noise in these frequency bands, as shown here: 
{{ :alias_fig6.jpg?200 |}}

You can see that the noise is fairly low in the transition band, so our looser specification should not result in too much additional noise aliasing, and shouldn't have too much impact on our final signal.

To evaluate this specification (Out of Band Rejection), use the freqz command again, but this time pass in a frequency array with values from 10 kHz up to 50 kHz. This will return the frequency response across the range of frequencies in your stop band. The maximum gain (absolute value of your frequency response) across this range of frequencies is your Out of Band Rejection. Similar to the Player Variation specification, you should also report this maximum gain in units of decibels. However, unlike Player Variation, in this case we want the number to be as small as possible.

For the long filter length, Out of Band Rejection should be around -58 dB, and for the short one it should be around -25 dB.

===== Actual Filter Design =====

You have now done a simple filter design, plotted the response of the filter in the time domain and frequency domain, and determined the specifications of a filter. You are now ready to design and characterize the low pass FIR filter that you are actually going to use.

In the sections above, we were able to get an FIR filter to meet all of the performance specifications. However, the filter is too long. With a filter length of 2001 you will use up too much of the processing power of the ZYBO board.

**We are limiting you to an FIR anti-aliasing filter length of 81.** To meet the performance specs with this length of a filter, you will need to optimize the windowing function, and possibly use a corner frequency of higher than 5 kHz.

As a review here are the design parameters.
  * Filter length (filter performance increases with length so the best performance will probably be with a length of 81)
  * Corner frequency
  * Windowing function 

Happy filter designing, and good luck! :-)

==== Noise Reduction ====

Now that you have finished the design of the anti-aliasing filter, your final job is to see how well it reduces the aliased noise.  Do the following:
  - Load in your optical noise sample (the one that is 500ms long with 50,000 samples)
  - Plot the noise spectrum (frequency domain) **without using your FIR filter** (just like you did before)
    - Down-sample to Fs = 10kHz
    - Plot the noise from 1kHz to 5 kHz
  - Plot the noise spectrum (frequency domain) **using your FIR filter before down-sampling**
    - Filter the noise signal with your FIR filter
    - Down-sample the signal to Fs = 10kHz
    - Plot the noise from 1kHz to 5 kHz
  - Compare the two signals

Hopefully you are getting pretty comfortable with this kind of MATLAB coding at this point, but here is some example code if you need some help:

<code>
load light % read optical noise data

% down-sample
y1 = y(1:10:length(y));
t1 = t(1:10:length(t));
Fs = 1/max(diff(t1)); % sampling frequency, which should be 10 kHz

% calculate frequency spectrum
Y1 = fft(y1);
len = length(Y1);
Y1 = Y1(1:floor((len)/2));
Y1 = 2*Y1/len;
freq = linspace(0, Fs/2, length(Y1));

b=[  % copy your filter coefficients here
];

y2 = filter(b,1,y); % filter the optical noise using the FIR filter

% down-sample
y3 = y2(1:10:length(y2));

% frequency spectrum
Y3 = fft(y3);
len = length(Y3);
Y3 = Y3(1:floor((len)/2));
Y3 = 2*Y3/len;

% plot the 2 noise spectra
subplot(2,1,1);
plot(freq*1e-3,abs(Y1));
ylim([0 .2]);
subplot(2,1,2);
plot(freq*1e-3,abs(Y3));
ylim([0 .2]);

</code>

===== What is Needed in the Lab Report =====
  - Description demonstrating an understanding of aliasing (including MATLAB plots that you generated)
  - Description of your FIR low-pass filter design including
    - Description of your design including bandwidth and windowing details
    - Filter coefficients
    - Frequency domain plot (i.e., transfer function of the filter)
    - Filter specifications (in-band uniformity, out of band rejection, filter length)
  - Optical Noise
    - Plot of optical noise down-sampled to Fs = 10 kHz (just taking every 10th sample)
    - Plot of optical noise decimated to Fs = 10 kHz (low pass filtered + down-sampled)

===== Pass Off =====
The following items need to be shown to the TAs for pass off:
  - MATLAB code portion showing design of FIR low pass anti-aliasing filter
  - List of FIR filter coefficients
  - List of achieved filter specifications
  - Frequency response plot of the low pass FIR filter
  - Plot of optical noise doing a straight down-sampling to 10 kHz (with no FIR filter)
  - Plot of optical noise decimated to 10 kHz (low pass filter + down-sampled)
ECEN 390 - Laser Tag Project

User Tools

Site Tools

Page Tools
