CMSIS DSP Biquad fixpoint error on Cortex M4

[E-T-R]

I have a problem by implementing a CMSIS Biquad lowpass filter an a Cortex M4 from TI.

The CMSIS FLOAT filter is working but not the fixpoint Q15 one.

It's just a simple 2nd order (1 stage) filter. The coeffs are calculated with Matlab.

Did anybody tried this allready and see an error in here?

The chrip function generate a chirp sinus with an amplitude of +-1.

Thanks in advance

René

// fs = 16000;
// fg = 1600 Hz, Butter
// [z,p,k] = butter(2,fg/fs,'low');
// sos2 = zp2sos(z,p,k);
// Hd = dfilt.df2tsos(sos2);
// fvtool(Hd,'Analysis','freq')
// coef = coeffs(Hd);
// coef.SOSMatrix
// ans = 0.0201 0.0402 0.0201 1.0000 -1.5610 0.6414

#define STAGES ((uint8_t)1)

arm_biquad_cascade_df2T_instance_f32 IIR; //stucture for CMSIS-DSP IIR DF2
arm_biquad_casd_df1_inst_q15 IIR_16;

// a11 und a12 müssen negiert werden
float fCoeffs[] = { 0.0200833659619093, 0.0401667319238186, 0.0200833659619093, 1.56101810932159, -0.641351521015167 };
float32_t IIRstate[2*STAGES];
BiquadLP BiQuadFloat(STAGES, fCoeffs, IIRstate); // this calls arm_biquad_cascade_df2T_init_f32(&IIR, numStages_, pCoeffs_, pState_);

// a11 und a12 müssen negiert werden
// b10, 0, b11, b12, a11, a12
q15_t n16Coeffs[STAGES * 6] = {329, 0, 658, 329, 25576, -10508};
q15_t n16StateA[STAGES * 4];
BiquadLP BiQuadQ15 = BiquadLP(STAGES, n16Coeffs, n16StateA); // this calls arm_biquad_cascade_df1_init_q15(&IIR_16, numStages_, pCoeffs_, pState_, 1);

float fFilterIn = chirp.calc_log();
fFilterOut1 = BiQuadFloat.filter(&fFilterIn);q15_t n16FilterIn = (q15_t)(fFilterIn * 16384.0f); // only scaled by Q14 for test
n16FilterOutQ15 = BiQuadQ15.filter(&n16FilterIn);

---

And this is what I get. (saved in an ring buffer)

The first chart is the chirp signal (fFilterIn).

The 2nd is the output of the float filter (fFilterOut1).

The 3rd is the output of the fixpoint filter (n16FilterOutQ15).

[llefaucheur]

Hi Rene, "BiQuadQ15" is not part of CMSIS-DSP, but is related to TI implementation. Can you share the related documentation (URL for example). Regards, Laurent.

[E-T-R]

Hi Laurent
I upload my biquad class
https://github.com/E-T-R/CMSIS-fixpoint-biquad-test
Kind regards
René

[llefaucheur]

Hi René, sorry I am not familiar with C++ syntax. I noticed "class BiquadF32", followed by fixed-point subroutine : arm_biquad_casd_df1_inst_q15. Could we have a mix of q15 and f32 here ? More specifically can you put a breakpoint at the time you enter CMSIS subroutine and give the parameter details ? Kind regards, Laurent.

[E-T-R]

Hello Laurent
I checked the memory for that idea and it looks ok to me. The pointer for the coeffs points to the right memory adress. The compiler choose the right function (because of the arguments).
Kind regards
René

[christophe0606]

Hello @E-T-R ,

A very experimental Python wrapper was recently introduced for the CMSIS-DSP. It is a good way to experiment and explore that kind of problems. Of course, the wrapper is not using the optimized version of the CMSIS-DSP so in case of bug there, it won't be visible from the wrapper.

Let's reproduce your example with scipy and the cmsisdsp module. It will be very close to MATLAB and the C API of the CMSIS.

First the standard import:

import cmsisdsp as dsp
import numpy as np
from scipy import signal
from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, grid, axes, show,semilogx, semilogy

Now, let's reproduce your matlab code to have the same coefficients.

fs = 16000
fg = 1600
[z,p,k] = signal.butter(2,fg/fs,btype='low',output='zpk')
sos2 = signal.zpk2sos(z,p,k)

You can check by printing sos2 that the same coefficients were computed

print(sos2)

gives:

[[ 0.02008337 0.04016673 0.02008337 1. -1.56101808 0.64135154]]

Now, let's generate a signal. I don't have your original signal but I've tried to generate something similar.

nb = 1000
T = 0.01
t = np.linspace(0, T, T*fs, endpoint=False)
sig = signal.chirp(t,f0=1000,f1=2000,t1=T)

Let's filter the signal with scipy and display the input and output of the filter:

res=signal.sosfilt(sos2,sig)
figure()
plot(sig)
figure()
plot(res)
show()

The input

The filtered signal

Now let's do the same with CMSIS-DSP using Q15 and df1 biquad filter.

Two things to notice:

• The coefficients are using a different convention in CMSIS-DSP for biquad Q15 as can be seen from the doxygen documentation
• Correct scaling of coefficients is important or there may be problems.

For first point, the a coefficient are the negative in CMSIS-DSP compared to Matlab / Scipy. The 1.0 coefficients (first 'a' coefficient) is ignored. And there is an additional 0 coefficient after the b10 coefficient.

Which means that you can't use the coefficients directly from sos2.
The following line is doing the conversion

numStages=1
coefs=np.reshape(np.hstack((sos2[:,:1],[[0.0]],sos2[:,1:3],-sos2[:,4:])),numStages*6)

Then those coefficients must be scaled and converted to Q15. After some experiment with this Python API and with my input signal, I have found that I need at least to scale by 2 the coefficients to get the right result. So the CMSIS-DSP postfix value must be set to 1 to compensate.

coefs = coefs / 2.0
coefsQ15 = toQ15(coefs)
postshift = 1

toQ15 is a small utility function I have written to convert to Q15. I also have the Q15toF32 to convert back and display the result.

Now, we are ready to use the CMSIS-DSP biquad. The Python API is very close to the C one

We create the instance variable

biquadQ15 = dsp.arm_biquad_casd_df1_inst_q15()

We initialize the instance

numStages=1
state=np.zeros(numStages*4)
dsp.arm_biquad_cascade_df1_init_q15(biquadQ15,numStages,coefsQ15,state,postshift)

We filter the signal

sigQ15=toQ15(sig)
res2=dsp.arm_biquad_cascade_df1_q15(biquadQ15,sigQ15)

Now we can convert it back to f32 and display the result

res2=Q15toF32(res2)
figure()
plot(res2)
show()

We have the same result than the scipy one:

So, I think your problem should be solved if you save the coefficients with the right CMSIS-DSP conventions and if you rescale them. Your scaling may be different from mine.

I hope it helps.

[E-T-R]

Hi Christophe

Thank you very much for this explanation.

What are your coefsQ15?
What is your toQ15() function doing - Is it just converting to integer without rounding?

Kind regards
René

[christophe0606]

@E-T-R toQ15 is the most simple (and wrong) way to convert coefficients. It is not doing any rounding. Just converting to integer (which is too simple so wrong). But I was not trying to design a filter. Just showing an example of the use of the API.

I no more have the coefsQ15 I used for my previous answer.

Q15 BiQuad may have a slight residual offset in the output depending on your coefficients. It is something that we will have to improve.