Magnitude and Phase Response of System

Magnitude and Phase Response of  System

Problem Statement : Write a C program to plot the magnitude and phase response of a given system ( given: h(n): impulse response of system S) (Observe the frequency response for different systems. Compare the frequency response of a system (filter) for different length h(n) i.e filter coefficients).

1:  #include<stdio.h>  
2:  #include<conio.h>  
3:  #include<math.h>  
4:  #include<graphics.h>  
5:  void main()  
6:  {  
7:    double RealNum,ImagNum,RealDen,ImagDen,Real,Imag;  
8:    double mag[640],phase[640],pi,w,wStep;  
9:    float static a[10],b[10];  
10:    int M,N,i,k,gd=DETECT,gm;  
11:    clrscr();  
12:    printf("  \nThis program displays the magnitude and phase"  
13:         " transfer function plots\nfrom given coefficients of"  
14:         " difference equation \n\n");  
15:    printf("\nenter the number of coefficients of x(n), M = ");  
16:    scanf("%d",&M);  
17:    for(i = 0; i < M; i++)  
18:    {  
19:      printf("b%d = ",i);  
20:      scanf("%f",&b[i]);  
21:    }  
22:    printf("\nenter the number of coefficients of y(n), N = ");  
23:    scanf("%d",&N);  
24:    a[0] = 1.0;  
25:    for(i = 1; i <= N; i++)  
26:    {  
27:       printf("a%d = ",i);  
28:       scanf("%f",&a[i]);  
29:    }  
30:    pi = 22.0/7.0;  
31:    wStep = pi/640.0;  
32:    w=0;  
33:    for(k = 0; k <=640; k++)  
34:    {  
35:     RealNum = b[0];  
36:     ImagNum = 0;  
37:     for(i = 1; i < M; i++)  
38:     {  
39:        RealNum = RealNum + b[i]*cos(i*w);  
40:        ImagNum = ImagNum + b[i]*sin(i*w);  
41:     }  
42:     ImagNum = ImagNum * (-1.0);  
43:     RealDen = 0;  
44:     ImagDen = 0;  
45:     for(i = 0; i <= N; i++)  
46:     {  
47:        RealDen = RealDen + a[i]*cos(i*w);  
48:        ImagDen = ImagDen + a[i]*sin(i*w);  
49:     }  
50:     ImagDen = ImagDen * (-1.0);  
51:     mag[k] = 20*log10(sqrt(RealNum*RealNum + ImagNum*ImagNum)/  
52:               sqrt(RealDen*RealDen + ImagDen*ImagDen));  
53:     phase[k]=(180.0*7.0*(atan(ImagNum/RealNum)-atan(ImagDen/RealDen)))/22.0;  
54:     w = w + wStep;  
55:    }  
56:    initgraph(&gd,&gm,"e:\\tc\\bgi");  
57:    cleardevice();  
58:    setlinestyle(DOTTED_LINE,1,1);  
59:    line(0,150,640,150);  
60:    line(0,350,640,350);  
61:    line(10,0,10,350);  
62:    for(k = 0;k<640;k++)  
63:    {  
64:       putpixel(k,-mag[k]*2+150,WHITE);  
65:       putpixel(k,-phase[k]+350,WHITE);  
66:    }  
67:    outtextxy(500,200,"Magnitude plot");  
68:    outtextxy(500,450,"Phase plot");  
69:    getch();  
70:    closegraph();  
71:  }