巴特沃斯滤波器

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

#define MAX 2500
#define PI 3.1415926

/****************************************************************************************
* biliner(int fc,fs,int N,float*x)
*  参数说明  fc 截止频率
*             fs    采样率
*     N     (1=<N<=3)  阶数
*     *x    传入的待滤波数组
*双线性法基于巴特沃斯的数字滤波器    
****************************************************************************************/

using namespace std;

void biliner(int fc,int N,int fs,float *x)
{
	 //定义变量
	 long i=0;   
	 float W,Wc,t,y[MAX];//MAX是3000个点
	 //转化为数字角频率
	 W=2*PI*fc/fs;
	 //预畸变
	 Wc=2*fs*tan(w/2);
	 //去归一化影响
	 t=2*fs/wc;
     switch(N)
     {
	  /*一阶情况*/
	  case 1:
		    //Hs=1/(p+1);
		    y[0]=x[0]/(1+t);
		    for(i=1;i<MAX;i++)
		     y[i]=(x[i]+x[i-1]-(1-t)*y[i-1])/(1+t);
		    break;
		    /*二阶情况*/
	   case 2: 
		    //Hs=1/(p*p+1.414*p+1);
		    y[0]=x[0]/(t*t+1.414*t+1);
		    y[1]=(x[1]+2*x[0]-(2-2*t*t)*y[0])/(t*t+1.414*t+1);
		    for(i=2;i<MAX;i++)
		     y[i]=(x[i]+2*x[i-1]+x[i-2]-(2-2*t*t)*y[i-1]-(1-1.414*t+t*t)*y[i-2])/(t*t+1.414*t+1);
		    break;
		    /*三阶*/
	   case 3:
		    y[0]=x[0]/(t*t*t+2*t*t+2*t+1);
		    y[1]=(x[1]+3*x[0]-(3-3*t*t*t-2*t*t+2*t)*y[0])/(t*t*t+2*t*t+2*t+1);
		    y[2]=(x[2]+3*x[1]+3*x[0]-(3-3*t*t*t-2*t*t+2*t)*y[1]-(3*t*t*t-2*t*t-2*t+3)*y[0])/(t*t*t+2*t*t+2*t+1);
		    for(i=3; i < MAX ;i++)
		     y[i]=(x[i]+3*x[i-1]+3*x[i-2]+x[i-3]-(3-3*t*t*t-2*t*t+2*t)*y[i-1]-(3*t*t*t-2*t*t-2*t+3)*y[i-2]-(1-t*t*t+2*t*t-2*t)*y[i-3])/(t*t*t+2*t*t+2*t+1);
		       break;
	    default:
			    return ;
				break;
	  }
	  
	  for(i=0;i<MAX;i++)
	    x[i]=y[i];  //输出x[] 
	    
	 return;
}
int main(int argc, char *argv[])
{
	 float x[MAX];
	 int i,lenth,m;
	 FILE *fpx=fopen("newdata.txt","w");
	 FILE *fp=fopen("ChannelData.txt","r");
	 if(fp==NULL||fpx==NULL)
	 {
	 printf("失败!!!\n");
	 return 0;
	 }
	 for(i=0;i<MAX;i++)
	 fscanf(fp,"%f ",&x[i]);//输入x[] 
	
	 biliner(2,2,200,x);
	 
	 for(i=0;i<MAX;i++)
	  fprintf(fpx,"%f\n",x[i]);//输出x[] 
	  fclose(fpx);
	  fclose(fp);

}