View source for 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''<sub>1</sub> and ''k''<sub>2</sub>, are hard coded in older firmware (before '''version 5.1.0'''), but with modern firmware, they can be programmed through software:

  wb rca {{param|rc|fltr_coeff}} ''b''<sub>11</sub> ''b''<sub>12</sub> ''b''<sub>21</sub> ''b''<sub>22</sub> ''k''<sub>1</sub> ''k''<sub>2</sub>

where ''b<sub>xy</sub>'' 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|''Filter Coefficients'']] section below when specifying the coefficients and truncation factors.

There are 3 filter-related MCE parameters/commands:
* {{param|rc|fltr_type}}: to determine whether filter coefficients are hard-coded or configurable
* {{param|rc|fltr_coeff}}: to specify filter coefficients
* {{param|rc|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''<sub>cutoff</sub>''' 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''<sub>samp</sub>''' is the sampling frequency, calculated as <math>\frac{50\;\mathrm{MHz}}{\mathrm{num\_rows} \times \mathrm{row\_len}}</math>. For example 50&nbsp;MHz / (33 &times; 100) = 15151.5&nbsp;Hz.

''Note that the filter ''f''<sub>3dB</sub> scales if ''f''<sub>samp</sub> changes.''

; Type 1 (hard-coded coefficients)
: ''supported in all rc '''except''' 5.0.7''
: DC amplification = 1217.9148 
: ''f''<sub>cutoff</sub> / ''f''<sub>samp</sub> = 122.226&nbsp;Hz / 15151&nbsp;Hz
: Gain @ 200&nbsp;Hz = 0.14189148 (wrt the DC gain) 

; Type 2 (hard-coded coefficients)
: ''supported '''only''' in rc version 5.0.7''
: DC amplification = 2044
: ''f''<sub>cutoff</sub> / ''f''<sub>samp</sub> = 75&nbsp;Hz / 30000&nbsp;Hz
: A plot of the transfer function generated by simulating MCE firmware (VHDL) 
[[File:mce_filter_type2_magnitude.png |150px]][[File:mce_filter_type2_phase.png|150px]]

; Type 255 (configurable coefficients)
: ''supported in rc 5.1.0+''
: Quantized filter coefficients ''b''<sub>11</sub>, ''b''<sub>12</sub>, ''b''<sub>21</sub>, ''b''<sub>22</sub> and truncation factors, ''k''<sub>1</sub>, ''k''<sub>2</sub> are parameterized and programmable by software. 

The fitler type can be determined by reading back the MCE parameter called {{param|rc|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): 
[[Image: BW_filter2.png]]

Here is the filter response for Scuba2 setting (''f''<sub>samp</sub> = 9.5&nbsp;kHz, ''f''<sub>readout</sub> = 200&nbsp;Hz): [[Media: Filter fs10kHz 200Hzdecimated response.ps ]]

== Filter-related Software changes ==

=== Determining the filter type ===

A {{param|rc|fltr_type}} parameter is introduced starting firmware 5.0.a:
* '''{{param|rc|fltr_type}} = 1''': ''f''<sub>cutoff</sub> / ''f''<sub>samp</sub> = 122.226&nbsp;Hz / 15151&nbsp;Hz
* '''{{param|rc|fltr_type}} = 2''': ''f''<sub>cutoff</sub> / ''f''<sub>samp</sub> = 75&nbsp;Hz / 30000&nbsp;Hz
* '''{{param|rc|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 {{param|rc|fltr_coeff}} MCE parameter.
** In [[mce.cin]], specify <code>$fw_rev["rc"] = 0x5010000</code> or greater.
* mce_script - use mce_script/trunk r764 or later
* In [[Mce config template system#experiment.cfg | 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''<sub>11</sub>, ''b''<sub>12</sub>, ''b''<sub>21</sub>, ''b''<sub>22</sub>, ''k''<sub>1</sub>, ''k''<sub>2</sub>],  e.g.:
  filter_params = [ 32092, 15750, 31238, 14895, 0, 11];  # type 1 filter, ''f''<sub>cutoff</sub>/''f''<sub>samp</sub> = 122.226&nbsp;Hz / 15.151&nbsp;kHz
  filter_params = [ 32295, 15915, 32568, 16188, 3, 14];  # type 2 filter, ''f''<sub>cutoff</sub>/''f''<sub>samp</sub> = 75&nbsp;Hz / 30&nbsp;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, ''f''<sub>cutoff</sub>/''f''<sub>samp</sub> = 75&nbsp;Hz / 12.97&nbsp;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 [[#Alternative to fdatool|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''<sub>samp</sub> = 12,195&nbsp;Hz = 50&nbsp;MHz/ (100 &times; 41)
*** ''f''<sub>cutoff</sub> = 100&nbsp;Hz
** For Type 2:
*** ''f''<sub>samp</sub> = 30,000&nbsp;Hz ≈ 50&nbsp;MHz / (50 &times; 33)
*** ''f''<sub>cutoff</sub> = 75&nbsp;Hz
** For Type 255:
*** ''f''<sub>samp</sub> = 50&nbsp;MHz / ({{param|sys|row_len}} &times; {{param|sys|num_rows}})
*** ''f''<sub>cutoff</sub>=''f''<sub>readout</sub> / 2 (= ''f''<sub>samp</sub> / {{param|cc|data_rate}} / 2, if using internal CC triggering)
* The attenuation at ''f''<sub>cutoff</sub> is fixed at 3&nbsp;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''<sub>1</sub> = 0.00065067508393319923 (not implemented)
:Section 2:
::Numerator: 1  2  1  
::Denominator: 1  -1.9066292518523014  0.90916270571237567
::Gain = 1/''g''<sub>2</sub>= 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''<sub>1</sub> &times; ''g''<sub>2</sub>)= 0.0000000037280516432624239  (not implemented)

; Type 255
: Plug in your desired ''f''<sub>samp</sub> and ''f''<sub>cutoff</sub> 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''<sub>1</sub>
:Section 2:
::Numerator: 1  2  1  
::Denominator: 1 <math>b^*_{21}</math> <math>b^*_{22}</math>
::Gain = 1/''g''<sub>2</sub>

=== 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<sub>xy</sub>'' 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<sup>14</sup> 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''<sub>1</sub> &lt; 16 and ''k''<sub>2</sub> &lt; 32
* Now, if 32 &minus; ''k''<sub>1</sub> + ''k''<sub>2</sub> &gt; 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''<sub>1</sub> and ''g''<sub>2</sub>, an inter-biquad truncation of ''k''<sub>2</sub> (implemented as a binary shift) and ''k''<sub>1</sub>, 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''<sub>2</sub> = 11, ''k''<sub>1</sub> = 0, the filter gain is estimated at 1184, but VHDL simulation results in a gain of 1216. 
; Type 2
:''k''<sub>2</sub> = 14, ''k''<sub>1</sub> = 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''<sub>11</sub> = 32092, ''b''<sub>12</sub> = 15750, ''b''<sub>21</sub> = 31238, ''b''<sub>22</sub> = 14895, ''k''<sub>1</sub> = 0, ''k''<sub>2</sub> = 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?!

[[Category:Readout Card Firmware]]