用窗函数法设计FIR滤波器
#include <stdio.h>
#include <dos.h>
#include <graphics.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#include <bios.h>
#define NFFT 512 /*length of fft*/
float ino(float x);
void cfft(float x[],float y[],int M,int N,int lc);
void plot(float h[],int n,int color);
int pow2(int k);
void main()
{
float hd[NFFT],h[NFFT],w[NFFT],wc,pi,him[NFFT],R[NFFT];
float a,b,beta,bes,g,q,p,wf;
int m,n,k,i;
char ch;
printf("this is a fir system design program.¥n");
printf(" m represents window function.¥n");
printf(" 1...... rectangular window function.¥n");
printf(" 2...... kaiser window function.¥n");
loop:printf("choose m, m is 1 or 2¥n m=¥n");
scanf("%d",&m);
printf("choose window length n ...21,51,101,201.¥n n=¥n");
scanf("%d",&n);
printf("choose wc,wc is between 0.1 to 0.99¥n wc=¥n");
scanf("%f",&wc);
a=(n-1)/2;
pi=4.0*atan(1.0);
for(i=0;i<n;i++)
{
if(i==a)
hd[i]=wc;
else
{b=i-a;
hd[i]=sin(pi*b*wc)/(pi*b);}
}
switch(m)
{
case 1:printf("rectangular window function.¥n");
for(i=0;i<n;i++)
w[i]=1.0;
break;
case 2:printf("kaiser window function.¥n");
printf("choose beta,beta is between 4 to 9.¥n beta=¥n");
scanf("%f",&beta);
bes=ino(beta);
for(i=0;i<n;i++)
{
g=1.0-pow(1.0-2.0*(float)i/(float)(n-1),2);
q=beta*sqrt(g);
w[i]=ino(q);
w[i]=w[i]/bes;
}
break;
default:printf("m is wrong,please choose m again.¥n");
goto loop;
}
for(i=0;i<n;i++)
{
h[i]=hd[i]*w[i];
printf("h(%2d)=%14.6f¥n",i,h[i]);
}
printf("¥nPlease Wait...¥n");
for(i=0;i<n;i++)
him[i]=0.0;
cfft(h,him,n,NFFT,0);
for(i=0;i<NFFT/2;i++)
{
R[i]=h[i]*h[i]+him[i]*him[i];
if(R[i]==0.0) R[i]=1.0e-16;
R[i]=10*log10(R[i]);
}
plot(R,NFFT/2,15);
getch();
closegraph();
}
/*********************************************/
float ino(float x)
{
float y,t,e,z,de,sde;
int i;
y=x/2.0;t=1.0e-08;e=1.0;de=1.0;
for(i=1;i<=25;i++)
{
de=de*y/i;sde=de*de;e=e+sde;
if(e*t<=sde) ;
else
z=e;
}
return (z);
}
/*------*/
/* FFT subroutine : lc=0 forword transform */
/* lc=1 ifft transform */
/* ------N is FFT's Length ,M is length of x----*/
void cfft(x,y,M,N,lc)
float x[],y[];
int M,N,lc;
{
float t1,t2,PI=3.1415926,u1,u2,e,a;
int NN,n1,n2,j,i,k,l,m;
j=0;
n1=N-2;
for(i=M;i<N;i++)
x[i]=y[i]=0.0;
if(lc==1)
for(i=0;i<N;i++)
y[i]=-y[i];
for(i=0;i<=n1;i++)
{ if(i<j)
{ t1=x[j];t2=y[j];
x[j]=x[i];y[j]=y[i];
x[i]=t1;y[i]=t2; }
k=N/2;
while(k<=j)
{ j=j-k;k=k/2; }
j=j+k;
}
/*main fft loop*/
n2=1;
NN=N;
m=pow2(NN);
for(k=1;k<=m;k++)
{ n1=n2;
n2=2*n2;
e=PI/n1;
a=0;
for(j=0;j<=n1-1;j++)
{
u1=cos(a);u2=-1*sin(a);
a=(j+1)*e;
for(i=j;i<=N-1;i+=n2)
{ l=i+n1;
t1=x[l]*u1-y[l]*u2;t2=x[l]*u2+y[l]*u1;
x[l]=x[i]-t1;y[l]=y[i]-t2;
x[i]=x[i]+t1;y[i]=y[i]+t2; }
}
}
if(lc==1)
for(i=0;i<N;i++)
{ x[i]=x[i]/(float)N;
y[i]=-y[i]/(float)N; }
}
int pow2(int k)
{ int k1=0;
while(k!=1)
{ k=k/2;k1=k1+1; }
return k1;
}
/*------*/
void plot(float h[],int n,int color)
{ int *Graphdriver=DETECT,*Graphmode;
int i,k,x,y,am;
int zbw[500];
char ss[10];
initgraph(Graphdriver,Graphmode,"c:¥¥borlandc¥¥bgi");
settextstyle(DEFAULT_FONT,HORIZ_DIR,0);
for(i=1;i<10;i++)
{ am=0-(i-1)*20;
itoa(am,ss,10);
if(am==0) outtextxy(184,78+(i-1)*40,ss);
else if(am>-100) outtextxy(168,78+(i-1)*40,ss);
else outtextxy(160,78+(i-1)*40,ss);
delay(2); }
for(i=1;i<=5;i++)
{ itoa(i,ss,10);
outtextxy(185+i*52,410,"0.");
outtextxy(198+i*52,410,ss); }
setlinestyle(DOTTED_LINE,0,1);
for(i=1;i<10;i++)
{ line(200,80+(i-1)*40,460,80+(i-1)*40);
delay(2); }
for(i=1;i<=11;i++)
{ line(200+(i-1)*26,80,200+(i-1)*26,400);
delay(2); }
setlinestyle(SOLID_LINE,0,1);
rectangle(200,80,460,400);
for(i=0;i<n;i++)
zbw[i]=(int)h[i];
setcolor(color);
moveto(200,80-2*zbw[0]);
for(k=1;k<n;k++)
{ x=200+k;
y=80-(2*zbw[k]);
if(y<450)
{
lineto(x,y);
lineto(x,y-1);
}
delay(10);
}
}