arm_biquad_cascade_df1_init_q15() - how to properly specify/scale coefs? about cmsis-dsp HOT 13 CLOSED

@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()

from cmsis-dsp.

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

Look for "scaling of coefficients"

from cmsis-dsp.

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

from cmsis-dsp.

@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.

from cmsis-dsp.

@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++.

from cmsis-dsp.

@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.

from cmsis-dsp.

@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.

from cmsis-dsp.

@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).

from cmsis-dsp.

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)

from cmsis-dsp.

@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.

from cmsis-dsp.

@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.

from cmsis-dsp.

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.

from cmsis-dsp.

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

from cmsis-dsp.