Fixed-Point Filter Design: Where comes the DC-Gain difference between the bode plot and the step response from?

Question

Introduction: What I wanted to achieve
I am currently trying to implement a digital fixed-point filter on a STM32-Microcontroller by using the ARM CMSIS-DSP Library Link to reference manual.

The function I'm most interested in is called arm_biquad_cascade_df1_fast_q15(), which implements a Direct Form I biquad filter as follows:

The corresponding transfer function for such a filter stage is: $$ H(z) = \frac{2^{postShift} \cdot b_0 \cdot z^2 + 2^{postShift} \cdot b_1 \cdot z + 2^{postShift} \cdot b_2}{z^2 - 2^{postShift} \cdot a_1 \cdot z - 2^{postShift} \cdot a_2} $$ EDIT: I originally entered plus signs in the denominator by mistake. This does not change anything in my calculations but might have confused you.

With the $a_x,b_x$ coefficients being in 1.15 signed format with a range of $[-1 +1)$.
A scaling parameter $postShift$ is implemented inside this structure to tackle filter designs whose coefficients would exceed this $[-1 +1)$ variable range.

The target filter to implement is a $4^{th}$-order Butterworth lowpass with a cut-off frequency of $f_c$ = 30kHz at a sampling frequency of $f_s$ = 1.2MHz.

My strange observation / my question

As far as I can tell, my fixed-point design was sucessful.
The bode plot of my fixed-point filter looks almost identical to the "original" double-precision, which is exactly what I wanted.

Also the simulated step-response of my fixed-point filter looks as I would expect it.

Unfortunately, there is something that I can't really explain. The bode plot of my fixed-point design says that my DC-Gain is ~1.0012. The settling value of the fixed-point filter step response however, says that my DC-Gain is only ~0.9738.

I know that this difference is rather small but I'm still wondering where it comes from and how I could possibly prevent/reduce it. Does anybody know what I'm observing here or what I'm doing wrong??

My design process in matlab for reference
The design process in matlab starts with the calculation of the poles, zeros and gain of my target filter in double-precision using the butter() function:

fs = 1200e3;    %Sampling frequency [Hz]
fc = 30e3;      %Butterworth filter cuttoff frequency [Hz]
n = 4;          %Butterworth filter order

%Get low-pass butterworth filter zeroes and poles and gain
[z, p, g] = butter(n, 2 * fc / fs, 'low');

The zeros are:
[-1, -1, -1, -1]
The poles are:
[0.862966697518470 + 0.052305361884381i,
0.862966697518470 - 0.052305361884381i,
0.931900246037618 + 0.136363232169535i
0.931900246037618 - 0.136363232169535i]
And the gain is:
3.123897691708262e-05

Further, these values are transformed into the coefficients of my second order section filter (also in double precision) by using tf2sos().
(Note: The additional parameters 'down' and 'two' should order the poles in my second order sections such that the peak-roundoff noise is minimized)

sos_matrix = zp2sos(z, p, g, 'down', 'two');

nCascadedFilterSections = length(sos_matrix(:,1));

b0_double =  sos_matrix(:, 1);
b1_double =  sos_matrix(:, 2);
b2_double =  sos_matrix(:, 3);
a1_double = -sos_matrix(:, 5);
a2_double = -sos_matrix(:, 6);

To implement this filter, I will need to cascade two second order stages. The double-precision coefficients for the first section are:
b0_double(1) = 0.000550374291902
b1_double(1) = 0.001100748583804
b2_double(1) = 0.000550374291902
a1_double(1) = 1.863800492075236
a2_double(1) = -0.887032999652695
And for the second section:
b0_double(2) = 0.056759513256195
b1_double(2) = 0.113519026512390
b2_double(2) = 0.056759513256195
a1_double(2) = 1.7259333950369416
a2_double(2) = -0.747447371907791

From here on, one can already see, that the coefficients of $a_1$ are out of my 1.15 fixed-point range. Hence, I will have to set the postShift value to 1, which allows me to scale all the fixed point coefficients by 0.5 such that the $[-1 +1)$ variable range is not exceeded:
postShift = 1;

The fixed-point coefficients can then be calculated using matlabs fixed-point variable type fi(). Herefore, I scale each coefficient by 2 to the power of -postShift in double precision and then perform the conversion to 1.15 fixed-point.

b0_fixed = fi(b0_double * 2^-postShift, true, 16, 15, fimath('RoundingMethod', 'Nearest'));  
b1_fixed = fi(b1_double * 2^-postShift, true, 16, 15, fimath('RoundingMethod', 'Nearest'));  
b2_fixed = fi(b2_double * 2^-postShift, true, 16, 15, fimath('RoundingMethod', 'Nearest'));  
a1_fixed = fi(a1_double * 2^-postShift, true, 16, 15, fimath('RoundingMethod', 'Nearest'));  
a2_fixed = fi(a2_double * 2^-postShift, true, 16, 15, fimath('RoundingMethod', 'Nearest'));  

The resulting coefficients for the first stage are:
b0_fixed(1) = 0.000274658203125
b1_fixed(1) = 0.000549316406250
b2_fixed(1) = 0.000274658203125
a1_fixed(1) = 0.931915283203125
a2_fixed(1) = -0.443511962890625
And for the second stage:
b0_fixed(2) = 0.028381347656250
b1_fixed(2) = 0.056762695312500
b2_fixed(2) = 0.028381347656250
a1_fixed(2) = 0.862976074218750
a2_fixed(2) = -0.373718261718750

The fixed-point filter is now designed. By calculating it's magnitude and phase, one is able to determine whether the 1.15 coefficient format can approximate the double-precision filter successful or not.

This calculation is done by evaluating the transfer function $H(z)$ for each filter stage at different frequencies $f$ with $z = e^{j \cdot 2 \cdot \pi \cdot \frac{f}{f_s}}$.
The magnitude of the complete filter is calculated by multiplying the magnitude of each stage. The total phase however, is calculated by adding each stage's phase.
(Note: The fixed-point coefficients are casted to a double-precision variable type and are also scaled by 2 to the power of postShift before passing them to the anonymous function)

%Anonymous Filter Transfer-Function
H = @(z, b0, b1, b2, a1, a2) ((b0 .* z.^2 + b1 .* z + b2) ./ (z.^2 - a1 .* z - a2));

f = 0 : 1 : fs / 2 - 1;         %Frequency vector [Hz]
z = exp(1i * 2 * pi * f ./ fs); %corresponding z vector
mag_fixed = 1;
ph_fixed = 0;

for (section = 1 : nCascadedFilterSections)
    mag_fixed = mag_fixed .* abs(H(z, ...
        double(b0_fixed(section)) * 2^postShift, ...
        double(b1_fixed(section)) * 2^postShift, ...
        double(b2_fixed(section)) * 2^postShift, ...
        double(a1_fixed(section)) * 2^postShift, ...
        double(a2_fixed(section)) * 2^postShift));
    
    ph_fixed = ph_fixed + rad2deg(unwrap(angle(H(z, ...
        double(b0_fixed(section)) * 2^postShift, ...
        double(b1_fixed(section)) * 2^postShift, ...
        double(b2_fixed(section)) * 2^postShift, ...
        double(a1_fixed(section)) * 2^postShift, ...
        double(a2_fixed(section)) * 2^postShift))));
end

The computed magnitudes and phases are visible in the bode plot.

Finally, the step response is calculated by simulating the fixed-point arithmetics exactly as it would be done on my microcontroller.

To do this, I set up the fixed point properties (precision and overflow-behaviour of the multiplier / adder) as they are documented in the CMSIS-DSP reference (link on top) using fimath().

fixedPointFilterProperties = fimath('CastBeforeSum', true, ...
    'OverflowAction', 'Wrap', ...
    'RoundingMethod', 'Floor', ...
    'ProductWordLength', 32, ...
    'ProductFractionLength', 30, ...
    'ProductMode', 'SpecifyPrecision', ...
    'SumWordLength', 32, ...
    'SumFractionLength', 30, ...
    'SumMode', 'SpecifyPrecision');

These math properties must be passed to the fixed point coefficient variables as follows:

b0_fixed = setfimath(b0_fixed, fixedPointFilterProperties);
b1_fixed = setfimath(b1_fixed, fixedPointFilterProperties);
b2_fixed = setfimath(b2_fixed, fixedPointFilterProperties);
a1_fixed = setfimath(a1_fixed, fixedPointFilterProperties);
a2_fixed = setfimath(a2_fixed, fixedPointFilterProperties);

My step signal is then created as a vector:

t = 0: 1 / fs : 1e-3; %Discrete Time vector [s]

stepHight = 0.5; %Height of the step input [-]
step = [zeros(1,20), stepHight * ones(1,length(t) - 20)]; %Step signal

Finally, the step response can be calculated according to the filter's difference equation:

in = fi(step, true, 16, 15, fixedPointFilterProperties);

%Define fixed-point filter output vector
out = fi(zeros(1, length(in)), true, 16, 15, fixedPointFilterProperties);

%Filter as it would be done by the CMSIS Library function
for (section = 1 : nCascadedFilterSections)
    for (i = 1 : length(in))
        if (i == 1)
            b0_prod = in(i)    * b0_fixed(section);
            b1_prod = 0        * b1_fixed(section);
            b2_prod = 0        * b2_fixed(section);
            a1_prod = 0        * a1_fixed(section);
            a2_prod = 0        * a2_fixed(section);
        elseif (i == 2)
            b0_prod = in(i)    * b0_fixed(section);
            b1_prod = in(i-1)  * b1_fixed(section);
            b2_prod = 0        * b2_fixed(section);
            a1_prod = out(i-1) * a1_fixed(section);
            a2_prod = 0        * a2_fixed(section);
        else
            b0_prod = in(i)    * b0_fixed(section);
            b1_prod = in(i-1)  * b1_fixed(section);
            b2_prod = in(i-2)  * b2_fixed(section);
            a1_prod = out(i-1) * a1_fixed(section);
            a2_prod = out(i-2) * a2_fixed(section);
        end
        accum = b0_prod + b1_prod + b2_prod + a1_prod + a2_prod;
        
        accum = accum * 2^postShift;
        
        out(i) = fi(accum, true, 16, 15, fimath('OverflowAction', 'Saturate', 'RoundingMethod', 'Floor'));
    end
    in = out;
end

The step response is then visible on the graph above.

Edit: DC-Gain calculation:

For the double precision coefficients, the evaluation of the DC-Gain from the transfer function is: $$ H(z=1) = \frac{b_{0,1} + b_{1,1} + b_{2,1}}{1 - a_{1,1} - a_{2,1}} \cdot \frac{b_{0,2} + b_{1,2} + b_{2,2}}{1 - a_{1,2} - a_{2,2}}$$

>> dcGain_double_stage1 = (b0_double(1) + b1_double(1) + b2_double(1)) / (1 - a1_double(1) - a2_double(1))

dcGain_double_stage1 =

   0.094759343573582

>> dcGain_double_stage2 = (b0_double(2) + b1_double(2) + b2_double(2)) / (1 - a1_double(2) - a2_double(2))

dcGain_double_stage2 =

  10.553049042846906

>> dcGain_double_total = dcGain_double_stage1 * dcGain_double_stage2

dcGain_double_total =

   0.999999999999990

The same thing with convoluting the double precision coefficients:

>> b_double = conv([b0_double(1), b1_double(1), b2_double(1)], [b0_double(2), b1_double(2), b2_double(2)])

b_double =

   1.0e-03 *

   0.031238976917083   0.124955907668331   0.187433861502496   0.124955907668331   0.031238976917083

>> a_double = conv([1, -a1_double(1), -a2_double(1)], [1, -a1_double(2), -a2_double(2)])

a_double =

   1.000000000000000  -3.589733887112174   4.851275882519412  -2.924052656162454   0.663010484385890

>> dcGain_double_total1 = sum(b_double) / sum(a_double)

dcGain_double_total1 =

   1.000000000000152

Unfortunately, I have a little difference between these two calculations which definitely is due to my machine precision..

>> dcGain_double_total - dcGain_double_total1

ans =

    -1.620925615952729e-13

The DC-Gain with my fixed-point coefficients is calculated similar as the double precision one, with the difference that the coefficients have to be scaled by $2^{postShift}$ with postShift being 1:

>> dcGain_fixed_stage1 = (2^postShift * double(b0_fixed(1)) + 2^postShift * double(b1_fixed(1)) + 2^postShift * double(b2_fixed(1))) / (1 - 2^postShift * double(a1_fixed(1)) - 2^postShift * double(a2_fixed(1)))

dcGain_fixed_stage1 =

   0.094736842105263

>> dcGain_fixed_stage2 = (2^postShift * double(b0_fixed(2)) + 2^postShift * double(b1_fixed(2)) + 2^postShift * double(b2_fixed(2))) / (1 - 2^postShift * double(a1_fixed(2)) - 2^postShift * double(a2_fixed(2)))

dcGain_fixed_stage2 =

  10.568181818181818

>> dcGain_fixed_total = dcGain_fixed_stage1 * dcGain_fixed_stage2

dcGain_fixed_total =

   1.001196172248804

Same thing with convoluting my fixed-point coefficients:

>> b_fixed = conv([2^postShift * double(b0_fixed(1)), 2^postShift * double(b1_fixed(1)), 2^postShift * double(b2_fixed(1))], [2^postShift * double(b0_fixed(2)), 2^postShift * double(b1_fixed(2)), 2^postShift * double(b2_fixed(2))])

b_fixed =

   1.0e-03 *

   0.031180679798126   0.124722719192505   0.187084078788757   0.124722719192505   0.031180679798126

>> a_fixed = conv([1, - 2^postShift * double(a1_fixed(1)), -2^postShift * double(a2_fixed(1))], [1, - 2^postShift * double(a1_fixed(2)), -2^postShift * double(a2_fixed(2))])

a_fixed =

   1.000000000000000  -3.589782714843750   4.851342819631100  -2.924055889248848   0.662994079291821

>> dcGain_fixed_total1 = sum(b_fixed) / sum(a_fixed)

dcGain_fixed_total1 =

   1.001196172248804

Here the difference between the two calculations is exactly 0:

dcGain_fixed_total - dcGain_fixed_total1

ans =

     0

Please note that it is very clear to me that the frequency response of the double-precision filter and the fixed-point filter will be different. The question is about the difference of the calculated fixed-point frequency response and the simulated step response, which show different DC-Gains.

The settling value of my fixed-point filter output is:

>> double(out(end))

ans =

   0.486877441406250

This settling value corresponds to a DC-Gain of:

dcGain_fixed_step = double(out(end)) / stepHight

dcGain_fixed_step =

   0.973754882812500

And this differs from what I would expect from my calculated frequency response:

>> dcGain_fixed_step - dcGain_fixed_total

ans =

  -0.027441289436304

Answer 1

The "DC Gain" for a Filter $H(z)$ is the magnitude of the filter with $z=1$. In the z transform the unit circle ($z=e^{j\omega}$) represents the "the frequency axis" with $z=1$ as "DC" and $z=-1$ as $f_s/2$ with $f_s$ as the sampling rate. The entire frequency response is thus the magnitude and phase of $H(e^{j\omega})$.

Thus for:

$$ H(z) = \frac{b_0 z^2 + b_1 z + b_2}{z^2 + a_1 z + a_2} $$

The DC Gain will be the sum of the coefficients:

$$H(1) = \frac{b_0+b_1 + b_2}{1+a_1+a_2}$$

(Note: As Robert has pointed out in the comments, there is a factor of 2 in the signal path in the implementation in the ARM CMSIS-DSP Library. The coefficients returned by MATLAB are consistent with the form above, and would each need to be divided by two to match the implementation coefficients. As long as that is done when mapping the solution to what is implemented, we would achieve the same results as described below.)

We then when working in fixed point must be very careful with the effects of rounding error effects due to how we scale and truncate the coefficients. When working with filters in general, even with floating point, there are lots of opportunities to accumulate rounding errors in precision, particularly as the order of the filter increases (consider in general a polynomial raised the Nth order, how a very small error can be amplified!).

For example with the OP's case, the floating coefficients returned by the butter function with the OP's parameters are:

>> [b,a] = butter(n, 2 * fc / fs, 'low');

>> b
b =

   0.000031239   0.000124956   0.000187434   0.000124956   0.000031239

>> a
a =

   1.00000  -3.58973   4.85128  -2.92405   0.66301

The DC gain out to the floating point precision given is

>> sum(b)/sum(a)
ans = 1

Keep in mind that the result of this in general is complex, meaning it has both the magnitude and phase of the DC gain. In this case the magnitude is exactly 1 and the phase is exactly 0 degrees.

(Also as noted above, all of the b's, and all of the a's except the first a=1 returned by MATLAB need to be divided by two prior to scaling and quantizing due to the factor of two in the implementation datapath. This does not change what is described here but would come into play when we determine our actual coefficient values.)

If we were to manually add the coefficients as reported to the console, we would already see the effects of rounding error, reduced if we enter format long first (which only effects what is displayed, not what is actually stored in the variable but suggests being very careful with coefficient precision and being aware of opportunities where rounding may occur). For example, observe the continuing precision in the reported numerator coefficients after changing the format to long.

>> format long

b =

 Columns 1 through 4:

   3.123897691708262e-05   1.249559076683305e-04   1.874338615024957e-04   
1.249559076683305e-04

 Column 5:

   3.123897691708262e-05

For purpose of the magnitude and phase of the frequency response, there is opportunity for accumulation of rounding error (even with the floating point precision) when arriving at the numerator and denominator coefficients indirectly, such as computing the second order sections, and then from cascading those which involves convolving their coefficients get to the combined $H(z)$ function, and similarly when we consider fixed point solutions, we can get different rounding results if we quantize the resulting coefficients after the cascade versus the coefficients of each second order section. Our implementation would be done with the second order sections, so to accurately predict the second order sections we should quantize those filter coefficients specifically.

Exploring both of these points, we can use tf2sos to get the coefficients for the second order sections:

>> [sos, g] = tf2sos(b,a)
sos =

   0.000031239 0.000062478 0.000031239 1.000000000 -1.725933395 
0.747447372 

   1.000000000 2.000000034 0.999999986 1.000000000 -1.863800492  
0.887033000  

We can convolve the coefficients for the second order sections to resynthesize the combined response, and in doing this even in floating point we can start to see the introduction of rounding errors that occurred in the translation between second order and direct transfer function forms. This will be reasonably small in this case but gets more pronounced as we increase the order of filters:

>> b2 = conv(sos(1,1:3), sos(2,1:3) * g

b2 =

   0.000031239   0.000124956   0.000187434   0.000124956   0.000031239

>> b-b2
ans =

   0.0000e+00  -2.7105e-19  -7.3184e-19  -7.5894e-19  -2.2362e-19

Similarly for the denominator:

>> a2 = conv(sos(1,4:6), sos(2,4:6))
a2 =

   1.00000  -3.58973   4.85128  -2.92405   0.66301

>> a-a2
ans =

   0.0000e+00  -2.2204e-15   5.3291e-15  -2.6645e-15   3.3307e-16

Since the DC gain is the summation of the coefficients, the result of determining the DC gain after the additional processing we used will be the summation of these errors. In this case with the small order of the filter, the accumulated error while we are still in floating point is acceptably small:

>> sum(b2)/sum(a2)
ans =  1.000000000001489

Fixed Point

All that said, as soon as we go to fixed point significant rounding error will be introduced, and our primary goal is to determine the effect of this on our actual implementation. For this reason it is important to scale and quantize the coefficients exactly as will be done in implementation (so on the second order sections not the composite transfer function), and then from that determine the DC gain (or complete frequency response) using the process above. Once the coefficients of the second order sections are determined, the composite transfer function can be determined by convolving the coefficients, OR the gain of each section can be determined and the results multiplied. Both should provide the same answer and will have a consistent DC error versus the floating point result. We also see from this the benefit of “decoupling the poles” by factoring the larger polynomial into a cascade of second order sections—- you will find the sensitivity of both approaches to the error introduced in fixed point is significantly reduced as a cascade of smaller order sections. As one large polynomial, a small error gets significantly amplified by the higher order of the polynomial.