Design filter coefficients of FIR and Biquad

Status
Not open for further replies.
W

Wicky21

Well-known member
Can somebody please tell me how to obtain Coefficients for IIR and FIR filters?? As noted in the FIR filter design tool, "note section", TFilter link is not properly working.Also if using MATLAB fdatool i don't know how to obtain Coefficients properly on Teensy 3.2. Somebody check this and let me know.
In fdatool, there's something called "Quantization" parameters section. in there we can set filter arithmetics as fixed point. Also something called "word length" and some other parameters. I'm not sure what to select and how to obtain proper Coefficients for teensy 3.2.
 
Last edited:
In File --> Examples --> Audio --> Effects --> Filter_FIR
There's an example. Does anybody know how these coefficients obtained ?
 
When using State variable filters and Biquad filters how to adjust the pass band range?
In biquad => setBandpass(stage,frequency,Q) --> frequency is just a one value. how to adjust a range?
In filter =>frequency(freq) --> Only can give one value for bandpass also
 
Wicky21 said:
Can somebody please tell me how to obtain Coefficients for IIR and FIR filters
Click to expand...

I made this small matlab program to plot an print FIR coefficients (here for single float and 32 bit integer) ready to be cut and pasted to header file
Code: 
nt = 128;
fsh=375/2;
f1=20; f2=120; f3=140;
h1=fir1(nt,f1/fsh)'; 
h2=fir1(nt,f2/fsh,'high')'; 
h3=fir1(nt,[f2,f3]/fsh)'; 
h4=fir1(nt,[f1,f2]/fsh,'stop')'; 
h5=fir1(nt,[f1,f2,f3]/fsh,'DC-1')';
yy=fft([h1,h2,h3,h4,h5],1024); 
figure(1)
plot((0:511)'/1014*375,20*log10(abs(yy(1:512,:))))
title(sprintf('fc = [%d, %d, %d]',f1,f2,f3))

fprintf('%17.10ef, %17.10ef, %17.10ef, %17.10ef, %17.10ef, %17.10ef,\n',h5)
fprintf('\n');
fprintf('%12d, %12d, %12d, %12d, %12d, %12d,\n',round(h5*(2^31-1)))
fprintf('\n');

adapt to your frequency and filter requirements

here comes a program that generates bi-quad coefficients for 'direct Form 1'
Code: 
xx=randn(10000,1);

f1=20; f2=120; f3=140;
fsh=375/2;
%[b,a]=butter(4,[0.1 0.2],'bandpass');
%[z,p,k]=butter(4,[0.1 0.2],'bandpass');
%[z,p,k]=cheby1(4,0.5,[120 140]/375*2,'bandpass');
if 0
[z1,p1,k1]=cheby1(8,0.5,f1/fsh,'low');
[z2,p2,k2]=cheby1(8,0.5,f3/fsh,'high');
[z3,p3,k3]=cheby1(4,0.5,[f2,f3]/fsh,'bandpass');
end
if 0
[z1,p1,k1]=cheby2(8,80,f1/fsh,'low');
[z2,p2,k2]=cheby2(4,80,[f2,f3]/fsh,'bandpass');
end
if 1
[z1,p1,k1]=butter(8,f1/fsh,'low');
[z2,p2,k2]=butter(8,f3/fsh,'high');
[z3,p3,k3]=butter(4,[f2,f3]/fsh,'bandpass');
end
[s1,g1]=zp2sos(z1,p1,k1);
[s2,g2]=zp2sos(z2,p2,k2);
[s3,g3]=zp2sos(z3,p3,k3);

yy(:,1)=g1*sosfilt(s1,xx);
yy(:,2)=g2*sosfilt(s2,xx);
yy(:,3)=g3*sosfilt(s3,xx);
figure(1)
pwelch(yy);

s1x=s1(:,[1,1,2,3,5,6]); s1x(:,2)=0;
s2x=s2(:,[1,1,2,3,5,6]); s2x(:,2)=0;
s3x=s3(:,[1,1,2,3,5,6]); s3x(:,2)=0;

lp=fix((s1x*((2^15)-1)/2))
hp=fix((s2x*((2^15)-1)/2))
bp=fix((s3x*((2^15)-1)/2))
to understand the '0's please consult the CMSIS manual e.g. https://www.keil.com/pack/doc/CMSIS...deDF1.html#gad54c724132f6d742a444eb6df0e9c731

This code is not suited for audio library biquad but is for CMSIS biquad
 
Last edited:
Wicky21 said:
Hello WMXZ,
Thanks for the code. But code is not much clear. How to adjust number of taps?Cutoff(bandpass,lowpass,hightpass)?? can u comment it and post it back? Also have you tried with fdatool of matlab? it's much easier. I prefer coefficients for Teensy 3.2
Click to expand...

Sorry, as you mentioned fdatool of matlab, I assumed you are used to program matlab.
the coefficients generated where used in Teensy 3.2 and for floating pint in Teensy 3.6

to change the parameter: hint: try and read the program, in the end you have to justify your choices.
 
Last edited:
Wicky21 said:
I want to find out coefficients of FIR bandpass filter,
lower cutoff = 350Hz, Upper cutoff = 2000Hz, Order = 100
Click to expand...

for q15 and fs 44.1 kHz:

Code: 
nt=100;
f1=350;
f2=2000;
fsh=44100/2;
hh=fir1(nt,[f1,f2]/fsh)'; 
ihh=round(hh*2^15-1)
 
I want to find out coefficients of FIR bandpass filter,
lower cutoff = 350Hz, upper cutoff = 2000Hz, Order = 100
Should design this filter teensy 3.2
What is "'fsh" ?
 
WMXZ said:
for q15 and fs 44.1 kHz:

Code: 
nt=100;
f1=350;
f2=2000;
fsh=44100/2;
hh=fir1(nt,[f1,f2]/fsh)'; 
ihh=round(hh*2^15-1)
Click to expand...

Here we go. now i understand. This is the basic structure how design coefficients. When your rounding up coefficients are you sure response want get change??? No need of quantization?

Also i need sample program for Biquad filter also. Same bandpass with cutoff frequencies. But with 4 stages( biquad cascaded filters)
 
Wicky21 said:
When your rounding up coefficients are you sure response want get change??? No need of quantization? It is with quantization! (2^15 for q15 or 1^31 for q31)

Also i need sample program for Biquad filter also. Same bandpass with cutoff frequencies. But with 4 stages( biquad cascaded filters)
Click to expand...
it is in second program I gave in post #4

eg:
Code: 
ns=4; % number of stages
[z3,p3,k3]=butter(ns,[f1,f2]/fsh,'bandpass');
[s3,g3]=zp2sos(z3,p3,k3);
s3x=s3(:,[1,1,2,3,5,6]); s3x(:,2)=0;
bp=round((s3x*((2^15)-1)/2))

you need to read CMSIS documentation on division by 2

Edit: this code is not suited for audio library
 
Last edited:
Wicky21 said:
Thanks alot. I check this things on teensy 3.2 and let you know if something wrong
Click to expand...

if you generate some random noise, you can assess the filter quality by looking into the spectrum before and after the filter
 
WMXZ said:
if you generate some random noise, you can assess the filter quality by looking into the spectrum before and after the filter
Click to expand...

Ok. Here WMXZ, Can you post some arduino code which was written for Biquad filter which has 4 stages with above requirements. Because hard to find out a example code. It will be a great help.

In audio system design tool simply needed to Get input from
ADC object --> biquad1 --> usb1 + dac1
 
Wicky21 said:
In audio system design tool simply needed to Get input from
ADC object --> biquad1 --> usb1 + dac1
Click to expand...

If you wanted to use the audio library, did you read the docu on biquad? esp. suggesting not to use lower frequencies below 400 Hz?
 
Is this an actual project (which really gets built with a Teensy), or a school or academic exercise that doesn't actually get implemented?
 
PaulStoffregen said:
Is this an actual project (which really gets built with a Teensy), or a school or academic exercise that doesn't actually get implemented?
Click to expand...
This is an actual implemented project. Paul just help me with figuring out proper coefficients. As WMAX mentioned do you have any other solution ? Can you post some example code?
 
My requirement is to make an IIR bandpass filter(which has higher order--> with four stages) which has lower cutoff frequency of 350Hz and upper cutoff frequency at 2000Hz by using teensy 3.2. Also frequncy response which has butter worth response.
 
Wicky21 said:
My requirement is to make an IIR bandpass filter(which has higher order--> with four stages) which has lower cutoff frequency of 350Hz and upper cutoff frequency at 2000Hz by using teensy 3.2. Also frequncy response which has butter worth response.
Click to expand...

here is a test case where you can see the filter in action
Transfer is given by 4 identical bandpass filters.

Code: 
#include "Audio.h"

// GUItool: begin automatically generated code
AudioInputAnalog         adc1;           //xy=149,215
AudioSynthNoiseWhite     noise1;         //xy=142,315
AudioFilterBiquad        biquad1;        //xy=337,418
AudioAnalyzeFFT1024      fft1024;        //xy=530,315
AudioConnection          patchCord2(noise1, biquad1);
AudioConnection          patchCord3(biquad1, fft1024);
// GUItool: end automatically generated code

// An array to hold the 16 frequency bands
float level[16];

void setup() {
  // Audio requires memory to work.
  AudioMemory(12);

  float fc = (2000.0+350.0)/2.0;
  float df = (2000.0-350.0);
  AudioNoInterrupts();
  for(int ii=0; ii<4;ii++)
  { biquad1.setBandpass(ii,fc,fc/df);
  }
  noise1.amplitude(0.7);

  AudioInterrupts();

}


void loop() {
  if (fft1024.available()) {
    // read the 512 FFT frequencies into 16 levels
    // music is heard in octaves, but the FFT data
    // is linear, so for the higher octaves, read
    // many FFT bins together.
    level[0] =  fft1024.read(0);
    level[1] =  fft1024.read(1);
    level[2] =  fft1024.read(2, 3)/(1+3-2);
    level[3] =  fft1024.read(4, 6)/(1+6-4);
    level[4] =  fft1024.read(7, 10)/(1+10-7);
    level[5] =  fft1024.read(11, 15)/(1+15-11);
    level[6] =  fft1024.read(16, 22)/(1+22-16);
    level[7] =  fft1024.read(23, 32)/(1+32-23);
    level[8] =  fft1024.read(33, 46)/(1+46-33);
    level[9] =  fft1024.read(47, 66)/(1+66-47);
    level[10] = fft1024.read(67, 93)/(1+93-67);
    level[11] = fft1024.read(94, 131)/(1+131-94);
    level[12] = fft1024.read(132, 184)/(1+184-132);
    level[13] = fft1024.read(185, 257)/(1+257-185);
    level[14] = fft1024.read(258, 359)/(1+359-258);
    level[15] = fft1024.read(360, 511)/(1+512-360);
    // See this conversation to change this to more or less than 16 log-scaled bands?
    // https://forum.pjrc.com/threads/32677-Is-there-a-logarithmic-function-for-FFT-bin-selection-for-any-given-of-bands

    for (int ii=0; ii<16; ii++) {
      Serial.printf("%4d",(int) (1000*level[ii]));
    }
    Serial.println();
    delay(100);
  }
}

you could set the coefficients yourself, but you have to understand the biquad implementation in the audio library for that.

I let the interpretation of the output (log bin to frequency) for you as exercise.
 
In audio library biquad implementation is done by using Direct form 1 structure or Direct form 11 transposed structure?
According to matlab code you have sent, biquad coefficients can be obtained right. Without getting mad on me can you post how these coefficients can be added to the arduino code? because i guess we need to convert the data into "Short" values right? (Which was mentioned in FIR note section)
 
Wicky21 said:
In audio library biquad implementation is done by using Direct form 1 structure or Direct form 11 transposed structure?
According to matlab code you have sent, biquad coefficients can be obtained right. Without getting mad on me can you post how these coefficients can be added to the arduino code? because i guess we need to convert the data into "Short" values right? (Which was mentioned in FIR note section)
Click to expand...
The structure of the biquad implementation is not clear to me, it is highly optimized for 16 bit arithmetic , maybe Paul can give a hint to the source.

conversion to short is simply done by multiplying real number with ((1<<15) -1) or (2^15-1) and taking the integer value (or rounding it to next integer)
 
WMXZ said:
maybe Paul can give a hint to the source.
Click to expand...

It's here:

https://github.com/PaulStoffregen/Audio/blob/master/filter_biquad.cpp#L32

Input and output are 16 bits. Internal state and coefficients are 32 bits.

conversion to short is simply done by multiplying real number with ((1<<15) -1) or (2^15-1) and taking the integer value (or rounding it to next integer)
Click to expand...

Coefficients use 1.31 format. Look at these conversion functions for examples.

https://github.com/PaulStoffregen/Audio/blob/master/filter_biquad.h#L42
 
WMXZ said:
The structure of the biquad implementation is not clear to me, it is highly optimized for 16 bit arithmetic , maybe Paul can give a hint to the source.

conversion to short is simply done by multiplying real number with ((1<<15) -1) or (2^15-1) and taking the integer value (or rounding it to next integer)
Click to expand...

Ok. according to matlab code you have sent coefficient results were as below
Stage 0 = 16384 0 32767 16384 -27508 11754
Stage 1 = 16384 0 32767 16384 -29313 14093
Stage 2 = 16384 0 -32767 16384 -31061 14748
Stage 3 = 16384 0 -32767 16384 -32270 15929

It has 6 numbers. According to biquad function --> setCoefficients(stage, array[5]); coefficients in order with B0, B1,B2,A1,A2.
So according to output results of matlab lets say Stage 0, B0=16384, B1 = 32767, B2 = 16384 , A1 = -27508 , A2 = 11754. Am i correct ? What is these "16384" value? Is this a Int value? Do we need to convert these into some other format? Or can directly assign this output values?
 
Paul check with the code which WMAX posted how to obtain coefficients.

ns=4; % number of stages
f1=350;
f2=2000;
fsh=44100/2;

[z3,p3,k3]=butter(ns,[f1,f2]/fsh,'bandpass');
[s3,g3]=zp2sos(z3,p3,k3);
s3x=s3:),[1,1,2,3,5,6]); s3x:),2)=0;

bp=round((s3x*((2^15)-1)/2))


If this code is correct simply what we have to do is
ex:-
Stage 0 = 16384 0 32767 16384 -27508 11754

Stage 0, B0=16384*(1073741824.0),
B1 = 32767*(1073741824.0),
B2 = 16384 *(1073741824.0),
A1 = -27508* (1073741824.0),
A2 = 11754 * (1073741824.0)

Am i right ?
 
Status
Not open for further replies.
Back
Top