北邮22级信通院DSP：用C++程序实现给定参数下四种滤波器的Butterworth模拟滤波器设计：给定上下截频和衰减系数求H（p）和H（s）

北邮22信通一枚~

跟随课程进度更新北邮信通院DSP的笔记、代码和文章，欢迎关注~

获取更多文章，请访问专栏：

北邮22级信通院DSP_青山入墨雨如画的博客-CSDN博客

一、核心算法

1.1判断滤波器类型

1.2 带通滤波器BP

1.3带阻滤波器BS

1.4综合四种滤波器算法

1.5展示函数show（）

1.6H（s）的显示（double_to_string）

1.6.1to_string方法

1.6.2 ostringstream方法

1.7自主输入函数input（）

二、代码部分

2.1总体代码

2.2 测试1

2.3测试2

一、核心算法

1.1判断滤波器类型

拿到通带和阻带的上下截频之后，我们首先应该判断滤波器的类型。

对带通带阻滤波器来说，如果s1<p1<p2<s2，则为带通滤波器；如果p1<s1<s2<p2，则为带阻滤波器。低通滤波器相当于通带下截频和阻带下截频都趋于负无穷，故对低通滤波器有p2<s2；同理对高通滤波器来说，相当于通带和阻带上截频趋于正无穷，故对高通滤波器有s1<p1。

同时需要考虑检测掉所有不合理的情况。

#define ERROR 0
#define IS_LP 1
#define IS_HP 2
#define IS_BP 3
#define IS_BS 4typedef long double  ld;//判断选通滤波器类型
int check_types(ld p1, ld p2, ld s1, ld s2)
{if (p1 > p2 || s1 > s2 || s1 == p2 ||s2 == p1 || s1 == s2 || p1 == p2)return ERROR;if (s1 * p1 != 0)return (p2 > p1 && s2 > s1 && p1 != s1 && p2 != s2) ?((p1 > s1) ? IS_BP : IS_BS) : (ERROR);elsereturn (s2 != p2) ? (s2 < p2 ? IS_HP : IS_LP) : (ERROR);
}

1.2 带通滤波器BP

首先检测对称情况，更新通阻带上下截频信息。之后计算butterworth滤波器的阶数N，之后展示H（p）的表达式。表达式的展示算法放在1.5展示函数show（）和1.6H（s）的显示（double_to_string）讲解。

#define ERROR 0
#define IS_LP 1
#define IS_HP 2
#define IS_BP 3
#define IS_BS 4typedef long double  ld;void BP(ld& p1, ld& p2, ld& Ap, ld& s1, ld& s2, ld& As)
{ld center = p1*p2;//保持通带cout << "通带中心频率为:" << sqrt(center) << "Hz;" << endl;//非中心对称改为中心对称if (p1 * p2 != s1 * s2)((center / s2) > s1) ? (s1 = (center / s2)) : (s2 = (center / s1));//计算N的值ld lambda_p = 1;ld lambda_s = (s2 - s1) / (p2 - p1);cout << "lambda_s的计算结果为:lambda_s = " << lambda_s << ";" << endl;double temp_1 = log10((pow(10, 0.1 * As) - 1) / (pow(10, 0.1 * Ap) - 1));double temp_2 = log10(lambda_s / lambda_p);int N = ceil(temp_1 / (2 * temp_2));//根据N的值查表得到低通滤波器归一化传输函数Hlp(p)分母表达式cout << "butterworth滤波器的阶数为:" << N << ";" << endl;show(N, "p");
}

1.3带阻滤波器BS

如上同理。需要注意参数修正位置，因为是带阻滤波器的设计，所以为了保证阻带，修正的是p1和p2的值。同时λs的计算公式也应该修改。

#define ERROR 0
#define IS_LP 1
#define IS_HP 2
#define IS_BP 3
#define IS_BS 4typedef long double  ld;void BS(ld& p1, ld& p2, ld& Ap, ld& s1, ld& s2, ld& As)
{ld center = s1 * s2;//保持通带cout << "通带中心频率为:" << sqrt(center) << "Hz;" << endl;//非中心对称改为中心对称if (p1 * p2 != s1 * s2)((center / p2) > p1) ? (p1 = (center / p2)) : (p2 = (center / p1));//计算N的值ld lambda_p = 1;ld lambda_s =  (p2 - p1)/ (s2 - s1);cout << "lambda_s的计算结果为:lambda_s=" << lambda_s << ";" << endl;double temp_1 = log10((pow(10, 0.1 * As) - 1) / (pow(10, 0.1 * Ap) - 1));double temp_2 = log10(lambda_s / lambda_p);int N = ceil(temp_1 / (2 * temp_2));//根据N的值查表得到低通滤波器归一化传输函数Hlp(p)分母表达式cout << "butterworth滤波器的阶数为:" << N << ";" << endl;show(N, "p");
}

1.4综合四种滤波器算法

如上同理。

首先根据1.1判断滤波器类型判断是何种滤波器，并计算中心频率center^2。对选通滤波器（高通滤波器和低通滤波器的合称）来说，由于侧重点不同，对中心频率的计算是不同的。带通滤波器（BP）center^2 = p1 * p2，而带阻滤波器（BS）应为center^2 = s1 * s2。

第二步：如果给定的通带阻带的上下截频不是自然几何对称的话，根据滤波器类型和中心频率修正相应的参数。

第三步，计算λs的值。由于带通滤波器λs的计算方式为λs = ((s2 - s1) / (p2 - p1))，而低通滤波器的为s/p，相当于低通滤波器是s1=p1=0的情况，λs的值可以用同一个式子处理；同理带阻滤波器和高通滤波器（s2=p2=0）的λs的值也可以共用带阻滤波器的式子处理λs = ((p2 - p1) / (s2 - s1))。

调用cmath头文件中的向上取整函数ceil，求得N值。

第四步，展示H（p）和H（s），详见1.5展示函数show（）和1.6H（s）的显示（double_to_string）

#define ERROR 0
#define IS_LP 1
#define IS_HP 2
#define IS_BP 3
#define IS_BS 4typedef long double  ld;
void compilation(ld& p1, ld& p2, ld& Ap, ld& s1, ld& s2, ld& As)
{ld center = 0;//判断滤波器类型if (check_types(p1, p2, s1, s2)){cout << endl << "该滤波器的类型为";switch (check_types(p1, p2, s1, s2)){case IS_LP:cout << "低通滤波器;" << endl; break;case IS_HP:cout << "高通滤波器;" << endl; break;case IS_BP:cout << "带通滤波器;" << endl; break;case IS_BS:cout << "带阻滤波器;" << endl; break;default:cout << "error" << endl; break;}center = (check_types(p1, p2, s1, s2) == IS_BP) ? p1 * p2 : s1 * s2;}if(sqrt(center))cout << endl << "通带中心频率为:" << sqrt(center) << "Hz;" << endl;//非几何对称改为几何对称if (p1 * p2 != s1 * s2)(check_types(p1, p2, s1, s2) == IS_BP) ?(((center / s2) > s1) ? (s1 = (center / s2)) : (s2 = (center / s1))) :(((center / p2) > p1) ? (p1 = (center / p2)) : (p2 = (center / p1)));//计算N的值ld lambda_p = 1;ld lambda_s = (check_types(p1, p2, s1, s2) == (IS_BP||IS_LP)) ?((s2 - s1) / (p2 - p1)) : ((p2 - p1) / (s2 - s1));	cout << endl << "归一化截频lambda_s的计算结果为:lambda_s = " << lambda_s << ";" << endl;double temp_1 = log10((pow(10, 0.1 * As) - 1) / (pow(10, 0.1 * Ap) - 1));double temp_2 = log10(lambda_s / lambda_p);int N = ceil(temp_1 / (2 * temp_2));//根据N的值查表得到低通滤波器归一化传输函数Hlp(p)分母表达式cout << endl << "butterworth滤波器的阶数为N = " << N << ";" << endl;cout << endl << "归一化传输函数为:H(p) = ";show(N, "p");string show_next_level = is_p(check_types(p1, p2, s1, s2), p2-p1, center);cout << endl << "传输函数H(s) = ";show(N, show_next_level);
}

1.5展示函数show（）

事先导入H（p）分母多项式系数的表。可以用二维数组存储。

输出表达式。

#define ERROR 0
#define IS_LP 1
#define IS_HP 2
#define IS_BP 3
#define IS_BS 4typedef long double  ld;
double box[7][8] =
{{1.0000,1.0000,0.00000,0.00000,0.00000,0.00000,0.0000,0.0000},{1.0000,1.4140,1.00000,0.00000,0.00000,0.00000,0.0000,0.0000},{1.0000,2.0000,2.00000,1.00000,0.00000,0.00000,0.0000,0.0000},{1.0000,2.6131,3.41420,2.61310,1.00000,0.00000,0.0000,0.0000},{1.0000,3.2361,5.23610,5.23610,6.23610,1.00000,0.0000,0.0000},{1.0000,3.8637,7.46410,9.14160,7.46410,3.86370,1.0000,0.0000},{1.0000,4.4940,10.0978,14.5918,14.5918,10.0978,4.4940,1.0000}
};void  show(int N, string p)
{cout << "1 / ( ";for (int i = 7; i >= 0; i--){if(box[N - 1][i])//系数不为0时有输出{if (box[N - 1][i] != 1)cout << box[N - 1][i] << "*";switch (i){case 0:cout << 1; break;case 1:cout << p << " + "; break;default:cout << p << "^" << i << " + "; break;}}}cout << " );" << endl;
}

1.6H（s）的显示（double_to_string）

用p=q（s）替代上述show函数中的“p”。带入中心频率。关键问题是浮点数如何转化为字符串。

给出两种解决方法：

1.6.1to_string方法

需要cstring头文件。

#include<cstring>
string is_p_1(int type, ld B, ld& center)
{string output;if (type == IS_LP)output = "(s/" + to_string(B) + ")";else if (type == IS_HP)output = "(" + to_string(B) + "/s";else if (type == IS_BP)output = "((s^2+" + to_string(pow(center, 2)) + ")/" +"(" + to_string(B) + "*s" + "))";else if (type == IS_BS)output = "((" + to_string(B) + "*s" + ")/" +"(s^2+" + to_string(pow(center, 2)) + "))";elseoutput = "error";return output;
}

1.6.2 ostringstream方法

需要sstream头文件。

#include <sstream>
#include<iomanip>
using namespace std;
auto format_doble_value(double val, int fixed) {std::ostringstream oss;oss << std::setprecision(fixed) << val;return oss.str();
}string is_p(int type, ld B, ld& center)
{string output;if (type == IS_LP)output = "(s/" + format_doble_value(B, define_setpercision) + ")";else if (type == IS_HP)output = "(" + format_doble_value(B, define_setpercision) + "/s";else if (type == IS_BP)output = "((s^2+" + format_doble_value(pow(center, 2), define_setpercision) + ")/" +"(" + format_doble_value(B, define_setpercision) + "*s" + "))";else if (type == IS_BS)output = "((" + format_doble_value(B, define_setpercision) + "*s" + ")/" +"(s^2+" + format_doble_value(pow(center, 2), define_setpercision) + "))";elseoutput = "error";return output;
}

1.7自主输入函数input（）

void input()
{cout << endl << "通带起始频率p1:"; cin >> p1;cout << endl << "通带截止频率p2:"; cin >> p2;cout << endl << "通带衰减系数(dB):"; cin >> Ap;cout << endl << "阻带起始频率s1:"; cin >> s1;cout << endl << "阻带截止频率s2:"; cin >> s2;cout << endl << "阻带衰减系数(dB):"; cin >> As;cout << endl;cout << "运算结果为:" << endl;
}

二、代码部分

2.1总体代码

#include<iostream>
#include<cmath>
#include<cstring>
#include <sstream>
#include<iomanip>
using namespace std;#define ERROR 0
#define IS_LP 1
#define IS_HP 2
#define IS_BP 3
#define IS_BS 4typedef long double  ld;
const double PI = acos(-1.0);
const int define_setpercision = 5;
double box[7][8] =
{{1.0000,1.0000,0.00000,0.00000,0.00000,0.00000,0.0000,0.0000},{1.0000,1.4140,1.00000,0.00000,0.00000,0.00000,0.0000,0.0000},{1.0000,2.0000,2.00000,1.00000,0.00000,0.00000,0.0000,0.0000},{1.0000,2.6131,3.41420,2.61310,1.00000,0.00000,0.0000,0.0000},{1.0000,3.2361,5.23610,5.23610,6.23610,1.00000,0.0000,0.0000},{1.0000,3.8637,7.46410,9.14160,7.46410,3.86370,1.0000,0.0000},{1.0000,4.4940,10.0978,14.5918,14.5918,10.0978,4.4940,1.0000}
};ld p1, p2, s1, s2, Ap, As;//判断选通滤波器类型
int check_types(ld p1, ld p2, ld s1, ld s2)
{if (p1 > p2 || s1 > s2 || s1 == p2 ||s2 == p1 || s1 == s2 || p1 == p2)return ERROR;if (s1 * p1 != 0)return (p2 > p1 && s2 > s1 && p1 != s1 && p2 != s2) ?((p1 > s1) ? IS_BP : IS_BS) : (ERROR);elsereturn (s2 != p2) ? (s2 < p2 ? IS_HP : IS_LP) : (ERROR);
}auto format_doble_value(double val, int fixed) {std::ostringstream oss;oss << std::setprecision(fixed) << val;return oss.str();
}string is_p(int type, ld B, ld& center)
{string output;if (type == IS_LP)output = "(s/" + format_doble_value(B, define_setpercision) + ")";else if (type == IS_HP)output = "(" + format_doble_value(B, define_setpercision) + "/s";else if (type == IS_BP)output = "((s^2+" + format_doble_value(pow(center, 2), define_setpercision) + ")/" +"(" + format_doble_value(B, define_setpercision) + "*s" + "))";else if (type == IS_BS)output = "((" + format_doble_value(B, define_setpercision) + "*s" + ")/" +"(s^2+" + format_doble_value(pow(center, 2), define_setpercision) + "))";elseoutput = "error";return output;
}string is_p_1(int type, ld B, ld& center)
{string output;if (type == IS_LP)output = "(s/" + to_string(B) + ")";else if (type == IS_HP)output = "(" + to_string(B) + "/s";else if (type == IS_BP)output = "((s^2+" + to_string(pow(center, 2)) + ")/" +"(" + to_string(B) + "*s" + "))";else if (type == IS_BS)output = "((" + to_string(B) + "*s" + ")/" +"(s^2+" + to_string(pow(center, 2)) + "))";elseoutput = "error";return output;
}//标准输出
void  show(int N, string p)
{cout << "1 / ( ";for (int i = 7; i >= 0; i--){if(box[N - 1][i])//系数不为0时有输出{if (box[N - 1][i] != 1)cout << box[N - 1][i] << "*";switch (i){case 0:cout << 1; break;case 1:cout << p << " + "; break;default:cout << p << "^" << i << " + "; break;}}}cout << " );" << endl;
}//带通滤波器算法
void BP(ld& p1, ld& p2, ld& Ap, ld& s1, ld& s2, ld& As)
{ld center = p1*p2;//保持通带cout << "通带中心频率为:" << sqrt(center) << "Hz;" << endl;//非中心对称改为中心对称if (p1 * p2 != s1 * s2)((center / s2) > s1) ? (s1 = (center / s2)) : (s2 = (center / s1));//计算N的值ld lambda_p = 1;ld lambda_s = (s2 - s1) / (p2 - p1);cout << "lambda_s的计算结果为:lambda_s = " << lambda_s << ";" << endl;double temp_1 = log10((pow(10, 0.1 * As) - 1) / (pow(10, 0.1 * Ap) - 1));double temp_2 = log10(lambda_s / lambda_p);int N = ceil(temp_1 / (2 * temp_2));//根据N的值查表得到低通滤波器归一化传输函数Hlp(p)分母表达式cout << "butterworth滤波器的阶数为:" << N << ";" << endl;show(N, "p");
}void BS(ld& p1, ld& p2, ld& Ap, ld& s1, ld& s2, ld& As)
{ld center = s1 * s2;//保持通带cout << "通带中心频率为:" << sqrt(center) << "Hz;" << endl;//非中心对称改为中心对称if (p1 * p2 != s1 * s2)((center / p2) > p1) ? (p1 = (center / p2)) : (p2 = (center / p1));//计算N的值ld lambda_p = 1;ld lambda_s =  (p2 - p1)/ (s2 - s1);cout << "lambda_s的计算结果为:lambda_s=" << lambda_s << ";" << endl;double temp_1 = log10((pow(10, 0.1 * As) - 1) / (pow(10, 0.1 * Ap) - 1));double temp_2 = log10(lambda_s / lambda_p);int N = ceil(temp_1 / (2 * temp_2));//根据N的值查表得到低通滤波器归一化传输函数Hlp(p)分母表达式cout << "butterworth滤波器的阶数为:" << N << ";" << endl;show(N, "p");
}void compilation(ld& p1, ld& p2, ld& Ap, ld& s1, ld& s2, ld& As)
{ld center = 0;//判断滤波器类型if (check_types(p1, p2, s1, s2)){cout << endl << "该滤波器的类型为";switch (check_types(p1, p2, s1, s2)){case IS_LP:cout << "低通滤波器;" << endl; break;case IS_HP:cout << "高通滤波器;" << endl; break;case IS_BP:cout << "带通滤波器;" << endl; break;case IS_BS:cout << "带阻滤波器;" << endl; break;default:cout << "error" << endl; break;}center = (check_types(p1, p2, s1, s2) == IS_BP) ? p1 * p2 : s1 * s2;}if(sqrt(center))cout << endl << "通带中心频率为:" << sqrt(center) << "Hz;" << endl;//非几何对称改为几何对称if (p1 * p2 != s1 * s2)(check_types(p1, p2, s1, s2) == IS_BP) ?(((center / s2) > s1) ? (s1 = (center / s2)) : (s2 = (center / s1))) :(((center / p2) > p1) ? (p1 = (center / p2)) : (p2 = (center / p1)));//计算N的值ld lambda_p = 1;ld lambda_s = (check_types(p1, p2, s1, s2) == (IS_BP||IS_LP)) ?((s2 - s1) / (p2 - p1)) : ((p2 - p1) / (s2 - s1));	cout << endl << "归一化截频lambda_s的计算结果为:lambda_s = " << lambda_s << ";" << endl;double temp_1 = log10((pow(10, 0.1 * As) - 1) / (pow(10, 0.1 * Ap) - 1));double temp_2 = log10(lambda_s / lambda_p);int N = ceil(temp_1 / (2 * temp_2));//根据N的值查表得到低通滤波器归一化传输函数Hlp(p)分母表达式cout << endl << "butterworth滤波器的阶数为N = " << N << ";" << endl;cout << endl << "归一化传输函数为:H(p) = ";show(N, "p");string show_next_level = is_p(check_types(p1, p2, s1, s2), p2-p1, center);cout << endl << "传输函数H(s) = ";show(N, show_next_level);
}void input()
{cout << endl << "通带起始频率p1:"; cin >> p1;cout << endl << "通带截止频率p2:"; cin >> p2;cout << endl << "通带衰减系数(dB):"; cin >> Ap;cout << endl << "阻带起始频率s1:"; cin >> s1;cout << endl << "阻带截止频率s2:"; cin >> s2;cout << endl << "阻带衰减系数(dB):"; cin >> As;cout << endl;cout << "运算结果为:" << endl;
}int main()
{system("color 0A");//输入六个参量input();//例子//p1 = 3.1e6; p2 = 5.5e6; Ap = 3; s1 = 3.8e6; s2 = 4.8e6; As = 20;//p1 = 0; p2 = 750; Ap = 3; s1 = 0; s2 = 1600; As = 7;p1 *= 2 * PI; p2 *= 2 * PI; s1 *= 2 * PI; s2 *= 2 * PI;compilation(p1, p2, Ap, s1, s2, As);return 0;
}