Digital 4-pole Butterworth Low-pass filter

Overview

A digital 4-pole Butterworth low-pass filter is implemented as 2 cascaded biquads (2-pole topology) in the Readout card firmware of the MCE. The transfer function of the IIR filter is: <math> H(z) = \frac {1+2 z^{-1}+z^{-2}}{1+b^*_{11} z^{-1}+b^*_{12} z^{-2}} \cdot 2^{-k_2} \cdot \frac {1+2z^{-1}+z^{-2}}{1+b^*_{21} z^{-1}+b^*_{22} z^{-2}} \cdot {2^{-k_1}} </math>

Filter coefficients <math>b^*_{11}</math>, <math>b^*_{12}</math>, <math>b^*_{21}</math>, <math>b^*_{22}</math>, and truncation factors, k₁ and k₂, are hard coded in older firmware (before version 5.1.0), but with modern firmware, they can be programmed through software:

 wb rca fltr_coeff b11 b12 b21 b22 k1 k2

where b_xy is a quantized version of <math>b^*_{xy}</math> computed as <math>b_{xy} = \left\lfloor\left|b^*_{xy}\right| \times 2^{14}\right\rfloor</math>, where <math>\lfloor\bullet\rfloor</math> indicates the floor function. Note that in order to accommodate the quantization effects, you have to follow the recipe prescribed in the Filter Coefficients section below when specifying the coefficients and truncation factors.

There are 3 filter-related MCE parameters/commands:

• fltr_type: to determine whether filter coefficients are hard-coded or configurable
• fltr_coeff: to specify filter coefficients
• fltr_rst: to reset the filter pipeline after changing coefficients

Filter Specification

Filter specification of the 4-pole Butterworth filter is determined by the filter coefficients <math>b^*_{11}</math>, <math>b^*_{12}</math>, <math>b^*_{21}</math>, <math>b^*_{22}</math>, described earlier.

Assume:

• f_cutoff is the frequency at which the signal after filtering is <math>\sqrt2</math> of its maximum value (nominally the read-out rate Nyquist frequency), and
• f_samp is the sampling frequency, calculated as <math>\frac{50\;\mathrm{MHz}}{\mathrm{num\_rows} \times \mathrm{row\_len}}</math>. For example 50 MHz / (33 × 100) = 15151.5 Hz.

Note that the filter f_3dB scales if f_samp changes.

Type 1 (hard-coded coefficients)

supported in all rc except 5.0.7

DC amplification = 1217.9148

f_cutoff / f_samp = 122.226 Hz / 15151 Hz

Gain @ 200 Hz = 0.14189148 (wrt the DC gain)

Type 2 (hard-coded coefficients)

supported only in rc version 5.0.7

DC amplification = 2044

f_cutoff / f_samp = 75 Hz / 30000 Hz

A plot of the transfer function generated by simulating MCE firmware (VHDL)

Type 255 (configurable coefficients)

supported in rc 5.1.0+

Quantized filter coefficients b₁₁, b₁₂, b₂₁, b₂₂ and truncation factors, k₁, k₂ are parameterized and programmable by software.

The fitler type can be determined by reading back the MCE parameter called fltr_type. (supported in firmware 5.0.a+)

The following plot shows the magnitude and the phase of Type 1 filter.

Here is Elia's original plot Media:BW_filter.ps. Then we see the same plot with the impulse-response also added (by Joe, Nov 29, 2007):

Here is the filter response for Scuba2 setting (f_samp = 9.5 kHz, f_readout = 200 Hz): Media: Filter fs10kHz 200Hzdecimated response.ps

Filter-related Software changes

Determining the filter type

A fltr_type parameter is introduced starting firmware 5.0.a:

• fltr_type = 1: f_cutoff / f_samp = 122.226 Hz / 15151 Hz
• fltr_type = 2: f_cutoff / f_samp = 75 Hz / 30000 Hz
• fltr_type = 255: configurable coefficient

All earlier firmware revisions are set to filter type 1 by default, except version 5.0.7 which is set to filter type 2.

Software requirements

• MAS - use mas/trunk r493 or later to get access to the fltr_coeff MCE parameter.
  • In mce.cin, specify $fw_rev["rc"] = 0x5010000 or greater.
• mce_script - use mce_script/trunk r764 or later
• In experiment.cfg:
  • specify one of:

 config_filter = 0;  #  do not write filter coefficients to the MCE
 config_filter = 1;  #  write filter coefficients to the MCE

and when config_filter = 1, the filter coefficients to be written to the MCE must be specified in the "filter_params" parameter, in the order [b₁₁, b₁₂, b₂₁, b₂₂, k₁, k₂], e.g.:

 filter_params = [ 32092, 15750, 31238, 14895, 0, 11];  # type 1 filter, fcutoff/fsamp = 122.226 Hz / 15.151 kHz
 filter_params = [ 32295, 15915, 32568, 16188, 3, 14];  # type 2 filter, fcutoff/fsamp = 75 Hz / 30 kHz
 filter_params = [ 32295, 15915, 32568, 16188, 5, 12];  # type 2 filter with more inter-biquad dynamic range (recommended)
 filter_params = [ 32297, 15934, 31683, 15320, 0, 11];  # SCUBA-2 filter after 2011-06, fcutoff/fsamp = 75 Hz / 12.97 kHz

Code for simulations

Note that as of mce_script r902, mce_data.py contains code for handling the MCE Butterworth filters. This can be used to load filter parameters from a runfile, and to get the complex frequency response of a filter. See mce_data.py#MCEButterworth_class.

A smaller, quicker python example: Media: mce_filt_py.txt

In this little IDL program you can find the functional form as a function of the frequency of the filter. First, Elia's original program, then as modified by Joe on November 29, 2007, then modified by Mandana on Jan. 14, 2009 for Scuba2 numbers:

• Media: filter_pro.txt
• Media: filter_pro2.txt
• Media: low_pass_filter_model_pro.txt

Digital, time domain simulation codes

Python-based model of a single biquad, in the time domain, for exploring digitization and dynamic range:

• Media: biquad_time_domain_py.txt

Filter Coefficients

Butterworth Coefficients

The filter coefficients are generated using fdatool (Filter-design & Analysis tool, part of DSP Toolbox) in MATLAB/Simulink. These coefficients are floating numbers and in order to feed them to MCE firmware, they need to be quantized by converting to signed binary fractional (SBF) 1.14 format. You may use alternate methods to generate coefficients. Here, we explain the method using fdatool.

Once you launch the fdatool, choose the following settings:

• Response Type: Low Pass
• Design Method: Butterworth
• Filter Order: 4

• Frequency Specifications:
  • For Type 1:
    • f_samp = 12,195 Hz = 50 MHz/ (100 × 41)
    • f_cutoff = 100 Hz
  • For Type 2:
    • f_samp = 30,000 Hz ≈ 50 MHz / (50 × 33)
    • f_cutoff = 75 Hz
  • For Type 255:
    • f_samp = 50 MHz / (row_len × num_rows)
    • f_cutoff=f_readout / 2 (= f_samp / data_rate / 2, if using internal CC triggering)
• The attenuation at f_cutoff is fixed at 3 dB (half the passband power)

Then click on Design Filter and you will get the following coefficients:

Type 1

Section 1:

Numerator: 1 2 1

Denominator: 1 -1.9587428340882587 0.96134553442399129

Gain = 1/g₁ = 0.00065067508393319923 (not implemented)

Section 2:

Numerator: 1 2 1

Denominator: 1 -1.9066292518523014 0.90916270571237567

Gain = 1/g₂= 0.00063336346501859835 (not implemented)

Type 2

Section 1:

Numerator: 1 2 1

Denominator: 1 -1.9711486088510415 0.97139181456687917

Section 2:

Numerator: 1 2 1

Denominator: 1 -1.9878047097960421 0.98804997058724808

Gain = 1/(g₁ × g₂)= 0.0000000037280516432624239 (not implemented)

Type 255

Plug in your desired f_samp and f_cutoff in MATLAB fdatool and you will get a set of coefficients:

Section 1:

Numerator: 1 2 1

Denominator: 1 <math>b^*_{11}</math> <math>b^*_{12}</math>

Gain = 1/g₁

Section 2:

Numerator: 1 2 1

Denominator: 1 <math>b^*_{21}</math> <math>b^*_{22}</math>

Gain = 1/g₂

Coefficient Quantization

The quantization format is signed-binary fractional (SBF) 1.14 and then 2's complement to be able to account for positive and negative values.

• Compute the quantised coefficients b_xy via: <math>b_{xy} = \left\lfloor\left|b^*_{xy}\right| \times 2^{14}\right\rfloor</math> (i.e. take the absolute value, multiply by 2¹⁴ and then drop the fractional part).
• <math>k_1 = \left\lfloor \log_2 g_2 \right\rfloor - 10</math> the output-truncation (10 is a historical constant)
• <math>k_2 = 1 + \left\lfloor \log_2 g_1 \right\rfloor</math> the inter-stage truncation

For example, if we use this formula to quantize Type 1 coefficients listed above:

<math>b_{11} = \lfloor 1.9587428340882587 \times 2^{14}\rfloor = 32092</math>

<math>b_{12} = \lfloor 0.96134553442399129 \times 2^{14}\rfloor = 15750 </math>

<math>b_{21} = \lfloor 1.9066292518523014 \times 2^{14}\rfloor = 31238 </math>

<math>b_{22} = \lfloor 0.90916270571237567 \times 2^{14}\rfloor = 14895 </math>

<math>k_1 = \left\lfloor \log_2 \left(1 / 0.00065067508393319923 \right)\right\rfloor - 10 = 0</math>

<math>k_2 = 1 + \left\lfloor \log_2 \left(1 / 0.00063336346501859835 \right)\right\rfloor = 11</math>

• k₁ < 16 and k₂ < 32
• Now, if 32 − k₁ + k₂ > 32, then you have to talk to UBC!

The following ONLY concerns the MCE firmware developers:

Prior to firmware 5.1.0, the coefficients were hardcoded in fsfb_calc_pack.vhd (FILTER_B11_COEF, FILTER_B12_COEF, FILTER_B21_COEFF, FILTER_B22_COEFF, FILTER_GAIN_WIDTH, FILTER_SCALE_LSB).

Calculating the DC Gain (k)

Note that in firmware, instead of g₁ and g₂, an inter-biquad truncation of k₂ (implemented as a binary shift) and k₁, number of bits dropped after the second biquad, are implemented.

Hence, the overall gain becomes <math>\frac {g_1 g_2} {2^{(k_1 + k_2)}}</math>.

(Note for firmware developers: see FILTER_GAIN_WIDTH and FILTER_SCALE_LSB in fsfb_calc_pack.vhd.)

Type 1

k₂ = 11, k₁ = 0, the filter gain is estimated at 1184, but VHDL simulation results in a gain of 1216.

Type 2

k₂ = 14, k₁ = 3, the filter gain is estimated at 2046, but VHDL simulation results in a gain of 2181.

The gain difference can be attributed to the rounding effects associated with fixed-width arithmetic.

To roughly compute the gain, including the effects from truncating the coefficients, reverse the quantization process to obtain the floating point values associated with your coefficients, and plug into the filter definition. Here is a rough python program which computes the true gain for our Type 1 filter (b₁₁ = 32092, b₁₂ = 15750, b₂₁ = 31238, b₂₂ = 14895, k₁ = 0, k₂ = 11):

from numpy import *
# Sum of k1 and k2 shifts:
k12 = 11+0
# Quantized coefficients b11, b12, b21, b22:
n = array([32092, 15750, 31238, 14895]).astype('float')
# Re-float
f = n/2**14
# Compute gain
gain = 16. / (2**k12 * (1 - f[0] + f[1]) * (1 - f[2] + f[3]))
print gain
# Result: 1217.8583043

Useful Links

http://www.phas.ubc.ca/~mce/mcedocs/hardware/Firmware_block_spec/reaout_card/fsfb_calculations_rev1.11.pdf

http://www.planetanalog.com/showArticle.jhtml?articleID=12802683

Alternative to fdatool

If you do not have access to fdatool or do not have DSP toolbox installed, you can use the following MATLAB functions (or equivalent in other packages) to get the coefficients:

• butter: is a matlab function to generate butterworth coefficients
• sos2tf: matlab function to break a transfer function to second-order sections
• bilinear: converts an s-domain (continuous) transfer function to z-domain (discrete) transfer function.
• Then run the following in matlab (or whatever syntax your tool has):

  num=[1 2 1]
  denum=[1 -1.9711486088510415  0.97139181456687917]
  max(dbode(num, denum, 1/fsamp) but also look at the bode-plot to make sure it is flat, otherwise read the max from the flat portion of the plot instead of the max function.

Note: When I used coefficients generated by [a,b] = butter(2, 50/20000); [SOS,G]=tf2sos(a,b), the filter was not as robust. not sure why?!