北邮22信通一枚~
跟随课程进度更新北邮信通院DSP的笔记、代码和文章,欢迎关注~
获取更多文章,请访问专栏:
北邮22级信通院DSP_青山入墨雨如画的博客-CSDN博客
承接上一篇博客
北邮22级信通院DSP:IIR_DF系统2.0版:用C++程序自主写题——第四章IIR_DF设计的大部分题目——给定参量和选用的方法,设计出IIR_DF,一直进行到得出H(s)-CSDN博客
在原有的基础上,应对通带衰减系数不是3dB的情况,提出了一种更为严谨精确的计算模式。
先来展示运行效果
一、更为严谨精确的complication()函数
因为对butterworth滤波器来说,表达式为:
所以如果用butterworth滤波器进行拟合,就需要确定参量N和Ωc,其中Ωc的意思是通带3dB截止频率。没有引入频率归一化定义时,对N的计算公式应为:
频率归一化的工作,主要是为了“查表法”操作的可实现性。λs定义为:
所以如果通带衰减不是3dB,不能笼统的认为Ωp=Ωc,而是应该应用公式计算Ωc。在查表法将H(p)转换为H(s)时,对低通滤波器应该将p换为s/Ωc,对高通滤波器应该将p换成Ωc/s。
对Ωc的计算有两种方法,分别是利用通带截频、通带衰减系数和阻带截频、阻带衰减系数来计算Ωc的值。
由于N向上取整,所以两种方法的计算结果会略有不同。
原码对所有测试用例均正确的原因是,由于衰减系数和阶数对Ωc和Ωp的近似程度影响较小,对λs的影响也就较小。带入归一化频率公式:
对由于分母进行取对数运算,导致误差影响更小,以至于对N的影响几乎为0。
但是将H(p)转为H(s)时候就“露馅”了,因为对低通滤波器和高通滤波器来说,我们需要将p分别替换的是s/Ωc和Ωc/s。
所以综上,对于通带截频不是3dB的情况,我们还是应该首先计算Ωc的值。
如果不是通带截频不是3dB的话,我们就需要确定Ωc——>确定λs——>从而用归一化频率的N计算式来确定N,但是我们在确定Ωc的时候就需要使用N这个值,所以——
对更为严谨精确的运算过程,我们应该采用的是非归一化频率的N计算式,得到N之后,判断通带衰减系数是不是3dB,如果不是,我们计算Ωc,并最终应用于H(p)转H(s)的过程中。
所以我们需要改动的函数是complication()。
原码:
void compilation()
{
ld center = 0;
input();
change();
//判断滤波器类型
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:
throw"error:参数不符合滤波器特征"; break;
}
cout << endl << "参量下滤波器存在性与对应方法局限性自检成功…" << endl;
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
|| check_types(p1, p2, s1, s2) == 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 compilation()
{
ld center = 0;
input();
change();
//判断滤波器类型
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:
throw"error:参数不符合滤波器特征"; break;
}
cout << endl << "参量下滤波器存在性与对应方法局限性自检成功…" << endl;
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的值
double temp_1 = log10((pow(10, 0.1 * As) - 1) / (pow(10, 0.1 * Ap) - 1));
double temp_2 = log10(abs(s2 - s1) / abs(p2 - p1));
//cout << temp_1 << endl;
//cout << temp_2 << endl;
int N = ceil(temp_1 / (2 * temp_2));
//计算3dB截频
Acp = (Ap == 3) ? abs(p2 - p1) : (abs(p2 - p1) / pow((pow(10, 0.1 * Ap) - 1), 0.5 / N));
Acs = (Ap == 3) ? abs(s2 - s1) : (abs(s2 - s1) / pow((pow(10, 0.1 * As) - 1), 0.5 / N));
cout << endl << "用通带截频和衰减系数计算Ωc得:" << Acp << endl;
cout << endl << "用阻带截频和衰减系数计算Ωc得:" << Acs << endl;
//根据N的值查表得到低通滤波器归一化传输函数Hlp(p)分母表达式
cout << endl << "butterworth滤波器的阶数为N = " << N << ";" << endl;
cout << endl << "归一化传输函数为:H(p) = ";
show(N, "p");
cout << "用Ωp计算得到的";
string show_next_level_Acp = is_p(check_types(p1, p2, s1, s2), Acp, center);
string show_next_level_Acs = is_p(check_types(p1, p2, s1, s2), Acs, center);
cout << endl << "传输函数H(s) = ";
show(N, show_next_level_Acp);
cout << "用Ωs计算得到的";
cout << endl << "传输函数H(s) = ";
show(N, show_next_level_Acs);
}
其中,新增的变量Acp和Acs分别对应:用通带参量计算出的Ωc和用阻带参量计算出的Ωc。
二、总体代码
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
#define UN_KNOWN 5
#define USING_IIM 6
#define USING_BLT 7
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}
};
int filter_type;
int operate_type;
ld p1, p2, s1, s2, Ap, As, fs, Acp, Acs;
//标准输出
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;
}
//判断选通滤波器类型
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);
}
void change()
{
p1 *= 2 * PI; p2 *= 2 * PI; s1 *= 2 * PI; s2 *= 2 * PI;
//如果是双线性变换法,需要进行预畸变处理。
//如果是低通或高通,原本是0的数字频率经预处理后还是0不变,归一化了。
if (operate_type == USING_BLT)
{
p1 = 2 * fs * tan(0.5 * p1 / fs);
p2 = 2 * fs * tan(0.5 * p2 / fs);
s1 = 2 * fs * tan(0.5 * s1 / fs);
s2 = 2 * fs * tan(0.5 * s2 / fs);
}
//cout << p1 << endl << p2 << endl << s1 << endl << s2 << endl;
//如果是冲激响应不变法,则不用进行预畸变处理。
//由于冲激响应不变法不能处理高通和带阻滤波器,
//所以当程序判断为高通或带阻滤波器时,应报错。
filter_type = check_types(p1, p2, s1, s2);
if (operate_type == USING_IIM)
if (filter_type == IS_HP || filter_type == IS_BS)
throw"error:冲激响应不变法不能用于设计高通或带阻滤波器";
}
auto format_double_value(double val, int fixed) {
std::ostringstream oss;
oss << std::setprecision(fixed) << val;
return oss.str();
}
string is_p(int type, ld bindwith, ld& center)
{
string output;
if (type == IS_LP)
output = "(s/" + format_double_value(bindwith, define_setpercision) + ")";
else if (type == IS_HP)
output = "(" + format_double_value(bindwith, define_setpercision) + "/s";
else if (type == IS_BP)
output = "((s^2+" + format_double_value(pow(center, 2), define_setpercision) + ")/" +
"(" + format_double_value(abs(p2 - p1), define_setpercision) + "*s" + "))";
else if (type == IS_BS)
output = "((" + format_double_value(abs(p2 - p1), define_setpercision) + "*s" + ")/" +
"(s^2+" + format_double_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 input()
{
string get_input;
cout << endl << "欢迎登录IIR_DF设计系统!" << endl;
cout << endl << "我们需要获取您的一些数据要求~" << endl;
cout << endl << "请问您想使用哪种方法设计IID_DF?" << endl;
cout << "如果想使用冲激响应不变法,请输入iim;\
\n如果想使用双线性不变法,请输入blt:" << endl;
cin >> get_input;
(get_input == "iim" || get_input == "blt") ?
(operate_type = (get_input == "iim") ? USING_IIM : USING_BLT) :
(throw "error:输入错误!");
cout << endl << "下面我们将获取您希望最终得到的数字指标。" << endl;
cout << "请输入您想设计的数字滤波器类型(包括LP、HP、BP、BS、UNKNOWN):";
cin >> get_input;
if (get_input == "LP")filter_type = IS_LP;
else if (get_input == "HP")filter_type = IS_HP;
else if (get_input == "BP")filter_type = IS_BP;
else if (get_input == "BS")filter_type = IS_BS;
else if (get_input == "UNKNOWN") filter_type == UN_KNOWN;
else { cout << "error:输入错误!"; exit(0); }
cout << endl << "下面是关于键入频率和衰减系数的指标的一些声明。\n\
声明:\n\
因为在模拟域中,低通滤波器的通带和阻带下截止频率均趋于负无穷,\n\
所以在数字域中,低通滤波器的通带和阻带下截止频率均趋于0,\n\
所以如果是低通滤波器,程序默认在通带和阻带下截频位置(p1、s1)键入0;\n\
同理,\n\
因为在模拟域中,高通滤波器的通带和阻带上截止频率均趋于正无穷,\n\
所以在数字域中,高通滤波器的通带和阻带上截止频率均趋于0,\n\
所以如果是低通滤波器,程序默认在通带和阻带上截频位置(p2、s2)键入0;" << endl;
cout << endl << "请输入取样频率(Hz):"; cin >> fs;
cout << endl << "请输入通带起始频率p1(Hz):";
if (filter_type == IS_LP) { p1 = 0; cout << p1; }
else cin >> p1;
cout << endl << "请输入通带截止频率p2(Hz):";
if (filter_type == IS_HP) { p2 = 0; cout << p2; }
else cin >> p2;
cout << endl << "请输入通带衰减系数(dB):"; cin >> Ap;
cout << endl << "请输入阻带起始频率s1(Hz):";
if (filter_type == IS_LP) { s1 = 0; cout << s1; }
else cin >> s1;
cout << endl << "请输入阻带截止频率s2(Hz):";
if (filter_type == IS_HP) { s2 = 0; cout << s2; }
else cin >> s2;
cout << endl << "请输入阻带衰减系数(dB):"; cin >> As;
cout << endl;
cout << "运算结果为:" << endl;
}
void compilation()
{
ld center = 0;
input();
change();
//判断滤波器类型
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:
throw"error:参数不符合滤波器特征"; break;
}
cout << endl << "参量下滤波器存在性与对应方法局限性自检成功…" << endl;
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的值
double temp_1 = log10((pow(10, 0.1 * As) - 1) / (pow(10, 0.1 * Ap) - 1));
double temp_2 = log10(abs(s2 - s1) / abs(p2 - p1));
//cout << temp_1 << endl;
//cout << temp_2 << endl;
int N = ceil(temp_1 / (2 * temp_2));
//计算3dB截频
Acp = (Ap == 3) ? abs(p2 - p1) : (abs(p2 - p1) / pow((pow(10, 0.1 * Ap) - 1), 0.5 / N));
Acs = (Ap == 3) ? abs(s2 - s1) : (abs(s2 - s1) / pow((pow(10, 0.1 * As) - 1), 0.5 / N));
cout << endl << "用通带截频和衰减系数计算Ωc得:" << Acp << endl;
cout << endl << "用阻带截频和衰减系数计算Ωc得:" << Acs << endl;
//根据N的值查表得到低通滤波器归一化传输函数Hlp(p)分母表达式
cout << endl << "butterworth滤波器的阶数为N = " << N << ";" << endl;
cout << endl << "归一化传输函数为:H(p) = ";
show(N, "p");
cout << "用Ωp计算得到的";
string show_next_level_Acp = is_p(check_types(p1, p2, s1, s2), Acp, center);
string show_next_level_Acs = is_p(check_types(p1, p2, s1, s2), Acs, center);
cout << endl << "传输函数H(s) = ";
show(N, show_next_level_Acp);
cout << "用Ωs计算得到的";
cout << endl << "传输函数H(s) = ";
show(N, show_next_level_Acs);
}
int main()
{
system("color 0A");
try
{
compilation();
}
catch (const char* cc)
{
cout << cc;
exit(0);
}
return 0;
}
2.2执行效果 例题4-10
除了H(s),其余均与 北邮22级信通院DSP:IIR_DF系统2.0版:用C++程序自主写题——第四章IIR_DF设计的大部分题目——给定参量和选用的方法,设计出IIR_DF,一直进行到得出H(s)-CSDN博客
中的执行效果相同。