matlab函数转C++(数字信号处理)

已于 2023-03-14 17:35:11 修改
阅读量1.4k 收藏 30
点赞数 10
分类专栏： matlab转C++ 文章标签： c++ 信号处理 matlab 矩阵
于 2023-02-16 00:35:26 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_45737966/article/details/129050338
版权

文章目录

前言

近期主要利用QT完成一个本科的通信教学软件，其中涉及大量matlab转C++的工作，本来是想利用matlab的Coder模块进行转换的，本人小白不太会用，还是自己按着matlab内置函数的代码进行转换，函数写的比较笨，希望大家能够多多指导.

matlab函数转C++

使用的是C++的armadillo矩阵库进行矩阵的运算，armadillo矩阵库内置许多信号处理算法，包括fft和ifft等运算等，但是一些matlab内置的函数还是没有的，这需要自己编写。

1.matlab的findpeaks函数
需求是找出一维矩阵的满足条件的谱峰数量，对应matlab的[fudu,~]=findpeaks(mat,'minpeakheight',1,'minpeakdistance',2);

整个算法主要包括两个部分：
1.1 计算二阶差分，得到峰值索引
1.2 将峰值与参数minpeakheight进行比对
1.3 计算筛选后的峰值之间的距离与minpeakdistance进行比对，得出满足条件的峰值个数

int peaks_num(mat sig, double minpeakdistance, double minpeakheight)
{
      /*
        只需要统计谱峰值即可
        minpeakeight-----最小峰值大小
        minpeakdistance---峰值之间的最小距离
      */
        int sig_length = sig.n_elem;
        vector<int>diff_v(sig_length-1);
        for (int i = 0; i < sig_length-1; i++)
        {
            if (sig(i+1)-sig(i)>0)
            {
                diff_v[i] = 1;
            }
            else if (sig(i+1)-sig(i)<0)
            {
                diff_v[i] = -1;
            }
            else
            {
                diff_v[i] = 0;
            }
        }
        //从尾部遍历向量，计算二阶差分运算
        for (int i = diff_v.size() - 1; i >= 0; i--)
        {
            if (diff_v[i] == 0 && i == diff_v.size() - 1)
            {
                diff_v[i] = 1;
            }
            else if (diff_v[i] == 0)
            {
                if (diff_v[i + 1] >= 0)
                    diff_v[i] = 1;
                else
                    diff_v[i] = -1;
            }
        }

        vector<int>peak_height;
        for (int i = 0; i < diff_v.size()-1; i++)
        {
            if ((diff_v[i+1]-diff_v[i]==-2)&&(sig(i+1)>=minpeakheight))
            {
                peak_height.push_back(i + 1);
            }
        }
        //对波峰之间的距离进行判断，从第二个到最后一个
        int peak_height_length = peak_height.size();
        int peak_num=0;  //计算波峰数量
        for (int i = 1; i < peak_height_length-1; i++)
        {
            if (peak_height[i]-minpeakdistance>=peak_height[i-1]&&peak_height[i]+minpeakdistance>=peak_height[i+1])
            {
                peak_num++;
            }
        }
        //判断第一个数据值
        if (peak_height[0]+minpeakdistance>=peak_height[1])
        {
            peak_num =peak_num+1;
        }
        //判断最后一个
        if (peak_height[peak_height_length-1]-minpeakdistance>=peak_height[peak_height_length-2])
        {
            peak_num = peak_num + 1;
        }
        return peak_num;
}

2.matlab的angle函数
需求是返回一维复数矩阵的位于区间【- π \pi π, π \pi π】相位角,实现原理非常简单。

mat angle(cx_mat sig)
{
    int sig_length = sig.n_elem;
    mat p = zeros(1, sig_length);
    for (size_t i = 0; i < sig_length; i++)
    {
        p(i) = atan2(imag(sig(i)),real(sig(i)));
    }
    return p;
}

3.matlab的fftshift函数
需求是需要进行频谱搬移 ，采用的还是for循环。

cx_mat fftshift(cx_mat fft_sig)
{
    /*
    分为奇数和偶数
    当矩阵长度为偶数时，前半部分为正频，后半部分为负频
    当矩阵长度为奇数时，前n/2+1为正频，后半部分为负频
    */
    mat fftshift_sig_real = zeros(1, fft_sig.n_elem);
    mat fftshift_sig_imag = zeros(1, fft_sig.n_elem);
    cx_mat fftshift_sig=cx_mat(fftshift_sig_real,fftshift_sig_imag);
    if (fft_sig.n_elem%2==0)
    {
        for (int i = 0; i < fft_sig.n_elem/2; i++)
        {
            fftshift_sig(i) = fft_sig(i + fft_sig.n_elem/2);
        }
        for (int i = fft_sig.n_elem/2; i < fft_sig.n_elem; i++)
        {
            fftshift_sig(i) = fft_sig(i - fft_sig.n_elem/2);
        }
    }
    else
    {
        for (int i = 0; i <(int) fft_sig.n_elem/2; i++)
        {
            fftshift_sig(i) = fft_sig(i + (int)fft_sig.n_elem/2+1);
        }
        for (int i = fft_sig.n_elem/2; i < fft_sig.n_elem; i++)
        {
            fftshift_sig(i) = fft_sig(i - (int)fft_sig.n_elem / 2);
        }
    }
    return fftshift_sig;
}

4.matlab的awgn函数
需求是对信号加噪，写了两个函数，分别是对实数信号和复数信号进行加噪

//实数信号的awgn函数
mat sig_add_noise(mat signal_type, Signal_parameter_T signal_parameter_T)
{
    double sigpower = 0;
    double noisepower = 0;
    mat nA = randn(1, signal_type.n_elem);  //考虑到都是一维数组

    for (int i = 0; i<signal_type.n_elem; i++)
    {
        sigpower += pow(signal_type(i), 2);
    }
    sigpower = sigpower / (signal_type.n_elem);    //得到信号强度
    noisepower = sigpower / (pow(10, (double)signal_parameter_T.sig_SNR/10));   //得到噪声强度
    mat noise = zeros(1, signal_type.n_elem);
    for (int i = 0; i<signal_type.n_elem; i++)
    {
        noise(i) = sqrt(noisepower)*nA(i);
    }
    return signal_type + noise;
}

//复数信号的awgn函数
cx_mat sig_add_noise2(cx_mat signal_type,int SNR)
{
    complex<double>J=(0,1);
    double sigpower = 0;
    double noisepower = 0;
    for(int j=0;j<signal_type.n_elem;j++)
    {
        sigpower+=abs(pow(signal_type(j),2));
    }
    sigpower = sigpower / (signal_type.n_elem);    //得到信号强度
    noisepower = sigpower / (pow(10, (double)SNR/10));   //得到噪声强度
    cx_mat noise=sqrt(noisepower/2)*(randn(1,signal_type.n_elem)+J*randn(1,signal_type.n_elem));
    return signal_type + noise;
}

5.matlab的希尔伯特变换

参考文章

cx_mat hilbert(mat sig) 
{
    int sig_length=sig.n_elem;
    cx_mat H(1,sig_length);
    cx_mat Y1;
    cx_mat Y2(1,sig_length);
    cx_mat Y2_temp;
    Y1=fft(sig,sig_length);
    H.col(0)=1;
    H.col(sig_length/2)=1;
    for(int i=1;i<sig_length/2;i++)
    {
        H.col(i)=2;
    }
    for(int j=sig_length/2+1;j<sig_length;j++)
    {
        H.col(j)=0;
    }
    Y2=Y1%H;  //矩阵点乘
    Y2_temp=ifft(Y2,sig_length);
    return Y2_temp;

}

6.matlab的firl1函数

参考文章

vector<double> sinc(vector<double> x)
{
    vector<double> y;
    for (int i = 0; i < x.size(); i++)
    {
        double temp = PI * x[i];
        if (fabs(temp) <= DBL_EPSILON) //if (temp == 0)
        {
            y.push_back(1.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;

}

当运算量比较大的时候，可以用动态数组对函数的运算速度进行改进

double sincEasy(double *x,int len,int index)
{
    double temp=PI*x[index];
    double y;
    if(fabs(temp) <= DBL_EPSILON)
    {
        y=1.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，但应注意传进一个数组中
                   
        参数输出介绍：
            fir1_improve: 动态数组，内存需要在在函数外进行销毁
    */
    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;
    }

    // 如果需要，现在可以处理缩放.
    bool scale=true;
    if (scale) {
        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;
}

7.matlab的filter函数

参考文章

mat filter(vector<double> b, int a, mat x)
{
    /*
    b:滤波器系数
    a:分母系数
    x:需滤波信号
    filter函数返回矩阵  长度跟x一致
    */

    vector<double> Y;
    for (int i = 0; i < b.size(); i++)
    {
        double real = 0.0;
        for (int j = 0; j <= i; j++)
        {
            real += b[j] * x(i - j);
        }

        Y.push_back(real);

    }

    for (int i = b.size(); i < x.n_elem; i++)
    {
        double real = 0.0;
        for (int j = 0; j < b.size(); j++)
        {
            real += (double)b[j] * x(i - j);
        }

        Y.push_back(real);
    }

    mat filter_result = zeros(1, x.n_elem);
    for (int i = 0; i < x.n_elem; i++)
    {
        filter_result(0, i) = Y[i];
    }

    return filter_result;
}

当运算两比较大时，可以通过动态数组进行改善

mat filter_darr(double lbf[],int lbflen,mat x)
{
    /*
       未写报错处理--当矩阵x的长度小于lbflen的时候，程序会崩
       参数输入介绍：
                    lbf: 滤波器系数
                    lbflen:滤波器系数的长度
                    x：输入矩阵
       参数输出介绍：
                   返回处理后的矩阵，长度与矩阵x一致
                   
    */
    double real=0.0;
    double *temp_iq=new double[x.n_elem];//存放结果的动态内存数组
    for(int i=0;i<lbflen;i++)
    {
        for(int j=0;i<i;j++)
        {
            real+=lbf[j]*x(i-j);
        }
        temp_iq[i]=real;
    }
    for(int i=lbflen;i<x.n_elem;i++)
    {
        double real=0.0;
        for(int j=0;j<lbflen;j++)
        {
            real+=(double)lbf[j]*x(i-j);
        }
        temp_iq[i]=real;
    }
    mat filter_result=zeros(1,x.n_elem);
    for(int i=0;i<x.n_elem;i++)
    {
        filter_result(0,i)=temp_iq[i];
    }
    delete[] temp_iq;
    return filter_result;
}
小河^^
关注
  • 10
    点赞
  • 踩
  • 30
    收藏
    觉得还不错? 一键收藏
  • 5
    评论
专栏目录
MATLAB转C/C++常见问题总结
ZYunfei的博客
05-05 5595
楔子 本指南试图帮助用户写出在转C过程中尽量少报错的Matlab m文件，其中更多是一些编程范式上的问题，在此记录下来，帮助后人写出无缝转C的代码。 Author: Y. F. Zhang Time: 2022.4.16 注意点 使用变量前声明变量 这一点在matlab中时可以容忍的，但对于一般编程语言是无法容忍的。使用一个变量前请预分配变量内存（注意预分配内存不是简单前面给个 =0就可以了，得是相应维度的预分配)： tmp = zeros(6, 6); % 预分配内存 for i = 1 : 6
将MATLAB程序代码转换成C++程序代码
09-06
将MATLAB的源程序转换成C++代码，不过转换之后，程序代码会增加很多，作辅助作用
如何将matlab代码生成为C/C++代码？
最新发布
weixin_49300040的博客
05-11 690
实际的matlab代码或者函数可能比较复杂，在根据coder命令提示进行问题检测中可能会遇到各种报错提示，此时无需惶恐，只有静下心去分析，大部分报错都是可以屏蔽或者修改的。采用Matlab自带的”coder”命令进行代码自动翻译，以达到省时省力的目的。
【MATLAB】MATLAB转C/C++语言并部署到VS2019
xiyang_DING
03-09 825
生成后会发现报错，这是因为不是所有语句都可以转换为C/C++，将不能转换的语句修改或者注释掉。Build type、Language以Toolchain根据自己需求选。More Setting中也根据自己需要选，注意最好选Debug。C/C++中需要给变量定义类型，选择指定类型。在APP中找到MATLAB Coder。打开VS2019，将生成代码文件导入。成功生成后C/C++的代码路径为。首先将你要转换的代码保存。
matlab代码转为c++流程
niuqq123的博客
10-04 1万+
这里写自定义目录标题欢迎使用Markdown编辑器新的改变功能快捷键合理的创建标题，有助于目录的生成如何改变文本的样式插入链接与图片如何插入一段漂亮的代码片生成一个适合你的列表创建一个表格设定内容居中、居左、居右SmartyPants创建一个自定义列表如何创建一个注脚注释也是必不可少的KaTeX数学公式新的甘特图功能，丰富你的文章UML 图表FLowchart流程图导出与导入导出导入 欢迎使用Markdown编辑器 你好！ 这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Mar
MATLAB转C语言的步骤与解决方案
热门推荐
weixin_44982036的博客
12-11 1万+
MATLAB转C语言的步骤 步骤 我们首先以求三阶矩阵的逆为例，首先在matlab中编写的代码如下： function[res] = inverse(src) res = inv(src) end 此为新建的一个函数，函数的保存名称为 inverse.m ,需要注意的是这个.m文件必须是函数文件，而不能是脚本文件，而且文件的名称和函数的名称必须一致。 在matlab中输入一个三阶矩阵与输出的结果如下所示： ` 在本文的末尾会验证在VS2017C语言代码的运行结果。 进入主题，在已经写好了函数的
Vc++下如何利用Matlab工具箱进行数字信号处理
qq_61141142的博客
09-25 157
编译器不支持图形函数，所以M文件中如有图形函数需注释掉。Matcom编译M文件，先将M文件按照与Matcom的cpp库的对应关系翻译为cpp源代码，然后用对应版本的c编译器将cpp文件编译成相应的exe或dll文件，所以第一次运行要指定c编译器的路径，否则无法编译，指定好的编译信息就写在Matcom\bin\matcom.ini文件中，不过这一步按装matcom时，它自动寻找编译器并将其写入matcom.ini文件中，matcom4.5版中使用TeeChart3.0 OCX控件，因而它支持图形操作。
如何将matlab代码转为C语言（2）--在C++中调用matlab的函数
weixin_44982036的博客
12-30 5651
如何将matlab代码转为C语言（2）–在C++中调用matlab的函数 在上一条博文中提供了一种直接在matlab操作中的方法，下面提供一种新的调用方法，即在C++中调用matlab中的dll文件。 首先就是环境的配置： 一、电脑的环境配置 1、右键我的电脑，点击属性，在控制面板主页中选择高级系统设置 2、点击环境变量->双击path->点击新建，添加路径为G:\MATLAB2019a\bin\win64（这是我电脑中matlab安装路径） 注意添加了路径之后一定要一路点击确定，这样才会配置成
3.6 使用matlab coder转换至c++代码的问题
South_Li的博客
03-06 311
出现(？：1)、数组溢出、维度不匹配等问题都可能是因为在标量扩展方面与MATLAB的不兼容性导致，可在官网查看原因，大部分问题在社区里还是有相似解决办法的
用matlab和visual C++开发信号分析处理的仿真与应用软件
07-12
例如，文件`digitalprocessusematlab.kdh`可能包含使用MATLAB编写的数字信号处理代码，如离散傅立叶变换（DFT）、快速傅立叶变换（FFT）以及滤波器设计等。 MATLAB的优势在于其交互性，但当需要高性能计算或者构建...
基于MATLAB的数字信号处理
08-01
"基于MATLAB的数字信号处理"这一主题涵盖了广泛的知识点，包括数字信号处理的基本概念、MATLAB编程基础以及如何在MATLAB环境中实现数字信号处理算法。以下是对这些知识点的详细阐述： 1. **数字信号处理基础**：...
matlab转C代码
05-04
matlab代码转为C语言在STM32上运行，该文档可以帮助了解过程及步骤
25种数字信号处理c++程序涉及到各种算法
05-15
25种数字信号处理c++程序涉及到各种算法如：离散傅里叶 实序数快速傅里叶 切比雪夫法 窗函数法等等
MATLAB数字信号处理仿真平台
10-05
MATLAB（Matrix Laboratory）是一款强大的数学计算和数据分析软件，尤其在数字信号处理（Digital Signal Processing, DSP）领域，它提供了全面的工具箱和强大的仿真环境。MATLAB数字信号处理仿真平台是工程师、科研...
使用MATLAB进行数字信号处理简介-第1部分-综合文档
05-20
在数字信号处理领域，MATLAB（Matrix Laboratory）是一款不可或缺的工具，它以其强大的计算能力、丰富的函数库和直观的编程环境深受广大工程师和科研人员喜爱。本教程将针对初学者，详细阐述如何利用MATLAB进行数字...
数字信号处理 王大伦.zip
07-10
MATLAB是一种交互式编程环境，特别适合于进行数值计算和数据可视化，因此在数字信号处理中被广泛应用。 MATLAB代码在数字信号处理中的作用主要体现在以下几个方面： 1. **信号生成**：MATLAB可以方便地生成各种...
matlab2c使用c++实现matlab函数系列教程-awgn函数
全栈工程师开发手册（原创）https://github.com/tencentmusic/cube-studio
09-06 6181
全栈工程师开发手册 （作者：栾鹏） matlab2c动态链接库下载 matlab库函数大全 matlab2c基础教程 matlab2c开发全解教程matlab2c调用方法：1、下载matlab2c动态链接库 2、将matlab2c.dll、matlab2c.lib和matlab2c.h放到项目头文件目录下 3、在cpp文件中引入下面的代码#include "Matlab2c.h" #pra
Matlab学习12——diff
i_can_do_nothing的博客
02-03 1225
Matlab学习12——diff
MATLAB实战：数字信号处理指南
本资源似乎是一本英文版的书籍，名为《Digital Signal Processing using MATLAB》，由 André Quinquis 编著，主要探讨了如何利用MATLAB进行数字信号处理，并提供了许多实际应用案例。 MATLAB（Matrix Laboratory）...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

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

余额充值