Digital BIQUAD filter calculations

TC · 2016-04-12 11:52

Hello all,

I am building my own bluetooth headphones. I am currently working on the DSP for a 5 band EQ. I was going to make a database for the biquad coefficients. Then just look up the coefficient value. But I then realized the amount of data that would be. I am going with the tlv320aic3268 from Texas Instruments. Each coefficient is 24bits wide, and each BIQUAD need 5 coefficients. That would be 5 longs for each band of the EQ (assuming left and right use same values).

5 longs * 5 bands = 25 longs

That's not a big deal, but I want to have the option to change the gain of each EQ band from +12db to -12db, in 0.1db steps. That is about 240(ish) steps.

5 coefficients (24bits each) * 5 bands * 240 steps = 6000 longs

That's a lot of data. So I thought, maybe the Prop can do the coefficient calculations on the fly. I just put in the EQ frequency, the "Q", and the gain, and the Prop could perform the calculation, and load the coefficients. I found papers that explains the biquad calculations (EQ cookbook), and one the gives Matlab equations (Texas). But when I do the equations in excel, I get incorrect results. I am using Texas's coefficients calculator to verify my results.

I am wondering if anyone would have some advice, or if anyone would be able to help me make the prop do coefficient calculations on the fly.

Thanks
TC

TC · 2016-04-13 12:23

Hello all,

I made a LPF biquad filter calculator code for the prop. This is just to see if I can get the same values as the Texas Instruments biquad calculator. The values I am getting are a little off from what the texas calculator is giving me. Could someone please take a look at it? I might of missed something, or I have something wrong.

user define values:
Low Pass Filter
sample freq = 48000 Hz
cutoff freq = 600 Hz
Q = 0.707

Values from Texas calculator

N0 = $002FD9
N1 = $002FD9
N2 = $002FD9
D1 = $78E5AB
D2 = $8D7541

Values from prop

N0 = $002FC0
N1 = $002FC0
N2 = $002FC0
D1 = $78E5D1
D2 = $8D755D

Thanks
TC

CON
        _clkmode = xtal1 + pll16x                                               'Standard clock mode * crystal frequency = 80 MHz
        _xinfreq = 5_000_000

CON

  baud          = 250000
  sample        = 48000
  range         = 8388608.0  '2^(number of bits - 1)-1
   
OBJ

  Math          : "Float32Full"
  pst           : "Parallax Serial Terminal"
  FtoS          : "FloatString"

pub start

  math.start
  pst.start(baud)

  LPF(600, 0.707, 0)


PUB LPF(freq, Q, gain)| w, alpha, cosW, sinW, b0, b1, b2, a0, a1, a2, N0_real, N1_real, N2_real, N0, N1, N2, D1, D2, factor, N0_hex, N1_hex, N2_hex, D1_hex, D2_hex, hex_temp1, hex_temp2, hex_temp3

  ' w = 2*PI*f0/Fs
  w := math.Fmul(math.Fmul(2.0, pi), math.Fdiv(freq, sample))
  
      pst.str(string("w = "))
      pst.str(FtoS.FloatToString(w))
      pst.newline  

  cosW := Math.Cos(w)

      pst.str(string("cosW = "))
      pst.str(FtoS.FloatToString(cosW))
      pst.newline
  
  sinW := Math.Sin(w)
  
      pst.str(string("sinW = "))
      pst.str(FtoS.FloatToString(sinW))
      pst.newline

  'alpha = sin(w)/(2*Q)
  alpha := math.Fdiv(sinW, math.Fmul(2.0, Q))
  
      pst.str(string("alpha = "))
      pst.str(FtoS.FloatToString(alpha))
      pst.newline

  'b0 = (1-cos(w))/2
  b0 := math.Fdiv(math.Fsub(1.0, cosW), 2.0)
  
      pst.str(string("b0 = "))
      pst.str(FtoS.FloatToString(b0))
      pst.newline

  'b1 = 1-cos(w)
  b1 := math.FSub(1.0, cosW)
  
      pst.str(string("b1 = "))
      pst.str(FtoS.FloatToString(b1))
      pst.newline

  'b2 = (1-cos(w))/2
  b2 := math.Fdiv(math.Fsub(1.0, cosW), 2.0)
  
      pst.str(string("b2 = "))
      pst.str(FtoS.FloatToString(b2))
      pst.newline

  'a0 = 1+alpha
  a0 := math.Fadd(1.0, alpha)
  
      pst.str(string("a0 = "))
      pst.str(FtoS.FloatToString(a0))
      pst.newline

  'a1 = -2*cos(w)
  a1 := math.Fmul(-2.0, cosW)
  
      pst.str(string("a1 = "))
      pst.str(FtoS.FloatToString(a1))
      pst.newline

  'a2 = 1-alpha
  a2 := math.Fsub(1.0, alpha)
  
      pst.str(string("a2 = "))
      pst.str(FtoS.FloatToString(a2))
      pst.newline

  'N0_real = (b0/a0)
  N0_real := math.Fdiv(b0, a0)
  
      pst.str(string("N0_real = "))
      pst.str(FtoS.FloatToString(N0_real))
      pst.newline

  'N1_real = b1/a0/2
  N1_real := math.Fdiv(math.Fdiv(b1, a0), 2.0)
  
      pst.str(string("N1_real = "))
      pst.str(FtoS.FloatToString(N1_real))
      pst.newline

  'N2_real = (b2/a0)
  N2_real := math.Fdiv(b2, a0)
  
      pst.str(string("N2_real = "))
      pst.str(FtoS.FloatToString(N2_real))
      pst.newline

  'scale the coefficients
  factor := math.Fmax(math.Fmax(N1_real, N2_real), N0_real)
  
    if math.Fcmp(factor, 1.0) == -1
      factor := 1.0
    
  pst.str(string("Factor = "))
  pst.str(FtoS.FloatToString(factor))
  pst.newline

  'convert to quantized coefficients

  'N0 = ((b0/a0)/factor)*range
  N0 := math.Fmul(math.Fdiv(math.Fdiv(b0, a0), factor), range)
  
      pst.str(string("N0 = "))
      pst.str(FtoS.FloatToString(N0))
      pst.newline

  'N1 = ((b1/a0)/(2*factor))/range
  N1 := math.Fmul(math.Fdiv(math.Fdiv(b1, a0), math.Fmul(2.0, factor)), range)
  
      pst.str(string("N1 = "))
      pst.str(FtoS.FloatToString(N1))
      pst.newline

  'N2 = ((b2/a0)/factor)*range
  N2 := math.Fmul(math.Fdiv(math.Fdiv(b2, a0), factor), range)
  
      pst.str(string("N2 = "))
      pst.str(FtoS.FloatToString(N2))
      pst.newline

  'D1 = (a1/(-2*a0))*range
  D1 := math.Fmul(math.Fdiv(a1, math.Fmul(-2.0, a0)), range)
  
      pst.str(string("D1 = "))
      pst.str(FtoS.FloatToString(D1))
      pst.newline

  'D2 = (a2/a0)*range
  D2 := math.Fmul(math.Fneg(math.Fdiv(a2, a0)), range)
  
      pst.str(string("D2 = "))
      pst.str(FtoS.FloatToString(D2))
      pst.newline

  N0_hex := math.Fround(N0)
  
      pst.str(string("N0_hex = "))
      pst.hex(N0_hex,6)
      pst.newline

  N1_hex := math.Fround(N1)
  
      pst.str(string("N1_hex = "))
      pst.hex(N1_hex,6)
      pst.newline

  N2_hex := math.Fround(N2)
  
      pst.str(string("N2_hex = "))
      pst.hex(N2_hex,6)
      pst.newline

  D1_hex := math.Fround(D1)
  
      pst.str(string("D1_hex = "))
      pst.hex(D1_hex,6)
      pst.newline

  D2_hex := math.Fround(D2)
  
      pst.str(string("D2_hex = "))
      pst.hex(D2_hex,6)
      pst.newline

Peter Jakacki · 2016-04-13 12:30

That may simply be a lack of precision in the Prop Floating point (only 32-bit) but personally I would simply load the 24k table in the top 32k of the EEPROM since all systems these days have a 64k EEPROM anyway.

TC · 2016-04-13 12:54

Thank you Peter,

I was wondering about the precision of the prop. The only main reason i didn't want to do a table of the values is because I am lazy. But I don't think I really have any real options. At least doing a table would make sure the correct values would be loaded.

Peter Jakacki · 2016-04-13 13:12

If you can do it on the Prop then what is stopping you from writing some script maybe even in a console then nabbing the results and then if necessary massage that into a form that you can load into the EEPROM. You might need a simple Spin program to load it for you.

ersmith · 2016-04-13 13:39

This isn't an issue of the precision of the propeller itself, or even of 32 bit math, it's of the particular math library you're using (Float32Full, in this case). My guess is that the trig functions in Float32Full are using the ROM lookup table and interpolating, and are not particularly accurate as a result. Also, many Spin floating point implementations trade accuracy for speed/simplicity. If you really want to do the calculations at run time it'd be worth looking at using a different library and/or rolling your own sin/cos functions.

I tried this in PropGCC and it worked fine; there are some 1 bit differences from the TI calculator, but the PropGCC values agreed with what I got on my PC. On the Prop I tried with both 32 bit and 64 bit doubles, and got the same results in both cases; 32 bit doubles are a lot faster, so that would be the sensible choice. The code is pretty big though (due to printf) so you'll definitely want to use CMM mode.

values from PropGCC
N0 = $002fda
N1 = $002fda
N2 = $002fda
D1 = $78e5ab
D2 = $8d7541

//
// biquad calculations
//
#include <stdio.h>
#include <math.h>

#define myround(x) ((int)((x)+0.5))

double sample = 48000;
double range = 8388608.0;

void LPF(int freq, double Q, int gain)
{
    double w = 2 * M_PI * freq / sample;
    double cosw = cos(w);
    double sinw = sin(w);
    double alpha = sinw / (2.0*Q);
    double b0 = (1.0-cosw) / 2.0;
    double b1 = 1.0 - cosw;
    double b2 = b1 / 2.0;
    double a0 = 1.0 + alpha;
    double a1 = -2.0*cosw;
    double a2 = 1.0-alpha;
    double N0_real = (b0/a0);
    double N1_real = (b1/a0)/2.0;
    double N2_real = (b2/a0);
    double factor;
    double N0, N1, N2;
    double D1, D2;
    int N0_hex, N1_hex, N2_hex;
    int D1_hex, D2_hex;
    
    factor = fmax(fmax(N1_real, N2_real), N0_real);
    if (factor < 1.0)
        factor = 1.0;

    N0 = ((b0/a0)/factor) * range;
    N1 = ((b1/a0)/(2*factor)) * range;
    N2 = ((b2/a0)/factor) * range;
    D1 = (a1/(-2.0*a0))*range;
    D2 = (-(a2/a0))*range;

    N0_hex = myround(N0);
    N1_hex = myround(N1);
    N2_hex = myround(N2);
    D1_hex = myround(D1);
    D2_hex = myround(D2) & 0x00ffffff;

    printf("N0 = $%06x\n", N0_hex);
    printf("N1 = $%06x\n", N1_hex);
    printf("N2 = $%06x\n", N2_hex);
    printf("D1 = $%06x\n", D1_hex);
    printf("D2 = $%06x\n", D2_hex);
    
}

int
main()
{
    LPF(600, 0.707, 0);
    return 0;
}

Peter Jakacki · 2016-04-13 13:44

So my suggestion would be to run ersmith's code with some other code that writes it to EEPROM. Then you can run your code as normal just looking up the values from EEPROM.

TC · 2016-04-14 10:39

ersmith wrote: »

This isn't an issue of the precision of the propeller itself, or even of 32 bit math, it's of the particular math library you're using (Float32Full, in this case). My guess is that the trig functions in Float32Full are using the ROM lookup table and interpolating, and are not particularly accurate as a result. Also, many Spin floating point implementations trade accuracy for speed/simplicity. If you really want to do the calculations at run time it'd be worth looking at using a different library and/or rolling your own sin/cos functions.

I tried this in PropGCC and it worked fine; there are some 1 bit differences from the TI calculator, but the PropGCC values agreed with what I got on my PC. On the Prop I tried with both 32 bit and 64 bit doubles, and got the same results in both cases; 32 bit doubles are a lot faster, so that would be the sensible choice. The code is pretty big though (due to printf) so you'll definitely want to use CMM mode.

Thank you ersmith. I never learned trig in school, and I am learning as I go. It makes me feel better knowing that the problem is not with my code, only the library I am using.

I'm going to do what Peter suggested, use ersmith's code and load a eeprom with the values. That would be a lot easier than me hand typing all the values.

Thank you so much for the help.
TC

Digital BIQUAD filter calculations

Comments