Overview

Implementing the narrow-band LPF using an IIR filter is probably the most difficult option to design but the most straightforward to implement. One reason for the design difficulty stems from the fact that in order to get such a sharp response poles close to the unit circle are needed. These poles can then drift outside the unit circle and the system can then become unstable when finite precision effects are added. Also, perfectly linear phase (constant group delay) cannot be realized using IIR filtering.

There are several approaches to designing an approximately linear phase IIR filter. For example, an IIR filter could be run on blocks of samples in both the forward and reverse direction and the results of each block added together; using the filter in both directions would cancel the nonlinear phase response of the filter. Also, iterative design methods exist to design filters that simultaneously minimize errors in magnitude and group delay. Yet another approach is to design an IIR filter which approximates the desired magnitude response (e.g., an elliptic filter using the ellip command in MATLAB) and then design an IIR all-pass filter which compensates for the nonlinear phase. This last approach is the one we will examine here.

MATLAB Filter Design Toolbox

The MATLAB Filter Design (FD) Toolbox contains algorithms used for optimal or near optimal design of filters subject to various constraints. You can view a description of the toolbox and the functions it contains at this link. The FD Toolbox will be installed on the white Dell machine nearest the door in the ECE 320 lab.

Using the Filter Design Toolbox

Although it is possible to design a very good LPF magnitude response using an elliptic filter, we do not have the advantage with the ellip command of being able to constrain the poles away from the unit circle to prevent instability. Fortunately, the FD Toolbox provides a command called iirlpnormc which allows us to keep the poles within a circle of a specified radius. Note that this command implements a least-pp'th algorithm. The term least-pp'th signifies that the algorithm attempts to minimize a LpLp-norm error. In the case of the magnitude response, the LpLp-norm error is given as

where HH is the actual frequency response, HdHd is the desired response, and WW is some weighting function. Most of the time the weighting function used is one which equals 1 over the passband and stop-band and 0 in the transition band. The role of WW in the above equation is similar to that used in FIR filter design (remez). The relative weights in the stop-band and pass-band are set by WW and control the relative magnitude of the ripples in these bands. Note that minimizing the L2L2-norm is equivalent to minimizing the RMS error in the magnitude. In contrast, the L∞L∞-norm is equivalent to minimizing the maximum error over the frequencies of interest (why?). In this lab we are concerned with minimizing the L∞L∞-norm. Of course, we cannot use infinity in any of our computations so using a large number (e.g. 128) must suffice.

Once our magnitude response has been selected, we need to perform group delay equalization to yield approximately constant group delay. This can be done using the iirgrpdelay command in the FD Toolbox. Again, note that a least-ppth algorithm is used and that we can constrain the radius of the poles. The resulting all-pass filter can be connected in series with the nonlinear phase low-pass filter created with iirlpnormc to complete the entire system.

The FD Toolbox can also aid in analyzing quantization effects. We suggest using FDATool, a convenient GUI for interacting with the FD Toolbox, to carry out the analysis. See the internet documentation for more information on the functions available. Of course, you may choose to do this by scaling and rounding as you have done in previous labs. Note that even though MATLAB uses high-precision arithmetic you may find that for long IIR filters MATLAB has difficulty rendering frequency responses, etc. Thus, you may find it useful to design a filter that has half the passband ripple, half the stop-band attenuation, etc. and implement it twice in your code to meet the specification.

Note that FDAtool and the filter design toolbox (qfilt function) can be used to analyze quantization effects on various filter structures, as well as on the FFT. The quantization parameters can be chosen and optimized in FDAtool. Also, FDAtool (with or without the filter design toolbox) can compute correct scaling to avoid overflow.

Implementation Structures

There are several ways to implement an IIR filter. One of these, a cascade of second-order systems, we have already seen. An alternative is placing these second-order sections in parallel. Another common implementation is a lattice structure (see any standard DSP textbook), which tends to be more resistant to finite word-length effects and may be more computationally efficient. To examine your choices, as a starting point you should examine the MATLAB functions listed as "Linear System Transformations" when you type help signal at the command prompt (does not require the FD Toolbox).

One of the difficult aspects of an IIR lattice is that although the lattice coefficients are in the interval (-1,1), the internals of the lattice can grow to be prohibitively large. To compensate for this, an m-file (latcfiltn.m) has been created that performs normalization after each lattice section to prevent overflow. If you are interested in exploring a lattice implementation you may want to copy this m-file to your own directory and modify it to suit your needs. Note that there are comments within the file to indicate where you might add checks for overflow conditions.

We suggest the use of FDAtool and dfilt for structure transformations. The function dfilt works also without the Filter Design Toolbox. It is also useful for evaluating cascade or parallel connections of sub-filters. The MATLAB command fvtool can be used to quickly evaluate frequency response of various filter structures.

Many extremely efficient structures for IIR filter implementations exist. Two of special note are the following:

Questions

Simulation

Simulate the IIR system in MATLAB. Compute the response of the system with appropriate test inputs. Make sure to include side-effects due to finite precision in your simulation.

Reading Exercise

Read through the following resources:

Design Exercises

Given the filter specification given in the filter specification, answer the following questions:

Matlab Simulation

Using your results from the previous part, simulate the two-stage multi-rate filter in MATLAB. Plot the frequency response of each stage's filter using freqz and determine the overall frequency response of your multi-rate system to verify that it meets the specifications. Since there is not a command for directly finding the frequency response plot of a multi-rate system in MATLAB, you will have to be a bit creative.

Additional Questions (optional, but for your benefit)

Questions

Be prepared to show all the necessary plots and MATLAB simulations as well as answers to all of the questions posed above to your T.A. as your prelab.

Implementation Issues

Due to the limitations of the core file, it is not possible to take in more than 64 input samples from the A/D converter at a time (unless the core file is rewritten to accomplish this task). Therefore, when implementing a Fourier-based filter, you should use the C skeleton from Lab 4 to perform the FFT on a large block of samples. All of your filtering operations (i.e., the multiplications of DFT samples, the additions of the overlap, the discarding of samples) and function calls must be performed in assembly. You will be graded on the number of cycles per input sample based on the portion of code in your assembly routine.

You should use the fft.asm routine provided in Lab 4 to perform the forward and reverse FFT's. You should study this file to determine how it works. If you need to change the length of the FFT, you will first need to change the relevant parameters in your assembly file (i.e., N, K_FFT_SIZE, K_LOGN, and other variables). You will also need to change the following parameters in the FFT file:

          K_TWID_TBL_SIZE
          K_TWID_IDX_3
        

K_TWID_TBL_SIZE is the size of the twiddle tables (how long should these be for a given FFT length?) and K_TWID_IDX_3 is the amount by which the program increments through the twiddle table during at the third stage of the FFT. What is this increment for a given N? Is fft.asm a decimation in time or decimation in frequency algorithm?

You will also need a modified twiddle table when you change the length of the FFT to use fft.asm as written. For a length 1024 FFT, the twiddle tables are length 512 each. TWIDDLE1 is a table of sine values from zero to ππ, and TWIDDLE2 is a table of cosine values from π 2 π 2 to 3π 2 3π 2 . The support for the cosine and sine is different because fft.asm code uses the fact that sin-θ=-sinθ θ θ when performing computations. If you want a length 64 FFT, you will need to ``decimate'' the twiddle table to length 32, or in other words, only keep one out of every 16 lines in the twiddle tables and discard the rest. We will provide a MATLAB function, edit_twiddle.m for this purpose. The function call in this example would be:

          edit_twiddle('TWIDDLE1','new_twiddle1',16)
        

You should verify that the new twiddle tables you generate indeed have 32 elements. To perform an inverse FFT, you can use the standard FFT algorithm and then appropriately scale and shift the outputs. Lecture 43 of the ECE 310 notes on the Discrete Fourier Transform suggests how this may be done (Property 3 of ``Properties of the DFT'').