search

ARM-software / CMSIS_5

CMSIS Version 5 Development Repository
http://arm-software.github.io/CMSIS_5/index.html
Apache License 2.0
1.33k stars 1.08k forks source link

arm_biquad_cascade_df2T_f32() accuracy issue #1320

Closed anhnhancao closed 3 years ago

anhnhancao commented 3 years ago

Dear all,

I'm developing a DSP project that need to design an IIR notch filter according to these coefficients value as coefficients array: notch_filter_coef[0] = 0.887839755; notch_filter_coef[1] = 1.436554901; notch_filter_coef[2] = 0.887839755; notch_filter_coef[3] = -1.436554901; notch_filter_coef[4] = -0.775679511;

However, when i call "arm_biquad_cascade_df2T_f32()" based on the above coefficients array to process input data, i did not get correct value of output data in compare to the "scipy.signal.iirnotch" in Python that i have implemented before. I got these coefficients with MATLAB and Python code based on some basic information about my notch filter (2nd order, fc = 50Hz, fs = 125, Q = 10). I import the DSP - CMSIS in STMCubeF7 (https://github.com/STMicroelectronics/STM32CubeF7/tree/master/Drivers/CMSIS/DSP)

I attach the link of images that i got results for my coefficients array.

https://scontent.fsgn5-4.fna.fbcdn.net/v/t1.15752-9/243764465_882154949360994_7550029049726538223_n.png?_nc_cat=104&ccb=1-5&_nc_sid=ae9488&_nc_ohc=P3lcTPwAenEAX9GVeBD&_nc_ht=scontent.fsgn5-4.fna&oh=82e9fd677b029330b52e66359d431153&oe=617C3334

https://scontent.fsgn5-2.fna.fbcdn.net/v/t1.15752-9/243654969_644007129915610_4194116538003354147_n.png?_nc_cat=107&ccb=1-5&_nc_sid=ae9488&_nc_ohc=luZ6d3ROCJMAX8KbLxm&_nc_ht=scontent.fsgn5-2.fna&oh=41260cb12875f794a7fab4742935f427&oe=617ACED8

I really appreciate for your help.

anhnhancao commented 3 years ago

I use this version of CMSIS - DSP

/**

christophe0606 commented 3 years ago

@anhnhancao The version 1.5.3 is quite old but the biquad df2T has not been modified for the ARM architecture used in a STMCubeF7.

So it should give the same result as a more recent version.

I have tested with a more recent version using the CMSIS-DSP Python wrapper:

import cmsisdsp as dsp
import numpy as np
from scipy import signal

# signal for tests
print("Random signal for tests")
sig=10*np.random.rand(10)-5
print(sig)

fs = 125.0  # Sample frequency (Hz)
f0 = 50.0  # Frequency to be removed from signal (Hz)
Q = 10.0  # Quality factor
# Design notch filter
b, a = signal.iirnotch(f0, Q, fs)

print("COEFS")
print(b)
print(a)
print("\n")

sos = np.hstack((b,a))
res1=signal.sosfilt(sos,sig)
print("Reference scipy")
print(res1)

numstages = 1
state=np.zeros(2*numstages)
coefs=np.reshape(np.hstack((b[:3],-a[1:])),5*numstages)

biquad = dsp.arm_biquad_cascade_df2T_instance_f32()
dsp.arm_biquad_cascade_df2T_init_f32(biquad,numstages,coefs,state)

res2=dsp.arm_biquad_cascade_df2T_f32(biquad,sig)

print("CMSIS-DSP result")
print(res2)

I am getting:

Random signal for tests
[ 3.81300478 -0.84438344  3.9325774  -3.0110439   3.98482643 -0.61073485
 -1.97300022 -1.75257331  1.78194616 -0.80930025]
COEFS
[0.88783976 1.4365549  0.88783976]
[1.         1.4365549  0.77567951]

Reference scipy
[ 3.38533723 -0.13530928  3.2322751  -2.31202232  3.51799345 -0.7515396
 -0.74038238 -3.28600826  2.60753755 -0.91165126]
CMSIS-DSP result
[ 3.385337   -0.13530916  3.2322748  -2.3120215   3.5179925  -0.751539
 -0.7403828  -3.286008    2.6075377  -0.9116513 ]

I don't see any problem.

Do you have a (short) signal I could test and which is demonstrating the problem ?

anhnhancao commented 3 years ago

I can see that there is a difference between 2 output, Let's say the output from Python code is tested and run well for many times with our application.

anhnhancao commented 3 years ago

Really appreciate for your previous answer, i will give more information on this topic:

And this is your code:

b, a = signal.iirnotch(f0, Q, fs)

print("COEFS") print(b) print(a) print("\n")

sos = np.hstack((b,a)) res1=signal.sosfilt(sos,sig)

So what is the difference between your code and mine in creating a notch filter?

christophe0606 commented 3 years ago

@anhnhancao I think the problem is that you are using filtfilt. Quote from the scipy documentation (I have put in bold):

This function applies a linear digital filter twice, once forward and once backwards. The combined filter has zero phase and a filter order twice that of the original.

scipy.signal.sosfilt is the function equivalent to CMSIS-DSP biquad.

We don't have a filtfilt in CMSIS-DSP.

Note that there is also a sosfiltfilt in scipy and we don't have it. In CMSIS-DSP we currently only have the standard single pass filtering.

We have another example of how to use the sosfilt with CMSIS-DSP here:

https://developer.arm.com/documentation/102463/latest/

but it is using DF1 Q31 biquad in this example.

anhnhancao commented 3 years ago

Thank you very much

I think I should discuss with our team to change filtfilt to sosfilt. I think maybe the effectiveness after the change won't drop much, doesn't it.

And one more question, Does CMSIS-DSP support any function like "scipy.signal.butter" and "scipy.signal.lfilter" in Python, since we still have a Band Pass filter using this API in Python?

For the "scipy.signal.butter" with second order, it returns a coefficient like below: Numerator: 0.0013487119483563447 Numerator: 0.0 Numerator: -0.0026974238967126894 Numerator: 0.0 Numerator: 0.0013487119483563447 Denominator: 1.0 Denominator: -3.8885442012895624 Denominator: 5.67620906434173 Denominator: -3.6865178016589466 Denominator: 0.8988589941552532

How can i put these value into the coefficient array of arm_biquad_cascade_df2T_f32()?

Finally, Really appreciate for your enthusiastic support.

christophe0606 commented 3 years ago

@anhnhancao You may implement the filtfilt with CMSIS-DSP using two filters. The tricky part will be the padding which must be done as in the original filtfilt.

If you use scipy.signal.butter with the option output="sos" then you get a second order sections representation of your IIR which can be used directly with the CMSIS-DSP biquads.

We don't have an equivalent to scipy.signal.lfilter because sosfilt should be used instead (it is even recommended in the documentation of scipy.signal.lfilter).

sosfilt is less sensitive to the rounding errors but otherwise is doing the same work as lfilter.

Both functions are just representing the IIR differently.

anhnhancao commented 3 years ago

We tried sosfilt for both notch filter and band pass filter, the results were quite promising for Python and CMSIS-DSP implementations. Thank you very much @christophe0606