I’m trying to implement a bandpass biquad filter in photon. The filter is designed in Matlab and the coefficients are passed to the filter. This introduces a lot of noise. I’m guessing that is due to the overflow. Is there any easy way to debug this?
Hi @ransinha
Without seeing your actual coefficients, I can’t really say but it is strange that you subtract the feed-forward section from the feedback section rather than add. It depends on the signs in the denominator. If you give us a bit more info like the actual coefficients and the SHIFT_POINT, we can help more.
I have access to MATLAB too and I would simulate the integer version there before going to the Particle device to make sure it works. It is easy to design a filter that is stable in doubles but not when scaled to integers.
Normally the numerator and denominator need different scaling (SHIFT_POINT in your code) since the numerator is often [0.5 1 0.5] or similar and for stability the denominator needs more integer precision.
Have you considered just doing the filter in double precision and then casting back to int16? It will be slower of course but it would be a good baseline step to help you debug. You only need the feedforward and backward delays in doubles, not the entire array, but you will need to write your code a bit differently.
1 Like
Thank you so much for the reply.
num =
0.0294 0.0294 0.0290 0.0290 0.0288
0 0 0 0 0
-0.0294 -0.0294 -0.0290 -0.0290 -0.0288
den =
-1.9637 -1.9854 -1.9335 -1.9604 -1.9368
0.9756 0.9881 0.9425 0.9639 0.9424
SHIFT_POINT = 14
The code written for matlab works just fine, even the C++ counterpart that I ran on my local seems to work fine.
As i need a bandpass filter, I’m using the following to design the filter:
d = fdesign.bandpass(‘Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2’,low_stop/ny, low_pass/ny, high_pass/ny, high_stop/ny,Astop1, Apass, Astop2)
biquad = design(d,‘butter’, ‘FilterStructure’,‘df1sos’)
It gives me a sosMatrix and a set of sclValues, that is used to calculate the num, den.
That’s not very many decimal places… I know you said you ran a local C++ version, but have you looked at the precision of your matlab printout?
i.e.
format long
numAre you saying, the coefficients are not correct?
I am saying that you may be using truncated values that do not do what you expect
For instance, the numerator values, which are normally [+1,0,-1] for Butterworth, you have rescaled to [0.0294 0 -0.0294] or [482 0 -482] when scaled by your SHIFT_POINT. This changes your signal-to-noise ratio in an unfavorable way. This is why I said above that the numerator and denominator are often scaled independently.
Can you please say what your design parameters are? That is the band edges and attenuations (low_stop/ny to Astop2 in the above).
sampling_rate = 32000;
low_pass = 300;
low_stop = 200;
high_pass = 500;
high_stop = 800;
ny = sampling_rate / 2;
Apass = abs(db(.98)) ;
Astop1 = abs(db(.05)) ;
Astop2 = abs(db(.05)) ;
Hi @ransinha
I looked your code over again and I didn’t see anything that is really wrong. You are multiplying by the center numerator coefficient which is always zero, so you can remove that for increased performance with these coefficients.
I did redesign your filter in MATLAB using the FDATool option to reorder and scale second-order section filters using the L2 norm (the setting in between avoid overflows and maintain SNR). This makes all of the scale values 1.0 except the first one:
SOS =
Columns 1 through 3
0.0135523894497101 0 -0.0135523894497101
0.0636935818924626 0 -0.0636935818924626
0.0224460111104291 0 -0.0224460111104291
0.0464504698939043 0 -0.0464504698939043
0.911807277234938 0 -0.911807277234938
Columns 4 through 6
1 -1.93677219228455 0.942386590630581
1 -1.98537817946402 0.988146918670487
1 -1.96374900333443 0.975591650233341
1 -1.96037733143248 0.963900018891779
1 -1.93346685873693 0.942511644255031
>> G
G =
0.0254164178771157
1
1
1
1
1
Based on your 32-bit coefficient data type I would recommend a SHIFT_POINT much larger than what you are using now. Your current value of 14 implies 18 integer bits and 14 fractional bits in your 32-bit word. I would try something like 20 to 28 for the binary point.
Are your input values really 16-bit or are you reading them from an analog input and they are 12-bits of precision?
2 Likes
We are using the photon’s ADC output as the input to the filter, which gives out unsigned integer of 12 bit precission. I have used an uint16_t buffer for the ADC output. The filter uses int16_t data. I also think there lies the problem.
So your data is intrinsically unsigned as it is now. Does that represent your actual data or do you intend for your data to be centered on 2048? If you want to make it signed and bipolar, you can subtract 2048 from your points to “center” it on the mid-rail value of the ADC.
There are certainly optimizations you could do based on the fact that you only have 12-bit data, but it doesn’t really change things.
Maybe it is time to try the double precision version I suggested earlier?