firwin.h

三种模拟调制与解调的算法

#include "head.h"
static double bessel0(double x)
{ 
	int i;
	double d,y,d2,sum;
	y=x/2.0;
	d=1.0;
	sum=1.0;
	for(i=1;i<=25;i++)
	{
		d=d*y/i;
		d2=d*d;
		sum=sum+d2;
		if(d2<sum*(1.0e-8))
			break;
	}
	return(sum);
}

static double kaiser(int i,int n,double beta)
{  
	double a,w,a2,b1,b2,beta1;
//	double bessel0(double x);
	b1=bessel0(beta);
	a=2.0*i/(n-1)-1.0;
	a2=a*a;
	beta1=beta*sqrt(1.0-a2);
	b2=bessel0(beta1);
	w=b2/b1;
	return(w);
}

static double window(int type,int n,int i,double beta)
{
	int k;
	double w;
	//double kaiser(int i,int n,double beta);
	w=1.0;
	switch(type)
	{
        case 1:
			{ w=1.0;
			  break;
			}
		case 2:
			{ k=(n-2)/10;
			  if(i<=k)
				  w=0.5*(1.0-cos(i*pi/(k+1)));
			  if(i>n-k-2)
				  w=0.5*(1.0-cos(n-i-1)*pi/(k+1));
			  break;
			}
		case 3:
			{ w=1.0-fabs(1.0-2*i/(n-1.0));
			  break;
			}
		case 4:
			{ w=0.5*(1.0-cos(2*i*pi/(n-1)));
			  break;
			}
		case 5:
			{ w=0.54-0.46*cos(2*i*pi/(n-1));
			  break;
			}
		case 6:
			{ w=0.42-0.5*cos(2*i*pi/(n-1))+0.08*cos(4*i*pi/(n-1));
			  break;
			}
		case 7:
			{ w=kaiser(i,n,beta);
			  break;
			}
	}
	return(w);
}


float *firwin(int n,int band,double fln,double fhn,int wn,double beta)
{
 //n	滤波器的阶数
 //band 滤波器类型，1--低通，2--高通,3--带通,4--带阻
 //fln,对于低通和高通滤波器：通带边界频率，对于带通和带阻
 //fln:通带下边界频率，fhn:通带上边界频率
 //wn  窗函数类型，1--矩形窗，2--图基窗，3--三角窗
 //4--汉宁窗，5--海明窗，6--布拉克曼窗，7--凯塞窗
 //beta凯塞窗的因子
 //h   一个指针，指向返回的h(n)
 //若失败，则返回一个空指针(NULL)
 int i,n2,mid;
 double s,wc1,wc2,delay;
 //3<beta<10;
 float *h;
 h=(float *) malloc((n+1)*sizeof(float));
 if(!h) return NULL;
 if((n%2)==0)
 {
	 n2=n/2-1;
	 mid=1;
 }
 else
 {
	 n2=n/2;
	 mid=0;
 }
 delay=n/2.0;
 wc1=2.0*pi*fln;
 if(band>=3)
	 wc2=2.0*pi*fhn;
 switch(band)
 {  
    case 1:
		{ 
			for(i=0;i<=n2;i++)
			{ s=i-delay;
			  *(h+i)=(sin(wc1*s)/(pi*s))*window(wn,n+1,i,beta);
			  *(h+n-i)=*(h+i);
			}
			if(mid==1)
		      *(h+n/2)=wc1/pi;
			break;
		}
	case 2:
		{ 
			for(i=0;i<=n2;i++)
			{ s=i-delay;
			  *(h+i)=(sin(pi*s)-sin(wc1*s))/(pi*s);
			  *(h+i)=*(h+i)*window(wn,n+1,i,beta);
			  *(h+n-i)=*(h+i);
			}
			if(mid==1)
		      *(h+n/2)=1.0-wc1/pi;
			break;
		}
    case 3:
		{ 
			for(i=0;i<=n2;i++)
			{ s=i-delay;
			  *(h+i)=(sin(wc2*s)-sin(wc1*s))/(pi*s);
			  *(h+i)=*(h+i)*window(wn,n+1,i,beta);
			  *(h+n-i)=*(h+i);
			}
			if(mid==1)
		      *(h+n/2)=(wc2-wc1)/pi;
			break;
		}
   case 4:
		{ 
			for(i=0;i<=n2;i++)
			{ s=i-delay;
			  *(h+i)=(sin(wc1*s)+sin(pi*s)-sin(wc2*s))/(pi*s);
			  *(h+i)=*(h+i)*window(wn,n+1,i,beta);
			  *(h+n-i)=*(h+i);
			}
			if(mid==1)
		      *(h+n/2)=(wc1+pi-wc2)/pi;
			break;
		}
 }
 return h;
}