arm_biquad_cascade_df1_init_q15() - how to properly specify/scale coefs?

[tdjastrzebski]

I try to put a simple 100Hz low-pass filter using coefs calculated with Nigel Redmon's online biquad calculator.
For 10kHz sample rate and Q=1 coefs are:
a0=0.0009566029866080141 (becomes b0)
a1=0.0019132059732160282 (b1)
a2=0.0009566029866080141 (b2)
b0=1.0 (-a0 - not used)
b1=-1.9352943868599919 (-a1)
b2=0.9391207988064236 (-a2)

This calculator uses a different convention, "a" values have to be swapped with "b" values and sign needs to be changed, but other than that, when used with arm_biquad_cascade_df2T_init_f32() as coefs, they work just great. No issues at all.

Now, instead of arm_biquad_cascade_df2T_init_f32() I try to use arm_biquad_cascade_df1_init_q15() function.
I normalize the above coesfs using this code:

double biquadCoeffs[5];
...
double maxValue;
arm_absmax_no_idx_f64(biquadCoeffs, 5, &maxValue);
for (int8_t i = 0; i < 5; i++) {
   biquadCoeffs[i] = (biquadCoeffs[i] / maxValue) * 32767.0;
}
q15_t coeffs_q15[6] = {(q15_t)biquadCoeffs[0], 0, (q15_t)biquadCoeffs[1], (q15_t)biquadCoeffs[2], (q15_t)biquadCoeffs[3], (q15_t)biquadCoeffs[4]};

which brings the following values: {16, 0, 32, 16, 32767, -15900} which, when supplied as arm_biquad_cascade_df1_init_q15() coefs, do not yield the expected results. The output signal has much smaller amplitude (1:285), but it seems like it is not filtered at all.
What am I missing? Are coefs need to be normalized somehow differently? Do they need to be anchored at a0? If so, at what value? I found no answer in docs.

[christophe0606]

@tdjastrzebski You need to use the postShift argument in arm_biquad_cascade_df1_init_q15 and only scale your coefficients by power of 2.
In your example, you probably only need to shift the coefficient by 1 and use a postShift of 1.

[christophe0606]

@tdjastrzebski

Here is an example with a postshift of 1 using Python package of CMSIS-DSP.

I have named the coefficients like above but in CMSIS-DSP and scipy the naming a,b is the opposite.

import cmsisdsp as dsp
import cmsisdsp.fixedpoint as f
import numpy as np
from scipy import signal
from pylab import figure, plot, show

a=[0.0009566029866080141
,0.0019132059732160282
,0.0009566029866080141] # zeros (numerator)
b=[1.0
,-1.9352943868599919
,0.9391207988064236] # poles (denominator)

z,p,k = signal.tf2zpk(a,b)
sos = signal.zpk2sos(z,p,k)

test_length_seconds = 0.1
signal_frequency = 100
sampling_freq = 8000
nbSamples = int(test_length_seconds*sampling_freq)
sig = 0.5*np.sin(2*np.pi*signal_frequency*(1+np.linspace(-0.2,0.2,nbSamples))*np.linspace(0,test_length_seconds,nbSamples))


ref=signal.sosfilt(sos,sig)
plot(ref)
show()

biquadQ15 = dsp.arm_biquad_casd_df1_inst_q15()

numStages=1
state=np.zeros(numStages*4)

# Convert to CMSIS-DSP format
# b10, 0, b11, b12, a11, a12
# where b is numerator (so naming different from above)
coefs=np.reshape(np.hstack((sos[:,:3],-sos[:,4:])),5*numStages)
coefs=np.hstack((coefs[0],0,coefs[1:]))

# Rescale coefficients so that they are all < 1
coefs = coefs / 2.0 
# postshift factor used in algorithm to compensate the
# coefficient shift
postshift = 1
# Convert to Q15
coefsQ15 = f.toQ15(coefs)
dsp.arm_biquad_cascade_df1_init_q15(biquadQ15,numStages,coefsQ15,state,postshift)
sigQ15=f.toQ15(sig)

resQ15=dsp.arm_biquad_cascade_df1_q15(biquadQ15,sigQ15)

res=f.Q15toF32(resQ15)
plot(res)
show()

[tdjastrzebski]

@christophe0606 Thank you for this example. I think it is time for me to learn Python:)
Anyway, the used Python functions are well documented, so I should be able to translate it to C++.

[tdjastrzebski]

@christophe0606, Is the way q15 coefs need to be scaled documented somewhere or it is just an industry standard?
If so, where/how is it defined? I need to update my references.

[christophe0606]

@tdjastrzebski On this page : https://arm-software.github.io/CMSIS-DSP/latest/group__BiquadCascadeDF1.html

Look for "scaling of coefficients"

[tdjastrzebski]

@christophe0606 Now I understand that what I was really missing was this transformation:

z,p,k = signal.tf2zpk(a,b)
sos = signal.zpk2sos(z,p,k)

Are you aware of any C++ implementation? I can calculate it in Python, but I may need to adjust Fc and Q dynamically starting with Nigel Redmon's or Robert Bristow-Johnson's formulae.

[christophe0606]

@tdjastrzebski Those functions are only useful if you need to express a high degree rational fraction as a product of biquads (second order sections).

In your example you start with degree 2 numerators and denominators so those functions are doing nothing.
If I print sos I get:

[[ 9.56602987e-04  1.91320597e-03  9.56602987e-04  1.00000000e+00
  -1.93529439e+00  9.39120799e-01]]

And this is your a and b coefficients.

The functions will be useful if you need more than one biquad in your implementation (for more complex filters).

[llefaucheur]

A general comment about your filter design. There are good reasons for cascading several biquads especially in your example: the transfer function of the filter without the numerator (plot below) gives a +50dB amplification. It means the truncation effect, after the sums from the direct path and feedback, will be amplified by 50dB in the recursive path of the filter. If for example your input data is -20dB below full-scale, it means you will have rougly 30dB SNR on the output (rule of thumb : full-scale(96dBm0) -20dB = 76dBm0 the noise is +50dB, your data will have a 76-50=26dB SNR, from an helicopter view).
When you design such a sharp low-pass filter, you may find some value in splitting the filtering effect with several biquads in cascade and have less gain in each of the recursive parts.
(to put in our low-priority-to-do-list : add a DF1-Q15 new filter structure with ESS / error-spectral-shaping to improve SNR)

[tdjastrzebski]

@llefaucheur, @christophe0606 Thank you, I think now I better understand q15 limits.
The reason why I would like to use q15 vs f32 or f64 is performance and limited RAM - it is an embedded STM32U5 device with 14bit oversampled ADC, 16bit effective. If 2 or 3 cascaded filters may be required to accomplish my goal, then the choice between fixed point and float is no longer so obvious.

[llefaucheur]

Hi Thomas, if your problem is to decimate then filter to extract a signal, I would propose you to use CIC filters (just add/sub, no coefficient) which are commonly used to decimate PDM or sigma-delta oversampled streams in A/D converters. Your last processing stage will be the IIR filter, but the processing will be on a decimated stream at FS/16 (for example) and consumming less CPU. Regards, Laurent.
https://en.wikipedia.org/wiki/Cascaded_integrator%E2%80%93comb_filter

[tdjastrzebski]

@llefaucheur Thank you, I also explore what STM32U5xx MPU has to offer, and there are three choices:

• Multi-function digital filter (MDF)
• Filter math accelerator (FMAC)
• Audio digital filter (ADF)

The only problem is that these are quite new interfaces and almost no examples/application notes exist yet.

[llefaucheur]

ADF and MDF can be configured as Sinc4/5 (this is the same as "CIC 4/5th-order"). The Application Note "AN5305" is for FMAC.

[christophe0606]

I close this ticket now since I think the CMSIS-DSP part has been covered.