#include "math.h"
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#define MAX 2500
#define PI 3.1415926
/****************************************************************************************
* biliner(int fc,int N,int fs,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,200,2,x);
for(i=0;i<MAX;i++)
fprintf(fpx,"%f\n",x[i]);//输出x[]
fclose(fpx);
fclose(fp);
}
巴特沃斯滤波器
最新推荐文章于 2024-08-12 21:26:33 发布