IIR Filter using Bilinear Transformation

 IIR Filter  using Bilinear Transformation

Problem Statement :Design of IIR filter for given specifications using Bilinear Transformation. (Programshouldwork for different types of filter specifications i.e LPF, HPF, BPF etc and for different
transfer functions of an analog filter).

1:                 // biliner.cpp  
2:  #include<stdio.h>  
3:  #include<conio.h>  
4:  #include<math.h>  
5:  void main()  
6:  {  
7:     float B[20],b[20][3],a[20][3],pReal,pImag;  
8:     float wc,N,pi,Theta,omegaC,den;  
9:     int k;  
10:     clrscr();  
11:     printf("\t\tButterworth filter design using Bilinear "  
12:                                 "transformation\n\n");  
13:     printf("Enter the order of the filter N = ");  
14:     scanf("%f",&N);               // order of the filter  
15:     printf("\nEnter the cutoff frequency of digital filter wc = ");  
16:     scanf("%f",&wc);   // cutoff frequency of digital filter i.e. wc  
17:     omegaC = tan(wc/2); // cutoff frequency of equivalent analog  
18:  //              filter by bilinear frequency relationship  
19:     pi = 22.0/7.0;                   // value of pi  
20:     for(k = 0; k < N/2; k++)// loop for computation of first N/2 poles  
21:     {  
22:      Theta = ((N+2*k+1)*pi)/(2*N);       // angle of k'th pole  
23:      pReal = -1*omegaC*cos(Theta); // real part of pole made positive  
24:      pImag = omegaC*sin(Theta);       // imaginary part of pole  
25:        den = 1+2*pReal+pReal*pReal+pImag*pImag;  
26:       B[k] = (omegaC*omegaC)/den;        // calculation of Bk  
27:      b[k][0] = 1;  //|  for all the sections the value of bk0 = 1,  
28:      b[k][1] = 2;  //|      bk1 = 2 & bk2 = 1. These values are  
29:      b[k][2] = 1;  //|        fixed(see theory for details).  
30:      a[k][0] = 1;             // value of ak0 is always 1  
31:      a[k][1] = (2*(pReal*pReal+pImag*pImag)-2)/den;      // ak1  
32:      a[k][2] = (1-2*pReal+pReal*pReal+pImag*pImag)/den;    // ak2  
33:               // the above two statements calculate ak1 and ak2  
34:     }  
35:     if((N/2) != k)  
36:     {  
37:      k--;       // recompute the pole without complex conjugate  
38:      den = omegaC+1;  
39:      B[k] = omegaC/den;             // calculation of Bk  
40:      b[k][0] = 1; //|  for first order section bk0 = 1, bk1 = 1 and  
41:      b[k][1] = 1; //|     bk2 = 0 always. These values are fixed  
42:      b[k][2] = 0; //|           and need not be calculated  
43:      a[k][0] = 1;             // value of ak0 is always 1  
44:      a[k][1] = (omegaC-1)/den;         // calculation of ak1  
45:      a[k][2] = 0;     // ak2 is always 0 for first order section  
46:     }  
47:     printf("\nThe coefficients of cascaded second order sections "  
48:                    "of digital\nfilter are as follows...\n");  
49:     for(k = 0; k < N/2; k++)  
50:     {  
51:      printf("\nB%d = %f\tb%d0 = %f\ta%d0 = %f"  
52:                            ,k,B[k],k,b[k][0],k,a[k][0]);  
53:      printf("\n\t\tb%d1 = %f\ta%d1 = %f",k,b[k][1],k,a[k][1]);  
54:      printf("\n\t\tb%d2 = %f\ta%d2 = %f\n",k,b[k][2],k,a[k][2]);  
55:     }  
56:  getch();  
57:  }  
58:  /*          Butterworth filter design using Bilinear transformation  
59:  Enter the order of the filter N = 4  
60:  Enter the cutoff frequency of digital filter wc = 100  
61:  The coefficients of cascaded second order sections of digital  
62:  filter are as follows...  
63:  B0 = 0.085426  b00 = 1.000000 a00 = 1.000000  
64:            b01 = 2.000000 a01 = -2.140140  
65:            b02 = 1.000000 a02 = 1.481843  
66:  B1 = 0.129408  b10 = 1.000000 a10 = 1.000000  
67:            b11 = 2.000000 a11 = -3.242013  
68:            b12 = 1.000000 a12 = 2.759646  
69:   */