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: }
