Simple Audio Filter

New to XMOS and XCore? Get started here.
Rammsteiner1988
Member
Posts: 12
Joined: Thu Dec 05, 2013 9:39 am

Post by Rammsteiner1988 »

OK I can here music!!! very nice!!!
But If I make a sound with Delphi (windows.Beep(15000,1000)) that are 15000 Hz I still can here.
So what is the cutoff frequency in your coeffs?


User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

[B,A]=butter(2,0.1); in MATLAB means 2:nd order butterworth at 0.1 * (½FS).
That would be 2205 Hz @ 44.1 kHz.

The next improvement to the filter is to replace

h<<2 with (h<<2) bitwise or (l>>30)
That will give your filter 12 dB better Dynamic range.

the 64 bit accumulator {h,l} should be divided with a0 or shifted >>30 => {h,l}/a0 = {h,l}>>30

h is already shifted <<32 so that will become h<<2 , and when finally add it together with or l>>30.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

You can use this calcFilt() to calculate your own coeffs. Not so "clean" code but anyway:
It works as long as all normalized |coefs| < 2, when larger you need to scale everything and change the shifting in the filter. As long that your filter has a maximum of unity-gain it will work.

Code: Select all

double norm(double a,double b,double c){
	double max;
	max=fabs(a);
	if(fabs(b)>max)
		max=fabs(b);
	if(fabs(c)>max)
		max=fabs(c);
	return max;
}


void calcFilt(unsigned f0, unsigned  mQ, int Gain_mdB , unsigned fs, unsigned Filtertype,int Bcoef[3],int Acoef[2],int &Shift , double masterGain){
	/*Calculate parametric EQ coef
    %f0 is the filter frequency
    %mQ is the Q value in 1/1000, milli
    %Gain_mdB is the gain in dB/1000
    %fs is the sampling frequency
    %Filtertype is one of
    % 1 LPF (Lowpass Filter)
    % 2 HPF (High pass filter)
    % 3 BPFs (Band pass filter constant skirt gain, peak gain = Q)
    % 4 BPFp (constant 0 dB peak gain)
    % 5 notch
    % 6 APF AllPass filter
    % 7 peakingEQ
    % 8 lowShelf
    % 9 highSelf*/

	double SCALE,max;
	const double pi=3.14159265358979323846;
    double w0 = 2 * pi * (double)f0/fs;
    double cos_w0 = cos(w0);
    double alpha=sin(w0)/(2*(double)mQ/1000);
    double A = pow(10,((double)Gain_mdB/40000));
    double sqA;
    double Gain=pow(10,((double)Gain_mdB/20000));
    double b0, b1, b2 ,a0 ,a1 ,a2;
    switch(Filtertype){
    	case LPF:
    		b0 = Gain *  (1 - cos_w0)/2;
            b1 = Gain *  2*b0;
            b2 = Gain *  b0;
            a0 =   1 + alpha;
            a1 =  -2*cos_w0;
            a2 =   1 - alpha;
        break;
    	case HPF:
            b0 = Gain *  (1 + cos_w0)/2;
            b1 = Gain *  -2*b0;
            b2 = Gain *  b0;
            a0 =   1 + alpha;
            a1 =  -2*cos_w0;
            a2 =   1 - alpha;
       break;
    case BPFs :
            b0 =   sin(w0)/2 ;
            b1 =   0 ;
            b2 =  -b0;
            a0 =   1 + alpha;
            a1 =  -2*cos_w0;
            a2 =   1 - alpha;
            break;
    break;
    case BPFp:
            b0 = Gain *  alpha;
            b1 =   0;
            b2 = Gain * -alpha;
            a0 =   1 + alpha;
            a1 =  -2*cos_w0;
            a2 =   1 - alpha;
            break;
     case NOTCH:
            b0 =   1;
            b1 =  -2*cos_w0;
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos_w0;
            a2 =   1 - alpha;
            break;
     case APF:
            b0 =   1 - alpha;
            b1 =  -2*cos_w0;
            b2 =   1 + alpha;
            a0 =   1 + alpha;
            a1 =  -2*cos_w0;
            a2 =   1 - alpha;
            break;
     case peakingEQ:
            b0 =   1 + alpha*A;
            b1 =  -2*cos_w0;
            b2 =   1 - alpha*A;
            a0 =   1 + alpha/A;
            a1 =  -2*cos_w0;
            a2 =   1 - alpha/A;
            break;
      case lowShelf:
            sqA=sqrt(A);
    	    b0 =    A*( (A+1) - (A-1)*cos_w0 + 2*sqA*alpha );
            b1 =  2*A*( (A-1) - (A+1)*cos_w0                   );
            b2 =    A*( (A+1) - (A-1)*cos_w0 - 2*sqA*alpha );
            a0 =        (A+1) + (A-1)*cos_w0 + 2*sqA*alpha;
            a1 =   -2*( (A-1) + (A+1)*cos_w0                   );
            a2 =        (A+1) + (A-1)*cos_w0 - 2*sqA*alpha;
      break;
      case highShelf:
    	    sqA=sqrt(A);
    	  	b0 =    A*( (A+1) + (A-1)*cos_w0 + 2*sqA*alpha );
            b1 = -2*A*( (A-1) + (A+1)*cos_w0                   );
            b2 =    A*( (A+1) + (A-1)*cos_w0 - 2*sqA*alpha );
            a0 =        (A+1) - (A-1)*cos_w0 + 2*sqA*alpha;
            a1 =    2*( (A-1) - (A+1)*cos_w0                   );
            a2 =        (A+1) - (A-1)*cos_w0 - 2*sqA*alpha;
      break;
    }
    SCALE=1<<30;
    Acoef[0]=(int) (SCALE* -a1/a0);
    Acoef[1]=(int) (SCALE* -a2/a0);

max = norm(b0/a0,b1/a0,b2/a0);
if(max<2){
	Shift=2;
	}
else{
	Shift=2+floor(log(max)/log(2));
	SCALE=1<<(32-Shift);
	printstrln("Shift !=2");
}
Bcoef[0]=(int) (masterGain*SCALE* b0/a0);
Bcoef[1]=(int) (masterGain*SCALE* b1/a0);
Bcoef[2]=(int) (masterGain*SCALE* b2/a0);

}
include print.h
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

For "PC" or X86, here is a free Delphi program http://tolvan.com/index.php?page=/rtsect/rtsect.php that you can use to check your transfer function "live"
Rammsteiner1988
Member
Posts: 12
Joined: Thu Dec 05, 2013 9:39 am

Post by Rammsteiner1988 »

THANK YOU!!! The filter is working!
One last question why scaling with 2^30 and with 2^8 or 2^16 ???
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

Take a look at the normalized coefs:

B={0.020083365564211 , 0.040166731128423 ,0.020083365564211}
A = { 1.000000000000000 , -1.561018075800718 , 0.641351538057563}

if we would use 2^8=256 and rounded the values to integer we would get


B = {5 , 10 , 5}
A = {256 , -400 , 164}

If we now normalize this again we would get.
B={0.01953125 , 0.0390625 , 0.01953125}
A={ 1 , -1.5625 , 0.640625}

This is in fact another filter compared to the one we started with. If the first filter have a large Q-value, the quantization error can make the filter unstable.

The rounding to integer introduces quantization error seen above, and the way to minimize this error is to use to largest possible scale value.
rmammina
Member
Posts: 9
Joined: Thu Mar 20, 2014 7:24 pm

Post by rmammina »

So are you suggesting that we use Q15 instead of floats or doubles with our biquad filter coefficients?
bearcat
Respected Member
Posts: 283
Joined: Fri Mar 19, 2010 4:49 am

Post by bearcat »

Thanks Lilltroll for the filter toolkit XC function. Been needing to write that for a while and had not seen it yet.

Would you happen to have a version using BW instead of Q? Save me even more time.
Post Reply