带阻滤波器.cpp

用c++编的巴特沃斯带阻滤波器原程序

#include<iostream>
#include<cmath>
#include<fstream>
using namespace std;
const double pi=4*atan(1);
void butterworthbandstop(double,double,double,double,double,double,double);
int main()
{
	//数字带阻滤波器性能要求
    double f=1000,f1=100,f2=400,fst1=200,fst2=350;//f:抽样频率；f1,f2:边带频率；fst1,fst2:阻带截止频率
    double q1=1,q2=30;//q1:通带允许的最大衰减；q2:阻带允许的最小衰减
    butterworthbandstop(f,f1,f2,fst1,fst2,q1,q2);
	return 0;}
void butterworthbandstop(double f,double f1,double f2,double fst1,double fst2,double q1,double q2)
{	double w1,w2,wst1,wst2;//数字频率
     w1=2*pi*f1/f;
	 w2=2*pi*f2/f;
	 wst1=2*pi*fst1/f;
	 wst2=2*pi*fst2/f;
    double Uc=1;
	double D1=Uc*tan((w2-w1)/2);
    double E1=2*cos((w2+w1)/2)/cos((w2-w1)/2);
    double Ust1=D1*sin(wst1)/(cos(wst1)-E1/2);
	double Ust2=D1*sin(wst2)/(cos(wst2)-E1/2);
	double Ust;
	if (fabs(Ust1)<fabs(Ust2))
		Ust=fabs(Ust1);
	else Ust=fabs(Ust2);
    double v;
    v=log10(pow(10,q2/10)-1)/(2*log10(Ust));//计算阶数
    double ceil(double x);
   int N=ceil(v);//阶数圆整
	int k;
	double b[20];
	for(k=1;k<=(N/2);k++)
	b[k]=-2*cos(pi/2+(2*k-1)*pi/(2*N));
	double A0[50],A1[50],A2[50],A3[50],A4[50],B0[50],B1[50],B2[50],B3[50],B4[50];//数字低通滤波器系统函数参数
	ofstream outfile("f1.dat",ios::out);//定义文件流对象，打开磁盘文件“f1.dat"
	if(!outfile)
	{ 
	   cerr<<"open error!"<<endl;
	exit(1);
	}
	for (k=1;k<=N/2;k++)
	{
		A0[k]=1;
		A1[k]=-2*E1;
		A2[k]=2+E1*E1;
		A3[k]=-2*E1;
		A4[k]=1;
		B0[k]=D1*D1+b[k]*D1+1;
		B1[k]=-b[k]*D1*E1-2*E1;
		B2[k]=-2*D1*D1+2+E1*E1;
		B3[k]=b[k]*D1*E1-2*E1;
		B4[k]=D1*D1-b[k]*D1+1;
	}
cout<<N<<endl;
		if (N%2==0)
		{
			for (k=1;k<=N/2;k++)
			outfile<<endl<<A0[k]<<" "<<A1[k]<<" "<<A2[k]<<" "<<A3[k]<<" "<<A4[k]<<" "<<B0[k]<<" "<<B1[k]<<" "<<B2[k]<<" "<<B3[k]<<" "<<B4[k]<<endl;
				}
		else 
		{
		outfile<<1<<" "<<-E1<<" "<<1<<" "<<D1+1<<" "<<-E1<<" "<<1-D1<<endl;
for (k=1;k<=N/2;k++)
outfile<<A0[k]<<" "<<A1[k]<<" "<<A2[k]<<" "<<A3[k]<<" "<<A4[k]<<" "<<B1[k]<<" "<<B2[k]<<" "<<B3[k]<<" "<<B4[k]<<endl;
		}
		outfile.close();}