高通.cpp

butterworth高通滤波器的原程序

#include<iostream>
#include<cmath>
#include<fstream>
using namespace std;
const double pi=4*atan(1);
void butterworthhighpass(double,double,double,double,double);
int main()
{
	//数字高通滤波器性能要求
    double f,fc,fst;//f:抽样频率；fc:通带截止频率；fst:阻带截止频率
    double q1,q2;//q1:通带允许的最大衰减；q2:阻带允许的最小衰减
cout<<"please input the parameter:"<<endl;
	cin>>f>>fc>>fst>>q1>>q2;
	butterworthhighpass(f,fc,fst,q1,q2);
	return 0;}
void butterworthhighpass(double f,double fc,double fst,double q1,double q2)
{double Wc,Wst;//数字频率
     Wc=2*pi*fc/f;
	 Wst=2*pi*fst/f;
    double Uc=1;
	double C1=Uc*tan(Wc/2);
    double Ust=C1*(1/tan(Wst/2));
    double log(double num);
    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[20],A1[50],A2[50],B0[50],B1[50],B2[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;
		A2[k]=1;
		B0[k]=C1*C1+b[k]*C1+1;
		B1[k]=2*C1*C1-2;
		B2[k]=C1*C1+1-b[k]*C1;
	}
cout<<N<<endl;
		if (N%2==0)
		{
			for (k=1;k<=N/2;k++)
			outfile<<endl<<A0[k]/B0[k]<<" "<<A1[k]/B0[k]<<" "<<A2[k]/B0[k]<<" "<<1<<" "<<B1[k]/B0[k]<<" "<<B2[k]/B0[k]<<endl;
				}
		else 
		{
		outfile<<endl<<1<<" "<<-1<<" "<<1+C1<<" "<<C1-1<<endl;
for (k=1;k<=N/2;k++)
outfile<<endl<<A0[k]/B0[k]<<" "<<A1[k]/B0[k]<<" "<<A2[k]/B0[k]<<" "<<1<<" "<<B1[k]/B0[k]<<" "<<B2[k]/B0[k]<<endl;
		}
		outfile.close();}