• 日期: 2024-07-31
  • 标签: 常规
#include 'stdafx.h'
#include 'D:\xhclgcyy\x_math.cpp'
#include 'D:\xhclgcyy\x_graph.cpp'


//变量声明
int order, L, N = 256;//滤波器阶数、级联个数、信号抽样点数
double b[10][2], c[2][3], H[10][2][5], a[1000] = { 0 }, b0[1000] = { 0 };
COMPLEX y[1000], x[1000], X[1000], Y[1000];


//二维坐标图工具
void plotgri2(COLORREF gridcolor, COLORREF linecolor, COMPLEX p[], int N)
{
	int i;
	HPEN pen1 = CreatePen(PS_SOLID, 1, gridcolor), oldpen = (HPEN)SelectObject(win3.hdc, pen1);
	HPEN pen2 = CreatePen(PS_SOLID, 1, linecolor);
	for (i = 0; i<N; i++)
	{
		moveto2(i, 0);
		lineto2(i, p[i].r);
	}
}


//相频曲线
double P(double f, double fs, double H[][2][5], int L)
{
	double B;
	int i;
	COMPLEX m;
	COMPLEX z_1(cos(8 * f*(1.0) / fs), -sin(8 * f*(1.0) / fs));
	for (B = 0, i = 0; i<L; i++)
	{
		m = polyval(&(H[i][0][0]), 5, z_1) / polyval(&(H[i][1][0]), 5, z_1);
		B = B + (atan2(m.i, m.r));
	}
	return B;
}


//级联数
int GET_L(double b[][2], int order)
{
	int k, M, L;
	M = order / 2; L = M;
	if (order != M * 2) { b[M][0] = 0; b[M][1] = 1; L++; }
	for (k = 0; k<M; k++)
	{
		b[k][0] = 1;
		b[k][1] = (-2)*cos((2 * k + order + 1) * 4 * atan(1.0) / (2 * order));
	}
	return L;
}


//滤波类型1234分别为LP、HP、BP、BS;a1,a2,fs,f1,f2,f3,f4
void IIR(int bandType, double a1, double a2, double fs, double f1, double f2, double f3, double f4)
{
	int i, j;
	order = btwOrder(bandType, a1, a2, fs, f1, f2, f3, f4);//阶数
	L = GET_L(b, order);//级连数
	btwC23(c, bandType, order, a1, fs, f1, f2, f3, f4);
	btwAf2Df(H, L, b, c);//模拟滤波转数字滤波

	printf('级联级数L=%d\n', L);
	for (i = 0; i<L; i++)//用printf()格式化把系统函数的系数H数组和级连数L显示出来
	{
		for (j = 0; j<2; j++)
		{
			if (j == 0)
			{
				printf('      %f+%f*z^(-1)+%f*z^(-2)+%f*z^(-3)+%f*z^(-4)', H[i][j][0], H[i][j][1], H[i][j][2], H[i][j][3], H[i][j][4]);
			}
			else
			{
				printf(' H(%d)=-------------------------------------------------------------------------', i + 1);
				printf('     %f+%f*z^(-1)+%f*z^(-2)+%f*z^(-3)+%f*z^(-4)\n\n', H[i][j][0], H[i][j][1], H[i][j][2], H[i][j][3]), H[i][j][4];
			}
		}
	}
	_getch();


	window2(L'幅频', -1., 10., 1000., -100., 'hz', 'db');
	xy2(BLUE);
	plotxy2(RED, 2, f, btw20lgHz(f,fs,H,L));


	//绘制通阻带参考线：
	if (bandType == LOWPASS || bandType == HIGHPASS)
	{
		line2(0, -a1, win2.x2, -a1);
		line2(0, -a2, win2.x2, -a2);
		line2(f1, 0, f1, win2.y1);
		line2(f2, 0, f2, win2.y1);
	}
	else {
		line2(0, -a1, win2.x2, -a1);
		line2(0, -a2, win2.x2, -a2);
		line2(f1, 0, f1, win2.y1);
		line2(f2, 0, f2, win2.y1);
		line2(f3, 0, f3, win2.y1);
		line2(f4, 0, f4, win2.y1);
	}
	_getch();
	frame2();

	window2(L'相频', -600, -30, 600., 30., 'hz', 'db');
	xy2(BLUE);
	plotxy2(RED, 2, f, P(f, fs, H, L));
	_getch();
	frame2();
}


//方波采样 
void chouyang(double fc)
{
	int i;
	double F;
	double t = 0;
	F = 2 * fc;//方波抽样 T=2
	for (i = 0; i<F; i++) { y[i] = 0; }
	for (i = 0; i<F; i++)
	{
		t = i * 1.0 / (1.0*fc);
		x[i] = COMPLEX((t >= 0 && t <= 1 ? 1 : 0.0), 0);//一个周期的方波
	}
}


//通用IIR滤波器 
void PUBLIC_IIR(COMPLEX input[], COMPLEX output[], double a[], double b[], int N, int Nmax)
{//N为数据点数，Nmax为0状态最大负输入
	int i, n;
	for (n = 0; n<N; n++)
	{
		for (i = 1; i<Nmax; i++)
		{
			output[n] = (n<i ? output[n] + a[i] * 0 : output[n] + a[i] * output[(n - i)]);
		}
		for (i = 0; i<Nmax; i++)
		{
			output[n] = (n<i ? output[n] + b[i] * 0 : output[n] + b[i] * input[(n - i)]);
		}
	}
}


//低通L级级联
void LP_JL(void)
{
	int i, j, n;
	window2(L'抽样图形显示', -1, -1, 200, 2, 'i', 'x[i]');
	xy2(RED);
	plotgri2(BLUE, RED, x, N);
	_getch();
	frame2();
	dft(X, x, N, 1);
	for (i = 0; i<N; i++)
		Y[i] = abs(X[i]);

	window2(L'DFT后频谱显示', -1, -1, 300, 30, 'i', 'X[i]');
	xy2(RED);
	plotgri2(BLUE, RED, Y, N);
	_getch();
	frame2();
	for (i = 0; i <= L - 1; i++)
	{
		a[0] = 0;
		for (j = 1; j<5; j++)
		{
			a[j] = -H[i][1][j] / H[i][1][0];
		}//归一化 x
		for (j = 0; j<5; j++)
		{
			b0[j] = H[i][0][j] / H[i][1][0];
		}//归一化 y

		PUBLIC_IIR(x, y, a, b0, N, 5);
		for (n = 0; n<N; n++)
			x[n] = y[n];
	}
	for (n = 0; n<N; n++)
		y[n] = x[n];

	dft(X, y, N, 1);
	for (i = 0; i<N; i++)
		Y[i] = abs(X[i]);

	window2(L'L级联后图形显示', -1, -200, 300, 1500, 'i', 'y[i]');
	xy2(RED);
	plotgri2(BLUE, RED, y, N);
	_getch();
	frame2();
	window2(L'函数图形显示', -1, -1, 300, 20000, 'i', 'X[i]');
	xy2(RED);//画xy轴。
	plotgri2(BLUE, RED, Y, N);
	_getch();
	frame2();
}


//对方波进行低通滤波
void LP_IIR(int bandType, double a1, double a2, double fs, double f1, double f2, double f3, double f4)
{
	order = btwOrder(bandType, a1, a2, fs, f1, f2, f3, f4);//阶数
	L = GET_L(b, order);//级连数
	btwC23(c, bandType, order, a1, fs, f1, f2, f3, f4);
	btwAf2Df(H, L, b, c);//模拟滤波转数字滤波
	chouyang(100.);//方波采样
	LP_JL();
	_getch();
}


void main()
{
	//IIR滤波器（参数:滤波类型1234分别为LP、HP、BP、BS;a1,a2,fs,f1,f2,f3,f4)
	IIR(1,3.0,35.0,2000.,200.,300.,400.,500.);//设计低通数字滤器fs=2000
	IIR(2,3.0,35.0,2000.,200.,300.,400.,500.);//设计高通数字滤器fs=2000
	IIR(3,3.0,40.0,2000.,200.,300.,400.,500.);//设计带通数字滤器fs=2000
	IIR(4,3.0,40.0,2000.,200.,300.,400.,500.);//设计带阻数字滤器fs=2000
	LP_IIR(1,3.0,35.0,2300.,200.,300.,400.,500.);//对方波进行低通滤波;fs=2300
}

代码分析

代码实现了多种IIR数字滤波器的设计与实现，并对方波信号进行滤波处理。

函数功能

  1. btwOrder函数：根据滤波类型和设计参数计算IIR滤波器的阶数。
  2. btwC23函数：根据滤波类型和设计参数计算C23滤波器的系数。
  3. btwAf2Df函数：将模拟滤波器的系数转换为数字滤波器的系数。
  4. P函数：计算相频曲线。
  5. GET_L函数：计算级联数。
  6. IIR函数：设计IIR滤波器，并绘制幅频曲线和相频曲线。
  7. chouyang函数：对方波进行采样。
  8. PUBLIC_IIR函数：通用IIR滤波器函数。
  9. LP_JL函数：低通L级级联函数。
  10. LP_IIR函数：对方波进行低通滤波函数，调用了LP_JL函数。

代码亮点

  • 使用了C++语言进行编写，代码结构清晰易懂。
  • 实现了多种类型的IIR滤波器，包括低通、高通、带通和带阻滤波器。
  • 提供了完整的代码示例，方便读者学习和使用。
  • 对代码进行了详细的注释，方便读者理解代码的功能和实现细节。
  • 使用了图表展示了滤波器的幅频特性和相频特性，更加直观地展示了滤波器的性能。

总结

本文提供了一个基于C++的IIR数字滤波器设计与实现的完整示例，并对其代码进行了详细的分析。读者可以通过学习本文的代码示例，快速掌握IIR数字滤波器的设计与实现方法。

基于C++的IIR数字滤波器设计与实现

原文地址: https://www.cveoy.top/t/topic/jQY5 著作权归作者所有。请勿转载和采集!

免费AI点我,无需注册和登录