Milestone 3, Task 1: Implement the FIR and IIR Filters

Now that you have designed your FIR and IIR filters using MATLAB, you will implement these filters using C-code that executes on the the ZYBO board. We will be using the queues that you designed previously to ease the overall implementation of these filters. You will verify the correct operation of your filters using provided test code in the ecen390 project directory.

Ultimately, for all of Milestone 3, you are writing the software - the detector - the thing that will detect when and what player frequency “hits” you when playing the game. For Task 1, you are implementing just the filtering and power-computation parts. The diagram below shows you the entire structure for the detector. You are implementing the part contained in the gray box. This consists of the following parts:

• xQueue: this contains the input data for the filters. Ultimately, inputs to the xQueue will be based upon digitized voltages obtained from the photo-diode circuitry that you constructed in ECEn 340. We will be using test data (provided) for this task.
• FIR Filter: the decimating anti-aliasing FIR filter.
• yQueue: the output from the FIR filter and the input to the 10 IIR bandpass filters.
• IIR Filter 0 - 9: the 10 IIR band-pass filters.
• zQueue 0 - 9: these are part of the IIR filters.
• outputQueue 0 - 9: these are the outputs of the 10 IIR filters.
• Compute Power 0-9: compute the amount of power in the signal output by the corresponding IIR filter.

Note: The IIR filters and associated queues are now numbered 0 - 9 in this task. This is because 'C' arrays are zero indexed and starting at 0 makes coding easier. In MATLAB for Milestone 2, the frequencies and IIR filters were numbered 1 - 10.

• You must implement all of the functions that are listed in filter.h.
• You must use the queue-based approach (described in the Background section) to implement your filters.
• You must demonstrate the behavior of your filters to the TAs, using the provided source code.
• Your FIR-filter must have a reasonably flat frequency response up through 5 kHz and then should quickly roll off.
• Your IIR-filters must have a very narrow frequency response. For example, for frequency 4, your power output should be high for IIR-filter 4, but should be very low for all of the other IIR-filters. A difference of 10X is preferred for good performance in the laser-tag game.
• You will test your filters and functions that compute power using the provided filter_runTest() function.

• Assume a sample-rate of 100 kHz for the input to the FIR filter.
• The decimation rate for the FIR-filter is 10.

1. A working, bug-free queue. If your queue code passed the provided tests from Milestone 1, you should be OK.
2. All necessary coefficients for the FIR and IIR filters. You will need one set of coefficients for the FIR filter and 10 sets of coefficients for the IIR filter code, e.g., for the 10 band-pass filters. Each of these coefficient sets consists of the a and b coefficients for the IIR filters.

1. Implement all of the functions shown in filter.h.
2. Use queues as described in this task. You must have an xQueue, a yQueue, an array containing 10 zQueues, and an array containing 10 outputQueues. You must declare these queues as static in filter.c.
3. Declare a 1-D array for the FIR coefficients, for example, const static double fir_coeffs[FIR_COEF_COUNT] = {0.25, 0.5, 0.75, 1.0}. See the examples in the code below.
4. Declare two 2-D arrays for the IIR filters, one for the A coefficients and one for the B coefficients, similar to how the FIR coefficients were declared, for example: const static double iir_a_coeffs[FREQUENCY_COUNT][IIR_A_COEFF_COUNT] = {{0.25. 0.5, 0.75, 1.0}, {1.0, 2.0, 3.0, 4.0}, …};
5. You must have a separate init function for each queue. For example, let's assume that the function that initializes my xQueue is called initXQueue(). Inside this function you would call queue_init() on xQueue and then would use queue_overwritePush() or queue_push() with a for-loop to fill the xQueue with zeros. You must write a corresponding function for each of your queues and you will call these functions in filter_init().
6. You must implement your filters as follows. filter_addNewInput(value) must push a value onto xQueue. filter_firFilter() reads input values from xQueue using queue_readElementAt() and pushes its output values onto yQueue. filter_iirFilter(filterNumber) reads input values from yQueue using queue_readElementAt() and pushes its output values onto zQueue[filterNumber].
7. Implement all of the power-computation functions using the comments to guide you.
8. You must test all of your filter.c code using the filter test code provided below. More information is provided below so you know what to expect when running the tests.
9. You must follow the software coding standard. Exception: clang-format is not required. Ignore any clang-format rules.

• Convert MATLAB Filter Coefficients to C Code
• Frequency Response Plots
• Video Demonstration of FilterTest code

Note that the following files are provided in your ecen390 project directory. The test code is used to check the correctness of your code.

• lasertag/main.c
• lasertag/filter.h
• lasertag/support/filterTest.h
• lasertag/support/filterTest.c
• lasertag/support/histogram.h
• lasertag/support/histogram.c

You are expected to create and implement the following file. See the provided header file (.h) for a description of each function.

• lasertag/filter.c

Declare and initialize 22 queues: an xQueue to store incoming inputs from the ADC, a yQueue that holds the history of outputs from the FIR filter, an array of 10 zQueues that hold the history of outputs from the corresponding IIR filters, and an array of 10 outputQueues that accumulate output values from each of the IIR filters for the purpose of power computations. Declare all of these queues as static queue_t variables in filter.c. Give each of these queues a meaningful name using the name argument in queue_init() function.

Here's an example of how to declare and initialize an array of queues. Note this code will not compile “as is.”

#define QUEUE_INIT_VALUE 0.0
#define Z_QUEUE_SIZE IIR_A_COEFFICIENT_COUNT
static queue_t zQueue[FILTER_IIR_FILTER_COUNT];	
 
void initZQueues() {
  for (uint32_t i = 0; i < FILTER_IIR_FILTER_COUNT; i++) {
    queue_init(&(zQueue[i]), Z_QUEUE_SIZE, "zQueue");
    for (uint32_t j = 0; j < Z_QUEUE_SIZE; j++)
     queue_overwritePush(&(zQueue[i]), QUEUE_INIT_VALUE);
  }
}

In filter_init(), initialize all of the queue's and fill them with zeros. Below is one example, that uses small helper functions to initialize each queue and fill it with zeros. These helper functions are then called from within filter_init().

// This must be called before invoking any filter functions.
void filter_init() {
  // Init queues and fill them with 0s.
  initXQueue();  // Call queue_init() on xQueue and fill it with zeros.
  initYQueue();  // Call queue_init() on yQueue and fill it with zeros.
  initZQueues(); // Call queue_init() on all of the zQueues and fill each z queue with zeros.
  initOutputQueues();  // Call queue_init() on all of the outputQueues and fill each outputQueue with zeros.
  ...

Note: For queues that are used to store signal histories, always fill them completely with zeros. This provides a stable starting condition.

Note: an empty queue is not the same as a queue that is filled with 0s.

At run-time, you “run” the filters by doing the following:

• As they become available, add new inputs to the xQueue. This is what the function filter_addNewInput() will do (once you implement it). When you finish the completion of your laser-tag system, you will use filter_addNewInput() to take the output from the digitized voltage from your photo-diode-based detector and provide it to the FIR filter.
• Each time you receive 10 new inputs, run the filters as follows:
  • invoke filter_firFilter(). This will use the appropriate FIR coefficients and the xQueue that contains the necessary signal history to compute the output of the FIR filter. Add this output computed by filter_firFilter() to the yQueue (filter_firFilter() also returns this computed value). The yQueue then contains the output from the FIR filter; this is the input for each of the 10 IIR-filters.
  • invoke the filter_iirFilter() function 10 times, once for each frequency. Return the computed z-value for each frequency. Also add the output from the IIR-filter to the corresponding outputQueue. You will use the values contained in the outputQueues to compute the total power of the signal that passes through the IIR band-pass filter.

For the IIR filters, note that the first a-coefficient (a0) is always 1.0 and is not used in the computation. As such, if your a-coefficient array is of size 11, your z-queue must be size 10 (1 less than the size of the array) because a0 is always ignored. Watch out for this!

In this task you will implement the FIR and IIR filters that you designed using MATLAB using 'C' code. Let's start with the FIR filter. Remember that the FIR filter is implemented as a weighted sum of some past number of inputs. Here's an example from Wikipedia:

It can be confusing to transition from the finite array-based approach used in MATLAB to the “infinite” approach that is required in the implementation of a signal-processing system. The inputs and outputs of a real-time signal-processing system are essentially infinite. As such, the array-based notation in the equation above fails us because the output is an indexed array y[n]. For example, at time=0, you start out computing y[0]. After playing the game for several minutes, n would be in the billions. And, it only goes up from there. Simply put, you want to eliminate the [n] part so that the output is simply y.

Note that when we read the English-based description from Wikipedia (see above), indexes were not discussed. Remember that the FIR-filter is implemented as a “weighted-sum of some past number of inputs”. All those indexes, the i, the k, etc., are just a way to keep the coefficients properly aligned with the data. As long as we can keep the incoming inputs properly aligned with the coefficients, we are good to go.

The idea is pretty simple and is based upon these ideas:

1. Create a data structure that will keep an ordered history of past values. The size of the data structure must match the order of the filter, e.g., a 50-tap FIR filter needs a history of 50 values.
2. At start-up time, fill the data structure with zeros.
3. As each new value arrives, throw away the oldest value.
4. Read the stored past values from the data structure and multiply them with the correct coefficients.

As you have probably guessed at this point, the queues that you implemented as a part of Milestone 1 are the perfect data structure for this purpose.

You can implement a FIR filter using the queues that you have already coded. Consider an example where the FIR filter uses 4 past values to compute its output. In the example code below, I have “pushed” four values onto the queue. Assume that these are 4 values that are based on values from the ZYBO's ADC. Note that for pedagogical purposes, the code below does not necessarily adhere to the coding standard.

#include "queue.h"
#include <stdio.h>
#define FIR_COEF_COUNT 4
 
int main() {
  // Initialization stuff.
  queue_t xQ;                           // x is the queue with the input history for the queue.
  queue_init(&xQ, FIR_COEF_COUNT);      // Size of history queue must equal coefficient count.
  for (uint32_t i=0; i<FIR_COEF_COUNT; i++)  // Start out with a queue full of zeros.
    queue_overwritePush(&xQ, 0.0);
 
  const double b[FIR_COEF_COUNT] = {0.25, 0.5, 0.75, 1.0};  // These coefficients are used for this example.
 
  // Add some example inputs to the queue.
  queue_overwritePush(&xQ, -0.1);  // Add a new input to the queue (oldest in input history).
  queue_overwritePush(&xQ, -0.4);  // Add a new input to the queue.
  queue_overwritePush(&xQ, 0.24);  // Add a new input to the queue.
  queue_overwritePush(&xQ,  0.54); // Add a new input to the queue (newest in input history).
 
  // Compute output of FIR-filter (y)
  // using a single lone statement (broken into 4 lines to keep it readable).
  // This is just for example. You will use a for-loop as shown below.
  double y;
  y = queue_readElementAt(&xQ, 0) * b[3] +
      queue_readElementAt(&xQ, 1) * b[2] +
      queue_readElementAt(&xQ, 2) * b[1] +
      queue_readElementAt(&xQ, 3) * b[0];
  printf("%lf\n\r", y);
 
  // Add new input.
  queue_overwritePush(&xQ, 0.33);
 
  // Compute the next y using a for loop (a better way).
  // += accumulates the result during the for-loop. Must start out with y = 0.
  y = 0.0;
  // This for-loop performs the identical computation to that shown above.
  for (uint32_t i=0; i<FIR_COEF_COUNT; i++) { // iteratively adds the (b * input) products.
    y += queue_readElementAt(&xQ, FIR_COEF_COUNT-1-i) * b[i];
  }
  printf("%lf\n", y);
 
}

Pictorially, implementing the FIR filter would appear as shown below. As shown, you can see the 4 values that were pushed onto the queue as well as the total computations.

The next computation of y is shown below. You can see that by adding a new value to the queue, all of the other values shifted over, relative to the b coefficients. Thus you can use the same code to compute y over and over again.

You can see that the purpose of the queue is to store past values in the order that they were received and make all of the queue-contained values accessible during the computation.

Decimation is really easy. In our laser-tag system we will be decimating by 10. All we do is invoke our FIR-filter each time we receive 10 new samples. As you add incoming samples to the FIR-filter input-queue, only invoke the FIR-filter each time you have received 10 new inputs. You will then invoke the IIR filters right after you invoke the FIR-filter. Decimation-wise, there is nothing required for this task. You will implement this later.

The equation for an IIR filter is shown below.

However, we can simplify this a bit because the first a coefficient is 1.0 and can be ignored.

We finally end up with:

The implementation of the IIR filter is similar to the FIR filter. However, the IIR filter relies on two signal histories: y and z, as shown in the equation above. As you can see from the equation, you would need two queues of different sizes (11 and 10) to keep the necessary signal histories. The only other difference is that the computed value (z) is also pushed onto the queue that keeps a history of z values. This is essentially what puts the “IIR” in the filter, e.g, feedback.

Implement all of the power-related functions (they all have the word “power” in their names). You will need to make sure to write filter_computePower() so that it does not take too much execution time. Carefully think about how you might be able to reuse computations performed in a previous invocation of filter_computePower() to reduce overall computation time. To initially debug your power code, you can fill the output queues with constant values, say 2.0 for example, and then compare the output from the function with your own calculation. The provided test code contained in the file filterTest.c provides a comprehensive test of the power functions.

To compute power, you must keep a running history of 200 ms of output data from each of the 10 IIR-based band-pass filters. To achieve this, do the following:

• Declare a static array of 10 power-output queues in filter.c (these were the output-queues that were discussed above). These will be indexed by filterNumber in the various functions that need to access them.
• Determine the size of the output queue so that it will hold 200 ms of output data.
• Initialize these queues and fill them with “0”s. Use a separate function for this initialization as was required in Task 1.
• Each time you invoke an IIR filter, push its output onto its respective output queue.
• Declare a static array of 10 doubles (in filter.c) that will contain the latest computed-power numbers for all 10 frequencies. Call this array currentPowerValue[].
• When you invoke filter_computePower(filterNumber), compute the power for a frequency by computing a sum of the squares of all values contained in the output queue for that frequency. Store this value in the currentPowerValue[] array at the index that corresponds to that frequency.
• Write the function filter_getCurrentPowerValue(uint16_t filterNumber) that simply returns the value stored in the array currentPowerValue[] at filterNumber.
• Write the function filter_getNormalizedPowerValues(double normalizedArray[], uint16_t* indexOfMaxValue). This function does two things:
  1. Normalizes the power values that are contained in currentPowerValue[] so that they are distributed between 0.0 and 1.0 and stores these normalized values in the array argument normalizedArray[]. To compute normalized values, simply find the largest value in the currentPowerValue[] array and divide all values stored in the array by the largest value.
  2. Stores the index of the currentPowerValue[] that contains the highest power value in the argument indexOfMaxValue.

You will need to make sure to write filter_computePower() so that it does not take too much execution time. Carefully think about how you might be able to reuse computations performed in a previous invocation of filter_computePower() to reduce overall computation time.

Verify the correct operation of filter_computePower(), filter_getCurrentPowerValue(), and filter_getNormalizedPowerValues() using your own test code. An easy way to do this is to fill the power queues with constant values, say 2.0 for example, and then compare the output from the function with your own calculation.

The code below provides context on how the filter functions will be used in a future task. Assume that this code is called whenever there is a new scaled ADC value available.

You don't need to write this code yet. For now, just enable the provided filter test code to verify that your filter functions work. The test code will call your filter functions with meaningful arguments.

// Constants and filter functions are declared in filter.h
 
...
filter_addNewInput(scaledAdcValue); // Add scaled ADC value to x-queue
sample_cnt++; // Count samples since last filter run
 
// Run filters and hit detection if decimation factor reached
if (sample_cnt == FILTER_FIR_DECIMATION_FACTOR) {
  uint16_t filterNumber;
  sample_cnt = 0; // Reset the sample count.
  filter_firFilter(); // Runs the FIR filter, output goes in the y-queue.
  // Run all the IIR filters and compute power in each of the output queues.
  for (filterNumber = 0; filterNumber < FILTER_FREQUENCY_COUNT; filterNumber++) {
    filter_iirFilter(filterNumber); // Run each of the IIR filters.
    // Compute the power for each of the filters, at lowest computational cost.
    // 1st false means do not compute from scratch.
    // 2nd false means no debug prints.
    filter_computePower(filterNumber, false, false);
  }
  ...
}
...

To pass off this task, you must run your filter and power code with the provided filterTest.c source code. The test code tests the arithmetic for all of your filters. It also plots the frequency response for the FIR and IIR filters on the TFT display. The provided test code comes in four pieces: histogram.h, histogram.c and filterTest.h and filterTest.c. In filter.h, you will see a section of code labeled “Verification-Assisting Functions”. These accessor functions provide a way for the test-code to access the various named variables in your filter.c code. These accessors are necessary because your naming schemes will likely differ from my test-code.

The filterTest code is a little over 1120 lines and performs several tests:

• checks to see that the FIR coefficients are properly aligned with the contents of the xQueue,
• checks to see that the FIR outputs are computed correctly,
• checks to see that the IIR B-coefficients are properly aligned with the yQueue,
• checks to see that the IIR A-coefficients are properly aligned with the zQueue,
• plots the frequency response of the FIR filter to the TFT display for frequencies ranging from 1.11 kHz to 50 kHz, and
• plots the frequency response of each IIR filter against each of the 10 user frequencies.
• checks to see that your power-computing functions generate correct values.

Note that “aligned” means that coefficient values are multiplied with corresponding queue values using correct indices.

To run this test code, uncomment filter_runTest() in the body of your main() code and uncomment #define RUNNING_MODE_TESTS. Along with various informational and perhaps error messages, the frequency response will be plotted out on the TFT display on the ZYBO carrier board. You can compare your results to those on this page. Your results should look similar. The informational messages that should appear in your console will look like this if everything passes. Numerical values won't be exact because your FIR filter will be different, but values should be roughly similar:

filter_runFirAlignmentTest passed.
filter_runFirArithmeticTest passed.
filter_runIirAAlignmentTest passed.
filter_runIirBAlignmentTest passed.
===== Starting filter_runPowerTest() =====
Testing to see that the power is computed correctly when forced.
Output queues are the correct size.
Power values were properly computed when forced.
Testing to see that the power is computed correctly incrementally over 3000 trials.
Power values were properly computed incrementally.
+++++ Exiting filter_runPowerTest +++++
running filter_runFirPowerTest() - plotting power values (frequency response) for frequencies 1.47 kHz to 50.00 kHz for FIR filter to TFT display.
freqCount:0, testPeriodPowerValue:1.760509e+03
freqCount:1, testPeriodPowerValue:1.707258e+03
freqCount:2, testPeriodPowerValue:1.655931e+03
freqCount:3, testPeriodPowerValue:1.628428e+03
freqCount:4, testPeriodPowerValue:1.609116e+03
freqCount:5, testPeriodPowerValue:1.588978e+03
freqCount:6, testPeriodPowerValue:1.538183e+03
freqCount:7, testPeriodPowerValue:1.491542e+03
freqCount:8, testPeriodPowerValue:1.417058e+03
freqCount:9, testPeriodPowerValue:1.298423e+03
freqCount:10, testPeriodPowerValue:1.119037e+03
freqCount:11, testPeriodPowerValue:3.572475e+02
freqCount:12, testPeriodPowerValue:5.366949e+02
freqCount:13, testPeriodPowerValue:2.239070e+02
freqCount:14, testPeriodPowerValue:4.032896e+01
freqCount:15, testPeriodPowerValue:1.090589e+00
freqCount:16, testPeriodPowerValue:3.327529e-03
freqCount:17, testPeriodPowerValue:1.363178e-03
freqCount:18, testPeriodPowerValue:1.627174e-03
freqCount:19, testPeriodPowerValue:3.659948e-03
freqCount:20, testPeriodPowerValue:3.264494e-03
Plotting response to square-wave input.
running filter_runFirPowerTest(0) - plotting power for all player frequencies for IIR filter(0) to TFT display.
running filter_runFirPowerTest(1) - plotting power for all player frequencies for IIR filter(1) to TFT display.
running filter_runFirPowerTest(2) - plotting power for all player frequencies for IIR filter(2) to TFT display.
running filter_runFirPowerTest(3) - plotting power for all player frequencies for IIR filter(3) to TFT display.
running filter_runFirPowerTest(4) - plotting power for all player frequencies for IIR filter(4) to TFT display.
running filter_runFirPowerTest(5) - plotting power for all player frequencies for IIR filter(5) to TFT display.
running filter_runFirPowerTest(6) - plotting power for all player frequencies for IIR filter(6) to TFT display.
running filter_runFirPowerTest(7) - plotting power for all player frequencies for IIR filter(7) to TFT display.
running filter_runFirPowerTest(8) - plotting power for all player frequencies for IIR filter(8) to TFT display.
running filter_runFirPowerTest(9) - plotting power for all player frequencies for IIR filter(9) to TFT display.

• You will show the TAs how your code behaves when running the provided test code. You will mostly show the TAs the information displayed on the LCD on the board.
• You will submit your source code by doing the following:
  1. From the top-level directory (e.g. ecen390), run the check_and_zip program: ”./check_and_zip.py 390m3-1”.
  2. The resulting .zip file will be in the top-level directory. Submit that to Learning Suite.
  3. Submit only one .zip file per team. Both team members will receive credit.

Please pay attention to the following:

1. Check to make sure that the filters pass all tests.
2. Check to make sure that the plots on the TFT display look correct, e.g., the FIR-filter is flat across the frequency range and that the bandpass filters have a narrow response.