C/C++知识点之IIR滤波器软件实现(Matlab+C++)
小标 2018-12-26 来源 : 阅读 1668 评论 0

摘要:本文主要向大家介绍了 C/C++知识点之IIR滤波器软件实现(Matlab+C++),通过具体的内容向大家展示,希望对大家学习C/C++知识点有所帮助。

本文主要向大家介绍了 C/C++知识点之IIR滤波器软件实现(Matlab+C++),通过具体的内容向大家展示,希望对大家学习C/C++知识点有所帮助。

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

本文由职坐标整理并发布,希望对同学们有所帮助。了解更多详情请关注职坐标编程语言C/C+频道!

本文由 @小标 发布于职坐标。未经许可,禁止转载。
喜欢 | 1 不喜欢 | 0
看完这篇文章有何感觉?已经有1人表态,100%的人喜欢 快给朋友分享吧~
评论(0)
后参与评论

您输入的评论内容中包含违禁敏感词

我知道了

助您圆梦职场 匹配合适岗位
验证码手机号,获得海同独家IT培训资料
选择就业方向:
人工智能物联网
大数据开发/分析
人工智能Python
Java全栈开发
WEB前端+H5

请输入正确的手机号码

请输入正确的验证码

获取验证码

您今天的短信下发次数太多了,明天再试试吧!

提交

我们会在第一时间安排职业规划师联系您!

您也可以联系我们的职业规划师咨询:

小职老师的微信号:z_zhizuobiao
小职老师的微信号:z_zhizuobiao

版权所有 职坐标-一站式IT培训就业服务领导者 沪ICP备13042190号-4
上海海同信息科技有限公司 Copyright ©2015 www.zhizuobiao.com,All Rights Reserved.
 沪公网安备 31011502005948号    

©2015 www.zhizuobiao.com All Rights Reserved

208小时内训课程