北邮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 4
typedef 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);
else
return (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 4
typedef 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 4
typedef 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 4
typedef 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 4
typedef 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)) + "))";
else
output = "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) + "))";
else
output = "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 4
typedef 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);
else
return (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) + "))";
else
output = "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)) + "))";
else
output = "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;
}
2.2 测试1
带阻滤波器:p1 = 3.1e6; p2 = 5.5e6; Ap = 3; s1 = 3.8e6; s2 = 4.8e6; As = 20;
2.3测试2
低通滤波器:p1 = 0; p2 = 750; Ap = 3; s1 = 0; s2 = 1600; As = 7;