C++实现matlab的fir1函数

已于 2022-04-27 18:14:16 修改
阅读量3.7k 收藏 47
点赞数 14
分类专栏： matlab转C FIR滤波器 文章标签： matlab c++
于 2021-08-02 11:25:28 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_47898198/article/details/119316275
版权

C++实现matlab的fir1函数

导言

最近在进行Qt开发，涉及大量的matlab转C的工作，其中包括插值滤波等，但遗憾就求滤波系数的函数fir1而言，多数都是直接用的matlab生成的系数进行滤波，很少有用C生成的，有少数用C进行实现也和matlab生成的系数相差甚远，因此这里对matlab生成滤波系数的fir1函数进行了研究并用C++进行了重写，以备进行FIR滤波。
经对比，结果与matlab的fir1函数生成的滤波系数完全一致。

函数需求分析

这里我对工作中调用matlab的fir1()函数的实例为例子进行讲解：

  1. 原matlab程序 ：由256（fir1()函数里进行了+1）阶，截止频率为0.25的低通滤波器（默认使用hamming窗）
   lbn=255;
   lbf=fir1(lbn,0.25,'low');

  1. fir1（）生成的结果
    matlab生成结果图

  2. 分析
    可以看到，输入参数为阶数n，截止频率Wn，选择低通滤波（也是默认的选择），输出参数为一个一维数组，即滤波参数。

  3. C++版fir1（）输入输出参数详情：

参数名介绍
n阶数，对应matlab的fir1里的阶数n
Wn对应matlab的fir1里的阶数Wn，但应注意传进来的数据应存在一个vector的double数组里。
hdouble类型的vector数组，里面存的是滤波器系数
  vector <double> fir1(int numtaps, vector<double> cutoff)
  /*
  	未写排错  检查输入有需要自己进行完善
  	原matlab函数fir(n, wn)	【函数默认使用hamming】

  	参数输入介绍：
  		n：  对应matlab的fir1里的阶数n
  		Wn:  对应matlab的fir1里的阶数Wn，但应注意传进
  				 来的数据应存在一个vector的double数组里。

  	参数输出介绍：
  				vector <double>的一个数组，里面存的是长度
  				为n的滤波器系数。
  */;
  1. C++版fir1函数生成结果 ：结果表明，生成结果与matlab生成结果一致。
    在这里插入图片描述

数学过程

数学小辣鸡的我，希望得到大家的指正

在这里插入图片描述

源码

那么废话不多说，直接上代码

#include <stdio.h>
#include <math.h>
#include <vector>
#include <iostream>
using namespace std;
#define PI acos(-1)
vector<double> sinc(vector<double> x)
{
   vector<double> y;
   for (int i = 0; i < x.size(); i++)
   {
   	double temp = PI * x[i];
   	if (temp == 0) {
   		y.push_back(0.0);
   	}
   	else {
   		y.push_back(sin(temp) / temp);
   	}
   }
   return y;
}

vector <double> fir1(int n, vector<double> Wn)
{

   /*
   	未写排错  检查输入有需要自己进行完善
   	原matlab函数fir(n, wn)	【函数默认使用hamming】

   	参数输入介绍：
   		n：  对应matlab的fir1里的阶数n
   		Wn:  对应matlab的fir1里的阶数Wn，但应注意传进
   				 来的数据应存在一个vector的double数组里。

   	参数输出介绍：
   				vector <double>的一个数组，里面存的是长度
   				为n的滤波器系数。
   */
   
   //在截止点的左端插入0（在右边插入1也行）
   //使得截止点的长度为偶数，并且每一对截止点对应于通带。
   if (Wn.size() == 1 || Wn.size() % 2 != 0) {
   	Wn.insert(Wn.begin(), 0.0);
   }

   /*
   	‘ bands’是一个二维数组，每行给出一个 passband 的左右边缘。
   	（即每2个元素组成一个区间）
   */
   vector<vector <double>> bands;
   for (int i = 0; i < Wn.size();) {
   	vector<double> temp = { Wn[i], Wn[i + 1] };
   	bands.push_back(temp);
   	i = i + 2;
   }

   // 建立系数
   /*
   	m = [0-(n-1)/2,
   		 1-(n-1)/2,
   		 2-(n-1)/2,
   		 ......
   		 255-(n-1)/2]
   	h = [0,0,0......,0]
   */
   double alpha = 0.5 * (n - 1);
   vector<double> m;
   vector<double> h;
   for (int i = 0; i < n; i++) {
   	m.push_back(i - alpha);
   	h.push_back(0);
   }
   /*
   	对于一组区间的h计算
   	left:	一组区间的左边界
   	right:  一组区间的右边界
   */
   for (int i = 0; i < Wn.size();) {
   	double left = Wn[i];
   	double right = Wn[i+1];
   	vector<double> R_sin, L_sin;
   	for (int j = 0; j < m.size(); j++) {
   		R_sin.push_back(right * m[j]);
   		L_sin.push_back(left * m[j]);
   	}
   	for (int j = 0; j < R_sin.size(); j++) {
   		h[j] += right * sinc(R_sin)[j];
   		h[j] -= left * sinc(L_sin)[j];
   	}

   	i = i + 2;
   }

   // 应用窗口函数，这里和matlab一样
   // 默认使用hamming，要用别的窗可以去matlab查对应窗的公式。
   vector <double> Win;
   for (int i = 0; i < n; i++)
   {
   	Win.push_back(0.54 - 0.46*cos(2.0 * PI * i / (n - 1)));	//hamming窗系数计算公式
   	h[i] *= Win[i];
   }

   bool scale = TRUE;
   // 如果需要，现在可以处理缩放.
   if (scale) {
   	double left = bands[0][0];
   	double right = bands[0][1];
   	double scale_frequency = 0.0;
   	if (left == 0)
   		scale_frequency = 0.0;
   	else if (right == 1)
   		scale_frequency = 1.0;
   	else
   		scale_frequency = 0.5 * (left + right);

   	vector<double> c;
   	for (int i = 0; i < m.size(); i++) {
   		c.push_back(cos(PI * m[i] * scale_frequency));
   	}
   	double s = 0.0;
   	for (int i = 0; i < h.size(); i++) {
   		s += h[i] * c[i];
   	}
   	for (int i = 0; i < h.size(); i++) {
   		h[i] /= s;
   	}
   }
   return h;
}

使用

//设置低通滤波器
   int lbn = 255;
   vector<double> Wn;
   Wn.push_back(0.25);
   vector<double> lbf = fir1(lbn + 1, Wn);

计算量大时性能提升（vector换成了动态数组）

   #include <stdio.h>
   #include <math.h>
   #include <vector>
   #include <iostream>
   using namespace std;
   #define PI acos(-1)
   
   double sincEasy(double *x, int len, int index) {
       double temp = PI * x[index];
       double y;
       if (temp == 0) {
           y = 0.0;
       }
       else {
           y = sin(temp) / temp;
       }
       return y;
   }

   //滤波器系数
   double *fir1(
           int lbflen,
           double Wn[],
           double lbf[])
   {

       /*
           未写排错  检查输入有需要自己进行完善
           原matlab函数fir(j, wn)	【函数默认使用hamming】

           参数输入介绍：
               j：  对应matlab的fir1里的阶数j
               Wn:  对应matlab的fir1里的阶数Wn，但应注意传进
                        来的数据应存在一个vector的double数组里。

           参数输出介绍：
                       vector <double>的一个数组，里面存的是长度
                       为j的滤波器系数。
       */

       //在截止点的左端插入0（在右边插入1也行）
       //使得截止点的长度为偶数，并且每一对截止点对应于通带。
       //if (Wn.size() == 1 || Wn.size() % 2 != 0) {
       //	Wn.insert(Wn.begin(), 0.0);
       //}

       double alpha = 0.5 * (lbflen - 1);
       double *m = new double[lbflen];
       for (int i = 0; i < lbflen; i++) {
           m[i] = i - alpha;
           lbf[i] = 0;
       }

       double *R_sin = new double[lbflen];
       double *L_sin = new double[lbflen];
       for (int i = 0; i < 2;) {
           double left = Wn[i];
           double right = Wn[i + 1];
           for (int j = 0; j < lbflen; j++) {
               R_sin[j] = right * m[j];
               L_sin[j] = left * m[j];
           }
           for (int j = 0; j < lbflen; j++) {
               lbf[j] += right * sincEasy(R_sin, lbflen, j);
               lbf[j] -= left * sincEasy(L_sin, lbflen, j);
           }

           i = i + 2;
       }

       // 应用窗口函数，这里和matlab一样
       // 默认使用hamming，要用别的窗可以去matlab查对应窗的公式。
       for (int i = 0; i < lbflen; i++)
       {
           double Win = 0.54 - 0.46*cos(2.0 * PI * i / (lbflen - 1));	//hamming窗系数计算公式
           lbf[i] *= Win;
       }

       // 如果需要，现在可以处理缩放.
       if (TRUE) {
           double left = Wn[0];
           double right = Wn[1];
           double scale_frequency = 0.0;
           if (left == 0)
               scale_frequency = 0.0;
           else if (right == 1)
               scale_frequency = 1.0;
           else
               scale_frequency = 0.5 * (left + right);

           double s = 0.0;
           for (int i = 0; i < lbflen; i++) {
               double c = cos(PI * m[i] * scale_frequency);
               s += lbf[i] * c;
           }
           for (int i = 0; i < lbflen; i++) {
               lbf[i] /= s;
           }
       }
       delete[] m;
       delete[] R_sin, L_sin;
       return lbf;
   }

使用

//设置低通滤波器
   int lbn = 255;
   double Wn[2] = { 0.0, 0.25 };
   double *lbf = new double[lbn + 1];
   lbf = fir1(lbn + 1, Wn, lbf);

注意

该文章仅个人学习使用，欢迎大家一起交流学习

宇智波盆
关注
专栏目录
matlab中窗函数的C++实现,以及fir1函数C++实现 经测试与matlab结果基本一致
09-09
窗函数包括：矩形窗(boxcar),三角窗(triang),图基窗(tukeywin),汉宁窗(hanning),海明窗(hamming)， 布拉克曼窗(blackman),凯塞窗(kaiser),切比雪夫窗(chebwin),高斯窗(gausswin),巴特里特窗(bartlett)。 fir1函数的实现
MATLAB函数fir1的C语言移值
weixin_39458727的博客
05-07 444
结果是完全一样的，OK！
【MATLAB】FIR滤波器的MATLAB实现
最新发布
aixdm的博客
09-06 2459
与IIR滤波器相比，FIR滤波器既有其优势也有其局限性。
C语言实现MATLAB中的fir1函数（含代码）
yangbupang的博客
08-23 5339
最近导师要求我将师兄写的一篇关于fsk调制解调的MATLAB代码用VS2013编成C语言，方便后续应用到工程实践中。由于本人编程水平有限，并且C语言只有在大学的时候水水的学过一遍，所以此次代码可能有许多疏漏，欢迎批评指正。 1.明确目标 本次要用C语言实现如下代码功能： recv_lowfilter = fir1(256,1.2*2*Br/fs); 这是一句MATLAB代码，该句代码的重点在于 ...
MATLAB函数fir1
热门推荐
芒果爱火锅的博客
05-03 9万+
fir1 基于窗函数的fir滤波器设计； 主要形式： b = fir1(n,Wn) b = fir1(n,Wn,ftype) b = fir1(___,window) b = fir1(___,scaleopt) 系数的含义 n — 滤波器阶数。 Wn — 截止频率， 0≤Wn≤1， Wn=1 对应于采样频率的一半。 当设计带通和带阻滤波器时， Wn=[W1W2],W1≤...
MATLAB函数——fir1
Josh_Gao
06-25 3万+
fir1：基于窗函数法的 FIR 滤波器设计
matlab设计FIR滤波器
hou19975418的博客
05-11 8482
B = fir1(N,Wn)设计FIR低通滤波器，返回的滤波器参数保存在长度为N+1的数组B中。%滤波函数，参数1为滤波器系数，参数2为行数（这里没有使用矩阵）故为1，参数3为输入信号。%N阶带通滤波器 B=fir1(N,Wn,'stop');%表示N阶高通滤波器 B=fir1(N,Wn,'low');Run一下生成的滤波器函数.m文件，会生成滤波器系数矩阵Hd，然后就可以调用filter函数对输入信号进行滤波了。
FIR1_fir_fir1c语言_FIR1_C语言实现fir1_matlab
09-11
对标Matlab中fir1函数的C语言实现
MATLAB中FIR滤波器的时延溢出问题详解：线性相位对信号造成的时延溢出及其消除方法，以及fir1等函数的使用
yhy0428的博客
03-06 1万+
用MATLAB进行滤波时，由于FIR滤波器线性相位特性造成的时延，信号会整体向后移动，造成信号序列溢出，本文将结合FIR滤波器性质及MATLAB各函数性质详细分析这一问题，并给出解决方案。
用C++实现 Matlab 中的Filtfilt
01-19
在本文中，我们将深入探讨如何使用C++来实现Matlab中的`Filtfilt`函数，这是一个用于零相位数字滤波的工具。`Filtfilt`的独特之处在于它通过两次滤波过程，一次正向，一次反向，来实现零相位滤波，从而在滤波过程中...
matlab filtfilt函数的c实现
11-02
在MATLAB中，`filtfilt`函数是基于FIR（Finite Impulse Response）或IIR（Infinite Impulse Response）滤波器的，它通过在正向和反向应用滤波器来实现零相位滤波。这种零相位特性使得滤波后信号的时间位置不会发生...
matlab smooth函数C++实现,matlab的smooth函数,C,C++源码.zip
09-30
然而，MATLAB的`smooth`函数在C++环境中并不能直接使用，需要程序员自己实现。 在C++中，我们可以采用类似的方法来实现数据平滑。以下是一些关键的知识点： 1. **移动平均**：移动平均是最简单的平滑方法之一，它...
matlab中fir1c语言,C语言实现fir1函数?
weixin_35940165的博客
05-19 621
#include #ifdef WIN32#include #endif#define SAMPLE double /* define the type used for data samples */void clear(int ntaps, SAMPLE z[]){int ii;for (ii = 0; ii < ntaps; ii++) {z[ii] = 0;}}SAMPLE fir_...
matlab 与C语言关于filter函数的实现
enjoy_learn的博客
12-05 7447
MATLAB实现滤波器，滤波器类似：低通、高通、带通、带阻。 实现方式：（1）利用Matlab代码生成DLL供其他模块调用方法。 优点：实现速度快，可利用Ｃ直接调用matlab的DLL快速实现算法功能，无需用C考虑算法实现过程。 缺点：运行需要电脑上具备matlab运行环境MCR（或电脑上具有完整的matlab版本存在），且配置调试麻烦。 如翻译bpfilter滤波器，C调用matlabD
信号处理常用matlab之数字滤波器及滤波函数
weixin_37315722的博客
06-13 1万+
1 freqz 2 fir1 3 fir2 4
Matlab中filter函数用法
jk_101的博客
01-01 6万+
filter函数的功能是对数据进行1 维数字滤波。
C语言实现MATLAB中的fir1函数(绝对一摸一样的系数矩阵)
thecentury的博客
10-29 1331
#include<iostream> #include<vector> #include<iomanip> #define _USE_MATH_DEFINES #include<math.h> #include<cmath> //#include<NumCpp.hpp> #define pi M_PI using std::vector; void create_hidea(int n, double fc, vector<d
FIR数字滤波器 C语言频域实现 matlab设计FIR参数
zhenhui1314的博客
04-11 3412
FIR数字滤波器 C语言频域实现 matlab设计FIR参数一、FIR实现步骤1、FIR参数做FFT2、输入数据做FFT3、FIR做复数乘积4、将复数乘积结果做IFFT输出二、FIR代码实现1、FIR参数的生成2、FIR系数做FFT计算3、输入数据做overlap4、输入数据做FFT5、数据做复数乘积6、数据做IFFT并输出三、代码说明以及待完善1、代码说明2、待完善 本文主要介绍用C语言实现频...
matlab的fir1的ftype,数字滤波器的matlab与FPGA实现读书笔记（二）FIR滤波器的matlab设计（窗函数）...
weixin_39717443的博客
03-19 816
FIR滤波器的matlab设计采样fir1函数设计b= fir1(n,wn),b= fir1(n,wn,’ftype’)b= fir1(n,wn,’ftype’,window)b= fir1(…,’noscale’)函数参数说明b:返回的FIR滤波器单位脉冲响应，脉冲响应为偶对称，长度为n+1n:滤波器的阶数，需要注意的是，设计出的滤波器长度为n+1.wn:滤波器的截至频率，；需要注意的是，wn的...
评论 15
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

打赏作者

宇智波盆

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

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

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

打赏作者

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

抵扣说明：

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

余额充值