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

已于 2024-05-29 18:57:09 修改
阅读量927 收藏 26
点赞数 15
分类专栏： 北邮22级信通院DSP 文章标签： c++ dsp开发 IIR 北邮 信息与通信
于 2024-05-25 20:58:39 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/bc202205/article/details/139202586
版权

北邮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

三、更改日志

2024-05-29 09:10:03  

一、 核心算法

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 
		|| 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);
}

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 
		|| 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 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;

 

三、更改日志

2024-05-29 09:10:03  

complication（）函数中对lambda_s的判断赋值有误。

原码：

	ld lambda_s = (check_types(p1, p2, s1, s2) == (IS_BP ||I S_LP) ? 
		((s2 - s1) / (p2 - p1)) : ((p2 - p1) / (s2 - s1));

改正后：

	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));
青山入墨雨如画
关注
  • 15
    点赞
  • 踩
  • 26
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
专栏目录
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

青山入墨雨如画

你的鼓励将是我创作的动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值