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