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 Read-out 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}} . 2^{-k_2} . \frac {1+2z^{-1}+z^{-2}}{1+b_{21} z^{-1}+b_{22} z^{-2}} . {2^{-k_1}} </math>

Filter coefficients b₁₁, b₁₂, b₂₁ b₂₂, and truncation factors, k₁ and k₂, are hard coded in older firmware and just recently (version 5.1.0+) they can be programmed through software.

 wb rca fltr_coeff b₁₁ b₁₂ b₂₁ b₂₂ k₁ k₂

Note that in order to accomodate the quantization effects, you have to follow the recipe prescribed here 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 b₁₁, b₁₂, b₂₁, b₂₂, 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 50MHz / (33 * 100) = 15151.5Hz.

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.226Hz / 15151Hz; Gain @ 200Hz=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 = 75Hz / 30000Hz; A plot of the transfer function generated by simulating MCE firmware (VHDL)

Type 255 (configurable coefficients): supported in rc 5.1.0+ Filter coefficients b₁₁, b₁₂, b₂₁, b₂₂ and truncation factors, k₁, k₂ are parametrized and programmable by software.

The fitler type can be determined by reading back the MCE parameter called fltr_type (rb rc1 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.5kHz, f_readout=200Hz): 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 fcutoff / fsamp = 122.226Hz / 15151Hz
fltr_type = 2 fcutoff / fsamp = 75Hz / 30000Hz
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.cfg, specify $fw_rev["rc"] = $5010000 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, fc/fs= 122.226Hz / 15.151kHz
 filter_params = [ 32295, 15915, 32568, 16188, 3, 14];  # type 2 filter, fc/fs= 75Hz / 30kHz
 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, fc/fs = 75Hz / 12.97kHz

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:

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 = 12195 = (50000/100*41), F_cutoff=100
- For Type 2: F_samp = 30000 ~= (50000/50*33), F_cutoff=75
- For Type 255: F_ssamp = You specify (50000/row_len*num_rows), F_cutoff=You specify (<F_readout / 2)
- (The attenuation at F_cutoff is fixed at 3dB (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_s and f_c in Matlab fdatool and you will get a set of coefficients:

Section 1:

Numerator: 1 2 1

Denominator: 1 B11 B12

Gain = 1/g₁

Section 2:

Numerator: 1 2 1

Denominator: 1 B21 B22

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.

Multiply each of the Bxx coefficients by 2¹⁴ and then if negative multiply by -1.
k2 = 1 + floor ( log₂ (g₁) ) the inter-stage truncation
k1 = floor ( log₂ (g₂) ) - 10 the output-truncation

(10 is a historical constant !!!)

For example, if we use this formula to quantize Type 1 coefficients listed above:
b11 = 32092
b12 = 15750
b21 = 31238
b22 = 14895
k1 = 0
k2 = 11

k1<16 and k2<32
Now if 32-k1+k2 >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: k2=11, k1=0, the filter gain is estimated at 1184, but vhdl simulation results in a gain of 1216.
Type 2: k2=14, k1=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:

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]))
# Answer: 1217.8583043
print gain

Useful Links

http://www.phas.ubc.ca/~mce/mcedocs/hardware/Firmware_block_spec/reaout_card/fsfb_calculations.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/fs) 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?!

Digital 4-pole Butterworth Low-pass filter

Contents

Overview

Filter Specification

Filter-related Software changes

Determining the filter type

Software requirements

Code for simulations

Digital, time domain simulation codes

Filter Coefficients

Butterworth Coefficients

Coefficient Quantization

Calculating the DC Gain (k)

Useful Links

Alternative to fdatool

Navigation menu

Views

Personal tools

Navigation

MCE

Search

Tools