Difference between revisions of "Digital 4-pole Butterworth Low-pass filter"
(→Calculating the DC Gain (k)) |
|||
(21 intermediate revisions by 3 users not shown) | |||
Line 1: | Line 1: | ||
== Overview == | == Overview == | ||
− | A digital 4-pole Butterworth low-pass filter is implemented as 2 cascaded biquads (2-pole topology) in the ''' | + | 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+ | + | <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 | + | 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 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> | + | 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> |
− | Note that in order to | + | 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: | There are 3 filter-related MCE parameters/commands: | ||
− | * fltr_type: to determine whether filter coefficients are hard-coded or configurable | + | * {{param|rc|fltr_type}}: to determine whether filter coefficients are hard-coded or configurable |
− | * fltr_coeff: to specify filter coefficients | + | * {{param|rc|fltr_coeff}}: to specify filter coefficients |
− | * fltr_rst: to reset the filter pipeline after changing coefficients | + | * {{param|rc|fltr_rst}}: to reset the filter pipeline after changing coefficients |
== Filter Specification == | == Filter Specification == | ||
− | Filter specification of the 4-pole Butterworth filter is determined by the filter coefficients | + | 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: | Assume: | ||
− | '''f<sub>cutoff</sub>''' is the frequency at which the signal is <math>\sqrt2</math> of its maximum value and | + | * '''''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 MHz / (33 × 100) = 15151.5 Hz. | ||
− | + | ''Note that the filter ''f''<sub>3dB</sub> scales if ''f''<sub>samp</sub> changes.'' | |
− | |||
− | ''Note that the filter f<sub>3dB</sub> scales if f<sub>samp</sub> changes.'' | ||
; Type 1 (hard-coded coefficients) | ; Type 1 (hard-coded coefficients) | ||
− | : | + | : ''supported in all rc '''except''' 5.0.7'' |
: DC amplification = 1217.9148 | : DC amplification = 1217.9148 | ||
− | : f<sub>cutoff</sub> / f<sub>samp</sub> = 122. | + | : ''f''<sub>cutoff</sub> / ''f''<sub>samp</sub> = 122.226 Hz / 15151 Hz |
− | : Gain @ | + | : Gain @ 200 Hz = 0.14189148 (wrt the DC gain) |
; Type 2 (hard-coded coefficients) | ; Type 2 (hard-coded coefficients) | ||
− | : | + | : ''supported '''only''' in rc version 5.0.7'' |
: DC amplification = 2044 | : DC amplification = 2044 | ||
− | : f<sub>cutoff</sub> / f<sub>samp</sub> = | + | : ''f''<sub>cutoff</sub> / ''f''<sub>samp</sub> = 75 Hz / 30000 Hz |
: A plot of the transfer function generated by simulating MCE firmware (VHDL) | : 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]] | [[File:mce_filter_type2_magnitude.png |150px]][[File:mce_filter_type2_phase.png|150px]] | ||
; Type 255 (configurable coefficients) | ; Type 255 (configurable coefficients) | ||
− | : ''supported in rc 5.1.0+'' | + | : ''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 | + | 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. | 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 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]] | [[Image: BW_filter2.png]] | ||
− | Here is the filter response for Scuba2 setting (f<sub>samp</sub>=9. | + | Here is the filter response for Scuba2 setting (''f''<sub>samp</sub> = 9.5 kHz, ''f''<sub>readout</sub> = 200 Hz): [[Media: Filter fs10kHz 200Hzdecimated response.ps ]] |
== Filter-related Software changes == | == Filter-related Software changes == | ||
Line 55: | Line 55: | ||
=== Determining the filter type === | === Determining the filter type === | ||
− | A | + | A {{param|rc|fltr_type}} parameter is introduced starting firmware 5.0.a: |
− | * fltr_type = 1 | + | * '''{{param|rc|fltr_type}} = 1''': ''f''<sub>cutoff</sub> / ''f''<sub>samp</sub> = 122.226 Hz / 15151 Hz |
− | * fltr_type = 2 | + | * '''{{param|rc|fltr_type}} = 2''': ''f''<sub>cutoff</sub> / ''f''<sub>samp</sub> = 75 Hz / 30000 Hz |
− | * fltr_type = 255 configurable coefficient | + | * '''{{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. | All earlier firmware revisions are set to filter type 1 by default, except version 5.0.7 which is set to filter type 2. | ||
Line 64: | Line 64: | ||
=== Software requirements === | === Software requirements === | ||
− | * MAS - use mas/trunk r493 or later to get access to the fltr_coeff MCE parameter. | + | * MAS - use mas/trunk r493 or later to get access to the {{param|rc|fltr_coeff}} MCE parameter. |
− | ** In [[mce. | + | ** In [[mce.cin]], specify <code>$fw_rev["rc"] = 0x5010000</code> or greater. |
* mce_script - use mce_script/trunk r764 or later | * mce_script - use mce_script/trunk r764 or later | ||
* In [[Mce config template system#experiment.cfg | experiment.cfg]]: | * In [[Mce config template system#experiment.cfg | experiment.cfg]]: | ||
Line 71: | Line 71: | ||
config_filter = 0; # do not write filter coefficients to the MCE | config_filter = 0; # do not write filter coefficients to the MCE | ||
config_filter = 1; # 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>] | + | 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 | + | filter_params = [ 32092, 15750, 31238, 14895, 0, 11]; # type 1 filter, ''f''<sub>cutoff</sub>/''f''<sub>samp</sub> = 122.226 Hz / 15.151 kHz |
− | filter_params = [ 32295, 15915, 32568, 16188, 3, 14]; # type 2 filter | + | filter_params = [ 32295, 15915, 32568, 16188, 3, 14]; # type 2 filter, ''f''<sub>cutoff</sub>/''f''<sub>samp</sub> = 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, ''f''<sub>cutoff</sub>/''f''<sub>samp</sub> = 75 Hz / 12.97 kHz | ||
== Code for simulations == | == Code for simulations == | ||
Line 86: | Line 88: | ||
* [[ Media: filter_pro2.txt ]] | * [[ Media: filter_pro2.txt ]] | ||
* [[ Media: low_pass_filter_model_pro.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 == | == Filter Coefficients == | ||
=== Butterworth Coefficients === | === Butterworth Coefficients === | ||
− | The filter coefficients are generated using ''fdatool'' ( | + | 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: | Once you launch the ''fdatool'', choose the following settings: | ||
Line 98: | Line 105: | ||
* Frequency Specifications: | * Frequency Specifications: | ||
− | ** For Type 1: | + | ** For Type 1: |
− | ** For Type 2: | + | *** ''f''<sub>samp</sub> = 12,195 Hz = 50 MHz/ (100 × 41) |
− | ** For Type 255: | + | *** ''f''<sub>cutoff</sub> = 100 Hz |
− | * | + | ** For Type 2: |
+ | *** ''f''<sub>samp</sub> = 30,000 Hz ≈ 50 MHz / (50 × 33) | ||
+ | *** ''f''<sub>cutoff</sub> = 75 Hz | ||
+ | ** For Type 255: | ||
+ | *** ''f''<sub>samp</sub> = 50 MHz / ({{param|sys|row_len}} × {{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 dB (half the passband power) | ||
Then click on Design Filter and you will get the following coefficients: | Then click on Design Filter and you will get the following coefficients: | ||
Line 108: | Line 121: | ||
::Numerator: 1 2 1 | ::Numerator: 1 2 1 | ||
::Denominator: 1 -1.9587428340882587 0.96134553442399129 | ::Denominator: 1 -1.9587428340882587 0.96134553442399129 | ||
− | ::Gain = 1/g<sub>1</sub> = 0.00065067508393319923 (not implemented) | + | ::Gain = 1/''g''<sub>1</sub> = 0.00065067508393319923 (not implemented) |
:Section 2: | :Section 2: | ||
::Numerator: 1 2 1 | ::Numerator: 1 2 1 | ||
::Denominator: 1 -1.9066292518523014 0.90916270571237567 | ::Denominator: 1 -1.9066292518523014 0.90916270571237567 | ||
− | ::Gain = 1/g<sub>2</sub>= 0.00063336346501859835 (not implemented) | + | ::Gain = 1/''g''<sub>2</sub>= 0.00063336346501859835 (not implemented) |
; Type 2 | ; Type 2 | ||
Line 121: | Line 134: | ||
::Numerator: 1 2 1 | ::Numerator: 1 2 1 | ||
::Denominator: 1 -1.9878047097960421 0.98804997058724808 | ::Denominator: 1 -1.9878047097960421 0.98804997058724808 | ||
− | :Gain = 1/(g<sub>1</sub> | + | :Gain = 1/(''g''<sub>1</sub> × ''g''<sub>2</sub>)= 0.0000000037280516432624239 (not implemented) |
; Type 255 | ; Type 255 | ||
− | : Plug in your desired f<sub> | + | : 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: | :Section 1: | ||
::Numerator: 1 2 1 | ::Numerator: 1 2 1 | ||
− | ::Denominator: 1 | + | ::Denominator: 1 <math>b^*_{11}</math> <math>b^*_{12}</math> |
− | :: Gain = 1/g<sub>1</sub> | + | :: Gain = 1/''g''<sub>1</sub> |
:Section 2: | :Section 2: | ||
::Numerator: 1 2 1 | ::Numerator: 1 2 1 | ||
− | ::Denominator: 1 | + | ::Denominator: 1 <math>b^*_{21}</math> <math>b^*_{22}</math> |
− | ::Gain = 1/g<sub>2</sub> | + | ::Gain = 1/''g''<sub>2</sub> |
=== Coefficient Quantization === | === 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. | 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 |
− | (10 | + | 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> < 16 and ''k''<sub>2</sub> < 32 | |
− | < | + | * Now, if 32 − ''k''<sub>1</sub> + ''k''<sub>2</sub> > 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).'' | |
− | |||
− | |||
− | 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) === | + | === 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. | + | 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>. | 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'''.) | + | :''(Note for firmware developers: see FILTER_GAIN_WIDTH and FILTER_SCALE_LSB in '''fsfb_calc_pack.vhd'''.) '' |
; Type 1 | ; 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 | ; 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. | 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: | + | 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 * | from numpy import * | ||
# Sum of k1 and k2 shifts: | # Sum of k1 and k2 shifts: | ||
Line 182: | Line 192: | ||
# Compute gain | # Compute gain | ||
gain = 16. / (2**k12 * (1 - f[0] + f[1]) * (1 - f[2] + f[3])) | gain = 16. / (2**k12 * (1 - f[0] + f[1]) * (1 - f[2] + f[3])) | ||
− | |||
print gain | print gain | ||
+ | # Result: 1217.8583043 | ||
== Useful Links == | == Useful Links == | ||
− | http://www.phas.ubc.ca/~mce/mcedocs/hardware/Firmware_block_spec/reaout_card/ | + | 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 | http://www.planetanalog.com/showArticle.jhtml?articleID=12802683 | ||
== Alternative to fdatool == | == Alternative to fdatool == | ||
− | If you do not have access to fdatool or do not have DSP toolbox installed, you can use the following | + | 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 | * butter: is a matlab function to generate butterworth coefficients | ||
* sos2tf: matlab function to break a transfer function to second-order sections | * sos2tf: matlab function to break a transfer function to second-order sections | ||
Line 198: | Line 208: | ||
num=[1 2 1] | num=[1 2 1] | ||
denum=[1 -1.9711486088510415 0.97139181456687917] | denum=[1 -1.9711486088510415 0.97139181456687917] | ||
− | max(dbode(num, denum, 1/ | + | 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?! | '''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]] |
Latest revision as of 13:01, 17 May 2017
Contents
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, k1 and k2, 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 bxy 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:
- fcutoff 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
- fsamp 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 f3dB scales if fsamp changes.
- Type 1 (hard-coded coefficients)
- supported in all rc except 5.0.7
- DC amplification = 1217.9148
- fcutoff / fsamp = 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
- fcutoff / fsamp = 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 b11, b12, b21, b22 and truncation factors, k1, k2 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 (fsamp = 9.5 kHz, freadout = 200 Hz): Media: Filter fs10kHz 200Hzdecimated response.ps
Determining the filter type
A fltr_type parameter is introduced starting firmware 5.0.a:
- fltr_type = 1: fcutoff / fsamp = 122.226 Hz / 15151 Hz
- fltr_type = 2: fcutoff / fsamp = 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.
- In mce.cin, specify
- 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 [b11, b12, b21, b22, k1, k2], 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:
Digital, time domain simulation codes
Python-based model of a single biquad, in the time domain, for exploring digitization and dynamic range:
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:
- The attenuation at fcutoff 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/g1 = 0.00065067508393319923 (not implemented)
- Section 2:
- Numerator: 1 2 1
- Denominator: 1 -1.9066292518523014 0.90916270571237567
- Gain = 1/g2= 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/(g1 × g2)= 0.0000000037280516432624239 (not implemented)
- Type 255
- Plug in your desired fsamp and fcutoff 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/g1
- Section 2:
- Numerator: 1 2 1
- Denominator: 1 <math>b^*_{21}</math> <math>b^*_{22}</math>
- Gain = 1/g2
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 bxy via: <math>b_{xy} = \left\lfloor\left|b^*_{xy}\right| \times 2^{14}\right\rfloor</math> (i.e. take the absolute value, multiply by 214 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>
- 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 g1 and g2, an inter-biquad truncation of k2 (implemented as a binary shift) and k1, 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 (b11 = 32092, b12 = 15750, b21 = 31238, b22 = 14895, k1 = 0, k2 = 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.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?!