Design filter coefficients of FIR and Biquad

Status
Not open for further replies.
Wicky21 said:
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 ?
Click to expand...

No it not correct:
Paul uses as coefficients (int) 2^31 and not short (2^15)
1073741824.0 = 2^30 = (2^31) /2
OK we can argue if we should subtract 1 to limit positive numbers or not.

it would be bp = s3x *2^30;
but better send to biquad.coefficients 's3x' as this function also allows float to be passed, as indicated by the documentation and the code you looked up.

the full code would be

Code: 
[z3,p3,k3]=butter(ns,[f1,f2]/fsh,'bandpass');
[s3,g3]=zp2sos(z3,p3,k3);
s3x=s3(:,[1,2,3,5,6]); % NOTE: this is what audio biquad expects (it does ot want a0 only: b0,b1,b2,a1,a2)

bp=floor((s3x*2^30)) % make sure no overflow by using floor

Note: that this code is not consistent with coefficient generation code used by Paul.
have not figured out why
 
According to Matlab code
bp=floor((s3x*2^30)); outputs below results

1073741824 2147483648 1073741824 -1802809978 770344455
1073741824 2147483648 1073741824 -1921103689 923625477
1073741824 -2147483648 1073741824 -2035696176 966563104
1073741824 -2147483648 1073741824 -2114896261 1043937547

But since internally multiply by 2^30,
coef[0] = coefficients[0] * 1073741824.0;
coef[1] = coefficients[1] * 1073741824.0;
coef[2] = coefficients[2] * 1073741824.0;
coef[3] = coefficients[3] * 1073741824.0;
coef[4] = coefficients[4] * 1073741824.0;

enough to take s3x as below right?

1 2 1 -1.6790 0.7174
1 2 1 -1.7892 0.8602
1 -2 1 -1.8959 0.9002
1 -2 1 -1.9697 0.9722

I wrote some code on arduino. But there will be errors. Please check and send me the proper code

Code: 
#include <Audio.h>
#include <Wire.h>
#include <SPI.h>
#include <SD.h>
#include <SerialFlash.h>

const float BP_coeff_stag0 [5]{1.0000, 2.0000, 1.0000, -1.6790, 0.7174 };
const float BP_coeff_stag1 [5]{1.0000, 2.0000, 1.0000, -1.7892, 0.8602 };
const float BP_coeff_stag2 [5]{1.0000, -2.0000, 1.0000, -1.8959, 0.9002 };  
const float BP_coeff_stag3 [5]{1.0000, -2.0000, 1.0000, -1.9697, 0.9722 };

// GUItool: begin automatically generated code
AudioInputAnalog         adc1;           //xy=68,274
AudioFilterBiquad        biquad1;        //xy=305,273
AudioOutputUSB           usb1;           //xy=559,160
AudioOutputAnalog        dac1;           //xy=565,313
AudioConnection          patchCord1(adc1, biquad1);
AudioConnection          patchCord2(biquad1, 0, usb1, 0);
AudioConnection          patchCord3(biquad1, 0, usb1, 1);
AudioConnection          patchCord4(biquad1, dac1);
// GUItool: end automatically generated code



void setup() {
  
 AudioMemory(12);
   
  
 float fc = (2000.0+350.0)/2.0;
 float df = (2000.0-350.0);

  for(int ii=0; ii<4;ii++)
  { biquad1.setBandpass(ii,fc,fc/df);
  }

 biquad1.setCoefficients(0,BP_coeff_stag0); 
 biquad1.setCoefficients(1,BP_coeff_stag1); 
 biquad1.setCoefficients(2,BP_coeff_stag2);
 biquad1.setCoefficients(3,BP_coeff_stag3);

  
}

void loop() {

}
 
The following code seems to work
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];

double gain = 1.4291e-04; // corrected 
double Coeff[4][5]=
{{1.0,  2.0, 1.0, -1.67899763, 0.71743918},
 {1.0,  2.0, 1.0, -1.78916723, 0.86019326},
 {1.0, -2.0, 1.0, -1.89588980, 0.90018204},
 {1.0, -2.0, 1.0, -1.96965063, 0.97224260}};
 
void setup() {
  // Audio requires memory to work.
  AudioMemory(12);

  AudioNoInterrupts();
  for(int ii=0;ii<2;ii++)
    for(int jj=0; jj<3;jj++) Coeff[ii][jj] *= sqrt(gain);
    
  for(int ii=0; ii<4; ii++)
  { biquad1.setCoefficients(ii,Coeff[ii]);
  }
  noise1.amplitude(0.7);

  AudioInterrupts();
}


void loop() {
  static uint32_t navg=0;
  if (fft1024.available()) {
    navg++;
    // 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
  }

  static uint32_t t0=0;
  if((millis()-t0)>1000)
  {
    for (int ii=0; ii<16; ii++) {
      Serial.printf("%4d",(int) (1e+4*level[ii]/navg));
    }
    Serial.println();
    for (int ii=0; ii<16; ii++) level[ii]=0;
    navg=0;
    t0=millis();
  }
}

note: averaging is in dB scale and therefore does not reflect mean intensity

For the record
the coefficients are estimated with matlab via
Code: 
f1=350;
f2=2000;
fsh=44100/2;

[z3,p3,k3]=butter(4,[f1,f2]/fsh,'bandpass');
[s3,g3]=zp2sos(z3,p3,k3);

s3x=s3(:,[1,2,3,5,6]);

fprintf('{%.8f, %.8f, %.8f, %.8f, %.8f},\n',s3x')
s3x are the coefficients to be copied to program
g3 is the gain
 
Last edited:
From where this value obtained from =>
double gain = 2.9137e-04;

Also whats exactly doing by the below part of the code?

for(int ii=0;ii<2;ii++)
for(int jj=0; jj<3;jj++) Coeff[ii][jj] *= sqrt(gain);

When using biquad1.setCoefficients(), No need to use biquad1.setBandpass() right??
 
Wicky21 said:
From where this value obtained from =>
double gain = 2.9137e-04;

Also whats exactly doing by the below part of the code?

for(int ii=0;ii<2;ii++)
for(int jj=0; jj<3;jj++) Coeff[ii][jj] *= sqrt(gain);

When using biquad1.setCoefficients(), No need to use biquad1.setBandpass() right??
Click to expand...
as explained in post g3 is the gain of the filter
this needs to be applied to filter
you can apply it full to first set of equations, but I tried it to apply it to the first two biquads (needs sqrt for this)

You either use setBandpass OR setCoefficients
 
Why exactly you selected 2.9137e-04 for the gain? In which range we can vary this gain?( Just asking because then i can vary this an obtain different results)

So if we'r using setBandpass() then it automatically adjust the frequency response of the filter right?We cannot select Butterworth, Chebyshev or other types of response.
 
Wicky21 said:
Why exactly you selected 2.9137e-04 for the gain? In which range we can vary this gain?( Just asking because then i can vary this an obtain different results)

So if we'r using setBandpass() then it automatically adjust the frequency response of the filter right?We cannot select Butterworth, Chebyshev or other types of response.
Click to expand...

the gain is result of matlab.
the selection of filter type is done in matlab.
that is why you use matlab, right?

If you use the setBandpass then you are stuck with the method Paul implemented, which comes from the text file quoted in the header file.
 
Alright. Now i understand well about biquad.
But when comes to fft part. My requirement is to analyze spectrum only between 350Hz - 2000Hz. If so how can i do that ?
 
Wicky21 said:
But when comes to fft part. My requirement is to analyze spectrum only between 350Hz - 2000Hz. If so how can i do that ?
Click to expand...

you cannot do that.
The FFT goes always from 0 to fsh (half sampling frequency)
Better if goes from -fsh to +fsh, but typically one ignores negative frequencies.

to be clear, for a 1024 point fft
-frequency bin 0 is DC
-frequency bins 1 to 511 are positive frequencies
-frequency bins 512 to 1023 are negative frequencies

the only thing you could do is
after biquad filter, decimate by 10 (take every 10th sample) so that your new sampling rate is 4.41 kHz and apply the fft on that.
difference is now the spectral resolution, with 44.1 kHz you have 44100/1023 Hz and with 4.41 kHz you get 4410/1023 Hz resolution

HOWEVER; I do not know if audio library allows that easily. (I'm not using Audio library for my projects)

For your application you may not need that fine resolution, so 44.1 kHz sampling should be fine
How to use fft, please read the different examples and the code.
 
WMXZ said:
Correction, re-doing matlab code I get for gain 1.4291e-04
sorry about that, I will correct in code
Click to expand...

Yes i figured it out. Also your multiplying gain from first 2 biquad, first 3 elemets only.
Code: 
for(int ii=0;ii<2;ii++)
    for(int jj=0; jj<3;jj++) Coeff[ii][jj] *= sqrt(gain);

It should be like

Code: 
for(int ii=0;ii<2;ii++)
    for(int jj=0; jj<5;jj++) Coeff[ii][jj] *= sqrt(gain);

Also for all four biquads

Code: 
for(int ii=0;ii<4;ii++)
    for(int jj=0; jj<5;jj++) Coeff[ii][jj] *= sqrt(gain);
 
Wicky21 said:
Am i right?
Click to expand...

No, you are not right with the other suggestions:

As you can read from biquad docu, the different stages are effectively multiplied together

so, if you wanted to distribute gain over all 4 stages, you have to take sqrt(sqrt(gain)).
I tied it, but results are not good, so I applied it only to the first two stages.

Most likely reason of failure is 16 bit arithmetic in biquad implementation (under- or over-flow).
But,applying gain to first, or sqrt(gain) to first two sections, provided similar data.

All this at change for real very low signal. I have no experience with that. You have to try.

Also, gain is only applied to b0,b1,b2, the numerator, and not the denominator. doing so would cancel the gain effect away.
 
Seems like only below configuration properly working as i desired
Code: 
for(int ii=0;ii<2;ii++)
for(int jj=0; jj<3;jj++) Coeff[ii][jj] *= sqrt(gain);

But in above case when sound levels get lower,frequencies get change slightly (Simply fundamental frquency- 400Hz disappear.At loud sounds, frequencies appears perfectly as Fundamental - 400Hz, 1st Har - 811Hz, 2nd Har - 1211, 3rd Har - 1611Hz)
 
Some additional pieces of advice: (I wrote audio software 20 years ago)

1) A bi-quad "band pass" filter is like an inverted "notch" filter. It falls off on both sides of the center. You don't get to select two different frequencies. You can change the "width" of the "peak" by changing Q, though.

2) A bi-quad implemented in 16/32-bit integer is not very stable at low frequencies. 350 Hz counts as "low" for this case when sampling rate is 44 kHz. (If you can halve the sampling rate, the relative frequencies shift accordingly.)

3) A bi-quad gives you 12 dB/octave filter slope outside of the "Q" affected area. If you want 24 dB/octave, you need to chain two of them.

4) To get selectable bottom and top frequencies for a wide band-pass, you need to chain two separate bi-quads; one low-pass, and one high-pass.

5) A state variable filter, or a simple "leaky integrator" filter, is stable down to whatever frequency you want. It only gives you 6 dB/octave response, though (single pole,) but may be a better option for the low-cut filter.

6) You can gang a number of single-pole filters, but you don't get control of Q.

7) For a specific bandpass filter with specifically chosen frequencies and strong (steep) response outside the pass band, you may be better off with a long FIR filter, plus a single-pole high-pass filter. For FIR to affect frequencies of H/Fs, the rule of thumb is to make the FIR length at least 2*Fs/H (and more is better.)

8) Long FIRs are much more efficiently implemented using FFT convolution than "brute force" convolution.

So, for the requirement of passband between 350 Hz and 2000 Hz, I would probably a convolution length of at least 256 samples (512 would be better) which leads to an overlap-add FFT window size of 512 samples (or 1024) and then couple with one or a few single-pole filter rolling off the low end below 350 Hz. You can make the top stopband very steep this way, and make the bottom end stopband okay. I don't think the 3.2 can keep up with 512 sample FFTs in real time (lack of FPU,) but I think the 3.6 should be able to do that just fine.
 
jwatte said:
Some additional pieces of advice: (I wrote audio software 20 years ago)


4) To get selectable bottom and top frequencies for a wide band-pass, you need to chain two separate bi-quads; one low-pass, and one high-pass.
Click to expand...

thanks jwatte, greatly appreciated.

in fact, considering the coefficients that seem to work, seems indeed to be LP followed by HP (sign of b1)

The method used in audio library biquad function is indeed a single narrowband bandpass filter.

I would agree with FIR, especially FFT-based as number of multiplications become comparably.
2*10*1024 for fft filter
4*5*1024 for 4 stage biquad
(OK, not really correct but only to give an idea)
 
If all you want to do is FFT analysis, and you don't actually need the processed waveform for encoding/playback/whatever, then things become simpler.
All you need to do is to run the FFT on the full spectrum signal at 44 kHz, and simply ignore the buckets (bins) that are outside the area you care about.
You can do this right now, without any need for decimation or anything else.
The resolution per bin will be about 43 Hz. ((44100/2)/(1024/2))
If you need better resolution, then you can indeed decimate. However, decimating by throwing away samples leaves you open to aliasing. You have to first low-pass filter, ideally with a brick wall. You will also still get a bunch of bins that represent frequencies < 350 Hz, that you have to discard.
A FFT is frequency independent, so if you have some FFT code that takes X points in and returns bins, it doesn't matter what the sampling frequency is; the bin width will be relative the sampling frequency, but the FFT itself just cares about number of samples.

You may be able to get to a reasonable trade-off between aliasing and quality and speed and resolution by doing something like, first filter with a 2 Khz corner frequency, using multiple Bi-quads, and then downsampling (decimating) to 11025 Hz sampling frequency, and then running 1024 bin FFT. This should give you a resolution of about 11 Hz per bin. Throw away the first 32 bins or so because they represent the frequencies from DC to 350 Hz. Then keep the next 150 bins or so (range from 350 to 2000 Hz) and then throw away the rest.

Another way to do band-pass is to do FFT convolution. So, use FFT convolution to brick-wall filter everything above 2 kHz. Then re-sample to a lower sample frequency (maybe even 4410 Hz.) Then run the FFT on that, with about 4.3 Hz bin sizes. This runs the stages of FFT -> multiplication -> IFFT -> decimation -> FFT, although the final FFT runs at a lower frequency than the earlier stages.

Why do you need to do this on a Teensy and not on a regular PC? What's your end result application?
 
Main reason is i'm trying to make some simple handy device which can do all these things. Easy to carry and user friendly device. Since teensy 3.2 is really small and power full, thought to use this.
 
jwatte said:
1) You may be able to get to a reasonable trade-off between aliasing and quality and speed and resolution by doing something like, first filter with a 2 Khz corner frequency, using multiple Bi-quads, and then downsampling (decimating) to 11025 Hz sampling frequency, and then running 1024 bin FFT. This should give you a resolution of about 11 Hz per bin. Throw away the first 32 bins or so because they represent the frequencies from DC to 350 Hz. Then keep the next 150 bins or so (range from 350 to 2000 Hz) and then throw away the rest.
Click to expand...

Seems like this method is a better one. But what if we can use some FIR at 2KHz corner frequency instead of multiple bi-quads. Because biquad seems like unstable. Also to obtain much better resolution can't we decimate 44100 by 10 which get 4410Hz with about 4.3Hz bin size. Don't know whether is this possible?

jwatte said:
2) Another way to do band-pass is to do FFT convolution. So, use FFT convolution to brick-wall filter everything above 2 kHz. Then re-sample to a lower sample frequency (maybe even 4410 Hz.) Then run the FFT on that, with about 4.3 Hz bin sizes. This runs the stages of FFT -> multiplication -> IFFT -> decimation -> FFT, although the final FFT runs at a lower frequency than the earlier stages.
Click to expand...

"Use FFT convolution to brick-wall filter everything above 2kHz??" How we can do this?

Simply my requirement is at any background environment just pick relevant sound frequencies (350Hz-2000Hz) without aliasing and observe the FFT. "jwatte" thanks for the bunch valuable information. If i start doing 2) method it will be fine right? Anymore suggestions?
 
Just asking for the additional knowledge.
When making coefficients for bi-quads earlier i have used Butterworth response which has a flat pass band as well as not much steeper with the cutoff edge. So according to my application if i have use some elliptic response but has ripples in both pass band and stop band. But has the steepest cutoff. So by using one stage of biquad can make steeper edge right? When there's ripples in the pass band is that going to attenuate the signal alot? What sort of response is ideal in my case?
 
Because biquad seems like unstable
Click to expand...

At 2 kHz, the bi-quads should be stable, unless you use a crazy Q value.

Decimating 10:1 can only be done AFTER you brick-wall filter. Else you will end up with aliasing mirror images of frequencies in the output signal, which will give you false readings.

"Use FFT convolution to brick-wall filter everything above 2kHz??" How we can do this?
Click to expand...

I haven't looked at the Teensy code for FFTs in detail. The textbook example is the "overlap-add" technique, where you determine a FIR response of block size (say, 512 samples,) and then FFT your FIR from a 1024 block where the second half is all zeros, once. Call this the "impulse response." Then, you input each successive 512-sample block into a 1024-sample buffer, FFT it, multiply the impulse response with the FFT-ed data, and IFFT back, which gives you 1024 samples out. Output 512 samples, and save 512 samples that you add to the next batch. (This is why it's "overlap/add")

https://en.wikipedia.org/wiki/Overlap–add_method

So, yes, if you brick-wall filter out anything above 2 Khz, using overlap-add FFT convolution, then you can decimate 10:1, after which you can analyze using FFT, and then throw away (ignore) all bins for frequencies below 350 Hz.

Also note that, depending on implementation, the FFT data bins may NOT be in numerical order (i e, index 3 is not necessarily bin 3) but instead may be in "butterfly" order. For the FFT - multiply - IFFT step, this is actually preferred, because it's faster, but for the final analysis/visualization, you want to de-butterfly the bin ordering to get a 1:1 relation between bin index and frequency.

It seems like the code for the Teensy only does forward FFT (in the AudioAnalyzeFFT class) so you'll want to find some other FFT library that does both FFT, multiplication, and IFFT for you. Back in the day, I found "djbfft" to be fast and accurate, but it's mainly aimed at x86, so maybe there's something else that does better on ARM with single precision float.
 
Status
Not open for further replies.
Back
Top