北邮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 7typedef 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);elsereturn (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) + "))";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 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博客
中的执行效果相同。