Filter State variable access with Python Wrapper

[LukeGary462]

I am in the process of prototyping some algorithms using the python wrapper for the biquad direct form two transposed filters, and notice that the state variables for a filter instance are not accessible via the wrapper.
I found this to be surprising since a list of state variables is used when initializing the filter, but according to Christophe, they are only used for initialization.

It would be useful to have access to these for things like system identification, goertzel or sliding DFT, or real-time filtering using the cmsis dsp library.

This would be a great feature to see!

below is the code I am using to test a float32 biquad for goertzel filtering, which requires state variable access for phase, and magnitude calculation. Please excuse the code, its more of a stream of consciousness 😄

BLOCK_SIZE_SAMPLES = 500
SAMPLING_FREQUENCY_HZ = 50e3
TARGET_FREQUENCY_HZ = 1.5e3

NUM_BLOCKS = 10
CSCALER = 1.0

import math
import cmsisdsp as dsp # see documentation of ARM-CMSIS_5 for installation procedure
import numpy as np
import pprint as pp
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams['figure.dpi'] = 100
plt.grid()

def calculate_coeff(Fs, Ft, N):
        ''' calculate goertzel filter coefficients to use in DF2T biquad IIR filter'''
        k = (N*Ft)/Fs
        w = 2.0*math.pi* (k/N)
        a1 = 2.0*math.cos(w)
        b1 = -1.0*math.exp(math.sin(w))
        return {
            'a': [a1, -1.0], # a0 omitted in ARM CMSIS DSP
            'b': [1.0, b1, 0.0],
            'coeffs': [1.0/CSCALER, b1/CSCALER, 0.0/CSCALER, a1/CSCALER, -1.0/CSCALER],
            'settling_time_ms': ((1.0/Fs) * N) * 1000,
            'bw': Fs/N,
            'k': k,
            'w': w,
            'states': np.zeros(2)
        }
    
goertzel_coeff = calculate_coeff(
    Fs=SAMPLING_FREQUENCY_HZ, 
    Ft=TARGET_FREQUENCY_HZ, 
    N=BLOCK_SIZE_SAMPLES
)

'''create a direct-form 2 transposed IIR filter cascade'''
goertzel = dsp.arm_biquad_cascade_df2T_instance_f32()

'''initialize the filter'''
dsp.arm_biquad_cascade_df2T_init_f32(
    goertzel,
    1, # single stage to implement goertzel
    goertzel_coeff['coeffs'],
    goertzel_coeff['states']
)

''' create some test data, two sinusoids 1kHz and 1.5kHz, and thier product. num_samples = NUM_BLOCKS x block_size '''
signal_1k0 = [np.sin(2.0*np.pi*1.0e3 *(i/SAMPLING_FREQUENCY_HZ)) for i in np.arange(BLOCK_SIZE_SAMPLES*NUM_BLOCKS)]
signal_1k5 = [np.sin(2.0*np.pi*1.5e3 *(i/SAMPLING_FREQUENCY_HZ)) for i in np.arange(BLOCK_SIZE_SAMPLES*NUM_BLOCKS)]
signal = np.asarray(signal_1k0) * np.asarray(signal_1k5)

'''actually run the filter across the product, emulating real-time sampling'''

def chunks(lst, n):
    """Yield successive n-sized chunks from lst."""
    for i in range(0, len(lst), n):
        yield lst[i:i + n]
        
signal_blocks = chunks(signal, BLOCK_SIZE_SAMPLES)
        
results = []
mags = []
for block in signal_blocks:
    results.append(dsp.arm_biquad_cascade_df2T_f32(goertzel, block))
    '''calculate magnitudes'''
    ## no access to state variables?  <-------------------------------------------------------------------
    
    pp.pprint(goertzel_coeff)
    '''reset the filter'''
    goertzel_coeff['states'] = np.zeros(2)
    dsp.arm_biquad_cascade_df2T_init_f32(
        goertzel,
        1, # single stage to implement goertzel
        goertzel_coeff['coeffs'],
        goertzel_coeff['states']
    )
    
result = []
for res in results:
    result += res.tolist()
    
'''plot raw input block'''
plt.plot(signal[0:BLOCK_SIZE_SAMPLES], label='input')
''' print the first 4 blocks'''
plt.plot(result[0:BLOCK_SIZE_SAMPLES], label='no-window')   
plt.legend(loc='best')
plt.show()

[christophe0606]

Thanks for the feedback.

I had never thought that it may be useful except for being able to call the filter several time. And, this is possible since the state is preserved and contained in the state variable arm_biquad_cascade_df2T_instance_f32.

But it is not readable from the Python.

So I need to add a Python API to enable this. I'll look at it. I cannot share any timeline because before looking at it I don't know how much time it will require.

[christophe0606]

@LukeGary462 I have pushed a new commit (d08d049) which should solve the issue.

Now you can do goertzel.state() to get the state array.

Note it is returning a copy of the internal state. You can't change the internal state array.

It is implemented only for arm_biquad_cascade_df2T_instance_f32. Other filters are not (yet) supported.

Tell me if it works.

[LukeGary462]

Thanks @christophe0606 ! I will give this a shot today

[LukeGary462]

@christophe0606 This seems to have done the trick, than you!

[christophe0606]

Don't hesitate to reopen the issue if there are problems with the solution.