fir.cpp

设计一FIR低通滤波器 0=<|w|<=0.25Pi,输入n代表FIR滤波器的长度,输出H(ejw)幅度值,N代表输出数,Hbs[N]存放H(ejw)的幅度值,h[n]存放H[z]的系数

/*====================程序描述=======================*/
/*  设计一FIR低通滤波器，0=<|w|<=0.25Pi              */
/*  输入n代表FIR滤波器的长度                         */
/*  输出H(ejw)幅度值,N代表输出数                     */
/*  Hbs[N]存放H(ejw)的幅度值                         */
/*  h[n]存放H[z]的系数                               */
/*  算法采用Hamming窗函数                            */
/*===================================================*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define  N   50

typedef struct complex
{
	double real;
	double imag;
}   COMPLEX;

COMPLEX   mul(COMPLEX,COMPLEX)                        ;  /*复数相乘*/
COMPLEX   mul_n(COMPLEX ,int )                        ;  /*计算复数的n次幂*/
double    mabs(COMPLEX)                               ;  /*求幅度*/
COMPLEX   e_jw(double)                                ;  /*求复指数*/


int main(void)
{
	int n,i,j;
	COMPLEX H[N+1],Ej,Ejw_n;
	double *h1,*h,*w,Pi,dw=0,Hbs[N+1];
	Pi=atan(1)*4;
	printf("please input the FIR's length:");
	scanf("%d$\n",&n);
    h1=(double*)calloc(n,sizeof(double));
	h=(double*)calloc(n,sizeof(double));
	w=(double*)calloc(n,sizeof(double));
	for (i=0;i<=n;++i)
	{   
		if(i==n/2)
			h1[i]=0.25;
		else 
			h1[i]=sin((i-n/2)*0.25*Pi)/(Pi*(i-n/2));
		w[i]=0.54-0.46*cos(2*Pi*i/n);
	}
	for (i=0;i<=n;++i)
		h[i]=h1[i]*w[i];
	
    printf("the h[n] is as follows:\n");
	for (i=0;i<=n;++i)
	{
		if(i%5==0)
			printf("\n");
		printf("%.4lf ",h[i]);
	}
	
	for(i=0;i<=N;++i,dw+=Pi/N)
	{
		H[i].real=0;
		H[i].imag=0;
		for(j=0;j<=n;++j)
		{
			Ej=e_jw(dw);
			Ejw_n=mul_n(Ej,j);
			H[i].real+=h[j]*Ejw_n.real;
			H[i].imag+=h[j]*Ejw_n.imag;
		}
		
	}
    		
    for(i=0;i<=N;++i)
		Hbs[i]=mabs(H[i]);

	printf("\nthe H[ejw] is as follows:\n");                 /*输出幅度值*/
	for (i=0;i<=N;++i)
	{
		if(i%5==0)
			printf("\n");
		printf("%.4lf ",Hbs[i]);
	}
    printf("\n");
	return 0;
}



COMPLEX   e_jw(double w)
{
	COMPLEX Ejw;
	Ejw.real=cos(w);
	Ejw.imag=-sin(w);
	return Ejw;
}

COMPLEX mul(COMPLEX a,COMPLEX b)
{
	COMPLEX c;
	c.real=a.real*b.real-a.imag*b.imag;
	c.imag=a.real*b.imag+a.imag*b.real;
	return c;
}

COMPLEX mul_n(COMPLEX a,int b)
{
	int i;
	COMPLEX tem,tem1;
	if(b==0)
		{
		   tem.real=1.0;
	       tem.imag=0  ;
		}
	else
	{
		tem.real=a.real;
		tem.imag=a.imag;
        for(i=1;i<b;++i)
		{
			tem1=mul(a,tem);
			tem.real=tem1.real;
			tem.imag=tem1.imag;
		}
	}
	return tem;
}

double mabs(COMPLEX a)
{
 double m;
 m=a.real*a.real+a.imag*a.imag;
 m=sqrt(m);
 return m;
}