Jedora

Happy blogging in my way...

Friday, June 08, 2007

FIR filter - from Matlab to C

How to implement Matlab fir1 function in C? Here is the answer.

//================== my_fir1.c ============================
#include
#include
#include

#define PI 3.14159265
#define L 256 /* Filter length */
#define Pr_L L
#define M 64

double * my_fir1(unsigned N, double Wn)
{
/* FIR filter
Return an array of 1 X 256
*/
unsigned odd, i, j, nhlf, i1;
double f1, gain, c1;
double wind[Pr_L],xn[Pr_L/2],b[Pr_L/2],c[Pr_L/2],c3[Pr_L/2];
double *bb;

bb = (double *) malloc(sizeof(double) * Pr_L);

gain = 0.0000000;
N = N+1;
odd = N - (N/2)*2; /* odd = rem(N,2) */

/*wind = hamming(N);*/
for (i=0; i < Pr_L; i++)
{
wind[i] = 0.54 - 0.46 * cos ((2 * PI * i) / (N-1));
}

f1 = Wn / 2.0;
c1 = f1;
nhlf = (N+1) / 2;
i1 = odd + 1;

/* Lowpass */

if(odd)
b[0] = 2 * c1;

for (i=0; i < nhlf; i++)
{
xn[i] = i + 0.5 * (1 - odd);
}

for (i=0; i < nhlf; i++)
{
c[i] = PI * xn[i];
}

for (i=0; i < nhlf; i++)
{
c3[i] = 2 * c1 * c[i];
}

/* b(i1:nhlf)=(sin(c3)./c) */
for (i=0; i < nhlf; i++)
{
b[i] = sin(c3[i]) / c[i];
}

/* bb = real([b(nhlf:-1:i1) b(1:nhlf)].*wind(:)') */
for (i=0,j=nhlf-1; i < nhlf; i++, j--)
{
bb[i] = b[j];
}
for (i=nhlf,j=0; i < Pr_L; i++,j++)
{
bb[i] = b[j];
}
for (i=0; i < Pr_L; i++)
{
bb[i] = bb[i] * wind[i];
}

/* gain = abs(polyval(b,1)); */
for (i=0; i < Pr_L; i++)
{
gain += bb[i];
}
/* b = b / gain */
for (i=0; i < Pr_L; i++)
{
bb[i] = bb[i] / gain;
}

return bb;

}

int main (void)
{
unsigned i;
double *Prot_filt;
Prot_filt = (double *) malloc(sizeof(double) * Pr_L);

Prot_filt = my_fir1(Pr_L-1,1 / (double)M); /* FIR filter */

for(i=0; i < Pr_L; i++)
{
Prot_filt[i] = Prot_filt[i] * (M/2);
}

for(i=0; i < Pr_L; i++)
printf("Prot_filt[%d]=%f\n",i,Prot_filt[i]); /*debug*/

}

======================
I test it under Fedora Core 6 Linux

$ gcc -o my_fir1 my_fir1.c
$ ./my_fir1


Outputs are:
Prot_filt[0]=-0.000156
Prot_filt[1]=-0.000473
Prot_filt[2]=-0.000797
Prot_filt[3]=-0.001132
Prot_filt[4]=-0.001480
Prot_filt[5]=-0.001844
Prot_filt[6]=-0.002227
Prot_filt[7]=-0.002631
Prot_filt[8]=-0.003059
Prot_filt[9]=-0.003513
Prot_filt[10]=-0.003996
Prot_filt[11]=-0.004510
Prot_filt[12]=-0.005056
Prot_filt[13]=-0.005637
Prot_filt[14]=-0.006254
Prot_filt[15]=-0.006907
Prot_filt[16]=-0.007598
Prot_filt[17]=-0.008327
Prot_filt[18]=-0.009095
Prot_filt[19]=-0.009901
Prot_filt[20]=-0.010744
Prot_filt[21]=-0.011623
Prot_filt[22]=-0.012537
Prot_filt[23]=-0.013484
Prot_filt[24]=-0.014462
Prot_filt[25]=-0.015466
Prot_filt[26]=-0.016495
Prot_filt[27]=-0.017543
Prot_filt[28]=-0.018608
Prot_filt[29]=-0.019682
Prot_filt[30]=-0.020762
Prot_filt[31]=-0.021841
Prot_filt[32]=-0.022913
Prot_filt[33]=-0.023971
Prot_filt[34]=-0.025008
Prot_filt[35]=-0.026015
Prot_filt[36]=-0.026985
Prot_filt[37]=-0.027908
Prot_filt[38]=-0.028777
Prot_filt[39]=-0.029580
Prot_filt[40]=-0.030309
Prot_filt[41]=-0.030954
Prot_filt[42]=-0.031504
Prot_filt[43]=-0.031948
Prot_filt[44]=-0.032275
Prot_filt[45]=-0.032476
Prot_filt[46]=-0.032538
Prot_filt[47]=-0.032451
Prot_filt[48]=-0.032203
Prot_filt[49]=-0.031783
Prot_filt[50]=-0.031181
Prot_filt[51]=-0.030384
Prot_filt[52]=-0.029382
Prot_filt[53]=-0.028165
Prot_filt[54]=-0.026723
Prot_filt[55]=-0.025044
Prot_filt[56]=-0.023120
Prot_filt[57]=-0.020940
Prot_filt[58]=-0.018497
Prot_filt[59]=-0.015782
Prot_filt[60]=-0.012786
Prot_filt[61]=-0.009504
Prot_filt[62]=-0.005928
Prot_filt[63]=-0.002052
Prot_filt[64]=0.002129
Prot_filt[65]=0.006618
Prot_filt[66]=0.011420
Prot_filt[67]=0.016538
Prot_filt[68]=0.021972
Prot_filt[69]=0.027725
Prot_filt[70]=0.033795
Prot_filt[71]=0.040183
Prot_filt[72]=0.046886
Prot_filt[73]=0.053900
Prot_filt[74]=0.061223
Prot_filt[75]=0.068849
Prot_filt[76]=0.076772
Prot_filt[77]=0.084985
Prot_filt[78]=0.093480
Prot_filt[79]=0.102248
Prot_filt[80]=0.111278
Prot_filt[81]=0.120560
Prot_filt[82]=0.130081
Prot_filt[83]=0.139829
Prot_filt[84]=0.149788
Prot_filt[85]=0.159945
Prot_filt[86]=0.170283
Prot_filt[87]=0.180786
Prot_filt[88]=0.191435
Prot_filt[89]=0.202213
Prot_filt[90]=0.213101
Prot_filt[91]=0.224079
Prot_filt[92]=0.235126
Prot_filt[93]=0.246222
Prot_filt[94]=0.257346
Prot_filt[95]=0.268474
Prot_filt[96]=0.279586
Prot_filt[97]=0.290658
Prot_filt[98]=0.301668
Prot_filt[99]=0.312592
Prot_filt[100]=0.323407
Prot_filt[101]=0.334091
Prot_filt[102]=0.344619
Prot_filt[103]=0.354969
Prot_filt[104]=0.365118
Prot_filt[105]=0.375043
Prot_filt[106]=0.384721
Prot_filt[107]=0.394131
Prot_filt[108]=0.403250
Prot_filt[109]=0.412058
Prot_filt[110]=0.420534
Prot_filt[111]=0.428659
Prot_filt[112]=0.436412
Prot_filt[113]=0.443777
Prot_filt[114]=0.450734
Prot_filt[115]=0.457267
Prot_filt[116]=0.463361
Prot_filt[117]=0.469001
Prot_filt[118]=0.474174
Prot_filt[119]=0.478866
Prot_filt[120]=0.483066
Prot_filt[121]=0.486764
Prot_filt[122]=0.489950
Prot_filt[123]=0.492618
Prot_filt[124]=0.494760
Prot_filt[125]=0.496371
Prot_filt[126]=0.497448
Prot_filt[127]=0.497987
Prot_filt[128]=0.497987
Prot_filt[129]=0.497448
Prot_filt[130]=0.496371
Prot_filt[131]=0.494760
Prot_filt[132]=0.492618
Prot_filt[133]=0.489950
Prot_filt[134]=0.486764
Prot_filt[135]=0.483066
Prot_filt[136]=0.478866
Prot_filt[137]=0.474174
Prot_filt[138]=0.469001
Prot_filt[139]=0.463361
Prot_filt[140]=0.457267
Prot_filt[141]=0.450734
Prot_filt[142]=0.443777
Prot_filt[143]=0.436412
Prot_filt[144]=0.428659
Prot_filt[145]=0.420534
Prot_filt[146]=0.412058
Prot_filt[147]=0.403250
Prot_filt[148]=0.394131
Prot_filt[149]=0.384721
Prot_filt[150]=0.375043
Prot_filt[151]=0.365118
Prot_filt[152]=0.354969
Prot_filt[153]=0.344619
Prot_filt[154]=0.334091
Prot_filt[155]=0.323407
Prot_filt[156]=0.312592
Prot_filt[157]=0.301668
Prot_filt[158]=0.290658
Prot_filt[159]=0.279586
Prot_filt[160]=0.268474
Prot_filt[161]=0.257346
Prot_filt[162]=0.246222
Prot_filt[163]=0.235126
Prot_filt[164]=0.224079
Prot_filt[165]=0.213101
Prot_filt[166]=0.202213
Prot_filt[167]=0.191435
Prot_filt[168]=0.180786
Prot_filt[169]=0.170283
Prot_filt[170]=0.159945
Prot_filt[171]=0.149788
Prot_filt[172]=0.139829
Prot_filt[173]=0.130081
Prot_filt[174]=0.120560
Prot_filt[175]=0.111278
Prot_filt[176]=0.102248
Prot_filt[177]=0.093480
Prot_filt[178]=0.084985
Prot_filt[179]=0.076772
Prot_filt[180]=0.068849
Prot_filt[181]=0.061223
Prot_filt[182]=0.053900
Prot_filt[183]=0.046886
Prot_filt[184]=0.040183
Prot_filt[185]=0.033795
Prot_filt[186]=0.027725
Prot_filt[187]=0.021972
Prot_filt[188]=0.016538
Prot_filt[189]=0.011420
Prot_filt[190]=0.006618
Prot_filt[191]=0.002129
Prot_filt[192]=-0.002052
Prot_filt[193]=-0.005928
Prot_filt[194]=-0.009504
Prot_filt[195]=-0.012786
Prot_filt[196]=-0.015782
Prot_filt[197]=-0.018497
Prot_filt[198]=-0.020940
Prot_filt[199]=-0.023120
Prot_filt[200]=-0.025044
Prot_filt[201]=-0.026723
Prot_filt[202]=-0.028165
Prot_filt[203]=-0.029382
Prot_filt[204]=-0.030384
Prot_filt[205]=-0.031181
Prot_filt[206]=-0.031783
Prot_filt[207]=-0.032203
Prot_filt[208]=-0.032451
Prot_filt[209]=-0.032538
Prot_filt[210]=-0.032476
Prot_filt[211]=-0.032275
Prot_filt[212]=-0.031948
Prot_filt[213]=-0.031504
Prot_filt[214]=-0.030954
Prot_filt[215]=-0.030309
Prot_filt[216]=-0.029580
Prot_filt[217]=-0.028777
Prot_filt[218]=-0.027908
Prot_filt[219]=-0.026985
Prot_filt[220]=-0.026015
Prot_filt[221]=-0.025008
Prot_filt[222]=-0.023971
Prot_filt[223]=-0.022913
Prot_filt[224]=-0.021841
Prot_filt[225]=-0.020762
Prot_filt[226]=-0.019682
Prot_filt[227]=-0.018608
Prot_filt[228]=-0.017543
Prot_filt[229]=-0.016495
Prot_filt[230]=-0.015466
Prot_filt[231]=-0.014462
Prot_filt[232]=-0.013484
Prot_filt[233]=-0.012537
Prot_filt[234]=-0.011623
Prot_filt[235]=-0.010744
Prot_filt[236]=-0.009901
Prot_filt[237]=-0.009095
Prot_filt[238]=-0.008327
Prot_filt[239]=-0.007598
Prot_filt[240]=-0.006907
Prot_filt[241]=-0.006254
Prot_filt[242]=-0.005637
Prot_filt[243]=-0.005056
Prot_filt[244]=-0.004510
Prot_filt[245]=-0.003996
Prot_filt[246]=-0.003513
Prot_filt[247]=-0.003059
Prot_filt[248]=-0.002631
Prot_filt[249]=-0.002227
Prot_filt[250]=-0.001844
Prot_filt[251]=-0.001480
Prot_filt[252]=-0.001132
Prot_filt[253]=-0.000797
Prot_filt[254]=-0.000473
Prot_filt[255]=-0.000156

==========================
It's equivalent to Matlab outputs.
Hope useful for someone who's looking for this function.