Cornell University ECE4760
IIR filter designer
Pi Pico RP2040

IIR filters in s1x14 fixed point with downsampler, second-order-sections, and design commands.

This project is an attempt to make an interactive infinite-impulse response (IIR) filter designer. The user specifies a filter type (lowpass, bandpass, highpass), the relevant frequencies and a parameter which sets the sharpness of cutoff. The actual filtering is done in an ISR running at 40 KHz, so we are designing audio-rate filters. The arithmetic is fixed point s1x14 format (16-bit total: a sign bit, one integer bit, 14 fractional bits). After designing a filter, the program may be set to do realtime filtering of audio signals using the ADC channel zero input. Actual filter input and output are displayed using a SPI DAC.

Fixed point s1x14 format filters are about 15 times faster to execute than floating point (about one second-order-section per microsecond) , but not accurate enough for multiple-pole, sharp cutoff, low bandwidth filters. There are two ways to fix this limitation. The first approach is downsample when you need low bandwidth, so that the ratio of bandwith to Nyquist rate is not too small. A CIC filter is used to drop the sample rate by 4, while avoiding aliasing. The second approach with s1x14 filters is to implement second-order-section filters for better numerical stability. This example implements 2,4, and 6 pole lowpass, 2 and 4 pole bandpass, and 2 pole highpass filters. There is optional 4x down-sampling.

All filters implement the "Direct Form II " implementation of a difference equation shown below (from the Matlab or Octave filter description).

    a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)   
                          - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

The y(n) term is the current filter output and x(n) is the current filter input. The a's and b's are constant coefficients, determined by some design process. Usually, a(1) is set to one. The filter design used here (getting a's and b's) is a slightly simplified version from Cookbook formulae for audio equalizer biquad filter coefficients by Robert Bristow-Johnson. You could also generate the second-order-section coefficients in Python, Octave, or Matlab.

All functions are controlled by a simple command sequence. In the commands below, the parameter Q is a measure of gain-change near the cutoff frequency, relative to the DC gain. In lowpass or high pass filters, it trades off passband flatness for faster cutoff. In bandpass filters it sets the sharpness of the filter. Some example Q values for a 2-pole lowpass are shown in the lp2 command below.

• lp2 cutoff_freq Q
  The lp2 command takes two parameters: the cutoff frequency in Hz, and a Q from 0.5-1.99.
  A Q=0.5 constructs a two-pole RC filter with slow cutoff.
  A Q=0.707 constructs a Butterworth filter, with optimally flat passband.
  A Q greater than ~0.8 becomes a Chebyshev filter with faster cutoff and some passband ripple. Range 1 to 1.3
  A Q of more than about 1.5 can saturate the fixed point numbers, unless used carefully.
  
• hp2 cutoff_freq Q
  Parameters have the same meaning as above.
• bp2 center_freq Q
  The frequency is the center frequency of the bandpass. Q sets the sharpness of the frequency range.
  Useful range of Q is about 1.0 to 20. 1/Q is approximately the fractional bandwidth. A Q of 10 means bandwidth is ~1/10 of the center frequency.
• lp4 cutoff_freq1 Q1 cutoff_freq2 Q2
  You need to specify two cutoffs and Qs, one for each second-order-section. The sections can be the same,
  or they can be optimized separately.
  This is a SOS 4-pole lowpass.
• lp6 cutoff_freq1 Q1 cutoff_freq2 Q2 cutoff_freq3 Q3
  You need to specify three cutoffs and Qs, one for each second-order-section.
  This is a SOS 6-pole lowpass.
• bp4 center_freq1 Q1 center_freq2 Q2
  A SOS 4-pole bandpass filter. Sharper than the 2-pole version, with steeper cutoffs.
• hp2lp2 hp_cutoff_freq1 Q1 lp_cutoff_freq2 Q2
  A 4-pole SOS bandpass filter which is flatter in the passband than bp2 or bp4
  Requires specifying a 2-pole highpass and lowpass filter for the band edges.
  --- display controls ---
• freq low high increment
  Control the plotting frequency range from low to high with ratio increment (>1, typically 1.1 to 2)
  For a very dense plot choose 1.05.
  
• color single_letter
  One of: g, r, b, c, m, w, y for the plot (green, red, blue, cyan, magenta, white, yellow),
  The chosen color will be used for the next Bode plot.
• downsample yes/no
  A yes performs a 4x sample slowdown before the filter. This shifts the filter frequencies lower.
• erase
  clears plot, draws axes
• dac yes/no
  A yes turns on a two channel SPI DAC displaying the filter input and output. The DAC is also turned on by the realtime command below.
• realtime yes/no
  A yes turns on channel zero of the ADC and applies it to the filter you just designed.
  Designing a new filter while realtime is turned on automatically turns it off!
  

Some Bode plot examples are shown below.
The first plot is for three bp4 filters with Q of 4 to 10(sharpest).
The second is the hp2lp2 4-pole bandpass filter with Q 0f 0.707 and cutoffs of (200,1000) and (300,700).
The third is 2, 4, and 6 pole lowpass filters all set to cutoff of 500 Hz and Qs chosen to be as flat as possible by educated guess of Qs.

Frequency grid lines are logarithmic, with decade, 2's (e.g. 20,200,200) and 5's (e.g. 50, 500, 5000) drawn. Amplitude in the upper plot is logarithmic, with decade lines drawn and labeled. Phase in the lower plot is a linear scale. At low frequency there may be some noise in the plot because the averaging time of the cross-correlation used to get the phase and magnitude did not average enough samples. To get phase and magnitude, the filter output was cross-correlated with the filter input sine wave, and a 90-degree phased shifted version of the sine wave. These are refered to as the I (in-phase) and Q (quatrature) channels of the correlation. You can think of the I and Q signals as the vector components of the output voltage. The average value of the magnitude of the vector is the output amplitide and the arctan of the ratio of the average of I and Q is the phase.

Realtime filtering example:
The filter is a 4-pole bandpass at 500 Hz with each SOS having a Q=10. It is quite narrow response.
First plot is the Bode plot. Note the peak at 500 Hz.
Second plot is the realtime response, the bottom trace is filter input, top trace is filter output.
The narrow response filters out all of the Fourier components of the square wave, except for the fundamental.
The actual commands:
filter cmd> color r
filter cmd> erase
filter cmd> freq 100 2000 1.1
filter cmd> bp4 500 10 500 10
filter cmd> realtime yes

The program is bulky because it handles the graphics to plot the frequency responses, supports a command language, computes the filter coeffients from user specifications, specifies several filter types, and runs a 40 KHz ISR to actually do the filtering. This ISR sweeps a frequency over the desired range using DDS, runs the desired digital filter using the DDS as input, then computes the cross-correlation of the filter output with the DDS input, and with a second DDS signal which is in quadrature with the filter input. From the average correlations, the graphics thread computes the Bode plots of filter phase and magnitude. The phase/magnitude calculations become inaccurate at very low frequency, or low filter output amplitude. You can see that in the bandpass plots above.

For everything to work, you need to attach a VGA monitor, a MCP4822 SPI DAC, and a signal source to ADC channel zero (source voltage must be between zero and Vdd volts). More hardware details are also in the code.
Photo of circuit.

Code, project ZIP

CIC downsampler.
For both Chebychev type 1 filters and Butterworth filters, the filters are stable down to bandwidth/Nyquist of 0.01 or so. With the full sample rate of 40KHz (Nyquist rate 20KHz), then downsampling by 4, the Nyquist rate is 5KHz. The 4th order CIC filter, plus Chebshev compensator has a rapid cutoff at 2 KHz, or 0.4 normalized to the new Nyquist rate. The first image below is the downsampler Bode plot. Note the fast cutoff at 2KHz, dropping to below 0.001 amplitude before the 5KHz Nyquest rate. The second plot is a Butterworth 2-pole lowpass feed by the downsampler. The cutoff of the lowpass is 0.025 of the downsampled Nyquist rate.

Copyright Cornell University May 27, 2023