使用 C++ 来写一个 IIR 滤波器
我们首先要在 MATLAB 中设计一个 IIR 滤波器, 并生成一个头文件, 这个头文件中反映了 IIR 滤波器的频率响应特性
理论支持
IIR 滤波叫做递归滤波器, 它是一种具有反馈的滤波器. 当阶数较大时一般采取多个二阶节滤波进行串联, 这样可以提高系统稳定性.
一个二阶节系数规律如图所示:
可以写出第 K 个二阶节的差分方程
N 个二阶节的级联结构如下图所示:
根据二阶节图, 把前一级的输出作为后一级的输入, 就可以通过软件实现 IIR 数字滤波的功能.
使用 Matlab 生成头文件
首先打开 MATLAB 中 Filter Design & Analysis Tool
这里我们先设计一个低通滤波器
Fs 代表采样频率, 采样频率必须大于原信号最高频率的两倍,
否则会产生频谱混叠.
Fpass 为通带频率, Fstop 为阻带截止频率
这些参数设置好就可以点击 Design Filter
生成的是一个二阶节滤波组合, 一共有 31 阶, 也就是多个二阶滤波器的组合
接着在 Target 选项中生成 C Header File
Numerator 为分子系数数组的命名, Numerator length 为分子系数数组的长度,
Denominator 为分母.
对生成头文件进行分析
以下以 Fpass 为 10K,Fstop 为 12K 的低通滤波器举栵
在使用头文件前需要根据情况将 Matlab 的数据类型转换为 C++ 支持的数据类型, 这里我们使用 double 类型
在分析头文件前先看下 Matlab 提供的第一节滤波参数
以第一个二阶节的数据举例:
- Numerator: 1 2 1
- Denominator: 1 -0.55930961405944157 0.92579835996619642
- Gain:0.34162218647668868
Numerator 为分子的系数, 分别为 b0,b1,b2
Denominator 为分母的系数, 分别为 a0,a1,a2.
Gain 为各节的增益, 此项为了稳定各节, 稳定信号大小
接着对照头文件, 以下为头文件主要部分的一段截取:
- #define MWSPT_NSEC 41
- int NL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,
- 3,1,3,1,2,1 };
- double NUM[MWSPT_NSEC][3] = {
- {
- 0.3416221864767, 0, 0
- },
- {
- 1, 2, 1
- },
- {
- 0.3180955154747, 0, 0
- },
- {
- 1, 2, 1
- },......
MWSPT_NSEC 为滤波器阶数, 具体的节数在头文件开头的注释中
NL[MWSPT_NSEC]这个数组定义了 NUM[MWSPT_NSEC][3]数组每一行的有用数据个数(可以不用)
在 NUM[MWSPT_NSEC][3]数组 (分子参数) 奇数行第一项都为增益项, 偶数行为 3 个系数, 分别为 b0,b1,b2.
由此可以找出规律, 定义 K 为目前所在的阶数, p 为数组的首指针, 则, 每一节的增益项为(p+6*K), 第一个系数为(p+3+6*K),
第二个系数为(p+3+6*K+1), 第三个系数为(p+3+6*K+2).
C++ 编程实现
在软件设计的过程中, 每个二阶节的延迟变量只取 和 , 作为中间变量在过程中直接赋给 . 这是因为对于下一个输入数据 n+1 的延迟变量即为上一个输入数据的 和 , 采用这种方式进行设计, 可以节省寄存器的空间.
为了提高处理速度, 程序中需要使用指针进行参数传递, 特别注意二维数组的首地址传递方式为 & a[0][0]->double* a.
滤波函数
double iir(double *a, double *b,double* w, double xin, int N_IIR)
- {
- int k;
- double temp = xin;
- for (k = 0; k<N_IIR; k++)
- {
- *(w+k*3) = temp - *(a + 3+6 * k + 1) *(*(w + k * 3+1)) - *(a + 3 + 6 * k + 2) *(*(w + k * 3+2));
- // 这里 temp 为本二阶节的输入, 也是上一个二阶节的输出
- temp = *(b + 3 + 6 * k )* (*(w + k * 3)) + *(b + 3 + 6 * k + 1) * (*(w + k * 3+1)) + *(b + 3+6 * k + 2)* (*(w + k * 3+2));
- // 这里 temp 为本二阶节的输出, 也是下一个二阶节的输入
- *(w + k * 3 + 2) = *(w + k * 3 + 1);
- *(w + k * 3 + 1) = *(w + k * 3);
- temp = temp*(*(b + 6 * k));// 放大倍数, 稳定信号
- }
- return temp;
- }
实际测试
测试 Fpass 为 10K,Fstop 为 12K 的低通滤波器
在程序中输入三个频率为 2K,11K,20K 的信号, 理论上 2k 完全通过, 11k 部分衰减, 20K 完全滤除.
上图为原信号, 下图为滤波后信号.
实际测试发现符合设计要求, 而且在过渡带信号也基本完全衰减.
测试用 C++ 程序
- void main()
- {
- const int N = 100;
- int i,j;
- double xn[N];
- double w[20][3];
- double yn[N];
- for (i = 0; i < 20; i++)// 初始化
- {
- for (j = 0; j < 3; j++)
- w[i][j] = 0;
- }
- for (i = 0; i < N; i++)
- {
- xn[i] = sin(2 * 3.1416 * 20 / 50 * i)+ sin(2 * 3.1416 * 2 / 50 * i)+ sin(2 * 3.1416 * 11 / 50 * i);
- yn[i]=iir(&DEN[0][0], &NUM[0][0], &w[0][0],xn[i], 20);
- }
- ofstream SaveFile_a("xn.txt");
- for (i = 0; i<N; i++)
- SaveFile_a << " " << xn[i] << endl;
- SaveFile_a.close();
- ofstream SaveFile_b("yn.txt");
- for (i = 0; i<N; i++)
- SaveFile_b << " " << yn[i] << endl;
- SaveFile_a.close();
- }
分析用 Matlab 程序
- xn1=fopen('xn.txt','r');
- [xn,count]=fscanf(xn1,'%f');
- fclose(xn1);
N = length(xn);% 求取抽样点数
xn_f = fft(xn);% 对信号进行傅里叶变换
- xn_f=abs(xn_f(1:N/2));
- f = 50000/N*(0:N/2-1);
- subplot(211);
- stem(f,abs(xn_f));
- xlabel('Frequency / (s)');ylabel('Amplitude');
- title('原信号频谱');
- grid;
- yn1=fopen('yn.txt','r');
- [yn,count]=fscanf(yn1,'%f');
- fclose(yn1);
yn_f = fft(yn);% 对信号进行傅里叶变换
- yn_f=abs(yn_f(1:N/2));
- subplot(212);
- stem(f,abs(yn_f));
- xlabel('Frequency / (s)');ylabel('Amplitude');
- title('滤波后信号频谱');
- grid;
来源: https://www.cnblogs.com/geek-wireless/p/IIR.html