#include "math.h" 
#include "stdio.h"
int		N=50;
int		j=4;						  				   /* IIR濾波器階數(shù) */
float   b[4]={0.0238018,0.0714055,0.0714055,0.0238018};/* IIR濾波器H(z)系數(shù)b */
float   a[4]={1,-1.6217269,1.0505308,-0.2383891};	   /* IIR濾波器H(z)系數(shù)a */
float   w=3.1415926/50;   
int 	i,k;
float   H,Ra=0,Rb=0,Ia=0,Ib=0;
float   real_H[50]; 
float   imag_H[50];    
float   db[50];

void main()
{  
	for(i=0;i<N;i++)               				       /* 由濾波器階數(shù)與系數(shù)求H(e^jw) */
	{ 	
		for(k=0;k<j;k++)               
		{ 	
		  Rb=Rb+b[k]*cos(-w*i*k);
		  Ib=Ib+b[k]*sin(-w*i*k);
		  Ra=Ra+a[k]*cos(-w*i*k);
		  Ia=Ia+a[k]*sin(-w*i*k);		
		}
		real_H[i]=(Ra*Rb+Ia*Ib)/(Ra*Ra+Ia*Ia);
		imag_H[i]=(Ib*Ra-Rb*Ia)/(Ra*Ra+Ia*Ia);
		Ra=Rb=Ia=Ib=0;		
	}
	
	for(i=0;i<N;i++)								   /* 求對數(shù)幅頻特性 */
	{
		H=sqrt(real_H[i]*real_H[i]+imag_H[i]*imag_H[i]);
		db[i]=20.0*log10(H);
	}
}