北邮22级信通院DSP：IIR_DF系统3.0版：从H（p）到H（s）：一种更为严谨精确的运算模式-CSDN博客

本文链接：https://blog.csdn.net/bc202205/article/details/139365494

北邮22信通一枚~

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

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

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

承接上一篇博客

北邮22级信通院DSP：IIR_DF系统2.0版：用C++程序自主写题——第四章IIR_DF设计的大部分题目——给定参量和选用的方法，设计出IIR_DF，一直进行到得出H（s）-CSDN博客

在原有的基础上，应对通带衰减系数不是3dB的情况，提出了一种更为严谨精确的计算模式。

先来展示运行效果

一、更为严谨精确的complication（）函数

因为对butterworth滤波器来说，表达式为：

$\frac{1}{1+(\frac{\Omega }{\Omega_{c}})^{2N}}$

所以如果用butterworth滤波器进行拟合，就需要确定参量N和Ωc，其中Ωc的意思是通带3dB截止频率。没有引入频率归一化定义时，对N的计算公式应为：

$N\geqslant \frac{lg(\frac{10^{0.1As}-1}{10^{0.1Ap}-1})}{2lg(\frac{\Omega_{s}}{\Omega_{p}})}$

频率归一化的工作，主要是为了“查表法”操作的可实现性。λs定义为：

$\lambda _{s}=\frac{\Omega_{s}}{\Omega_{c}}$

所以如果通带衰减不是3dB，不能笼统的认为Ωp=Ωc，而是应该应用公式计算Ωc。在查表法将H（p）转换为H（s）时，对低通滤波器应该将p换为s/Ωc，对高通滤波器应该将p换成Ωc/s。

对Ωc的计算有两种方法，分别是利用通带截频、通带衰减系数和阻带截频、阻带衰减系数来计算Ωc的值。

$\Omega _{c}=\frac{\Omega _{p}}{\sqrt[2N]{10^{0.1Ap}-1}}$

$\Omega _{c}=\frac{\Omega _{s}}{\sqrt[2N]{10^{0.1As}-1}}$

由于N向上取整，所以两种方法的计算结果会略有不同。

原码对所有测试用例均正确的原因是，由于衰减系数和阶数对Ωc和Ωp的近似程度影响较小，对λs的影响也就较小。带入归一化频率公式：

$N\geqslant \frac{lg(\frac{10^{0.1As}-1}{10^{0.1Ap}-1})}{2lg(\frac{\lambda _{s}}{\lambda _{p}})}$

对由于分母进行取对数运算，导致误差影响更小，以至于对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;
}