1. Toy模板网首页
  2. > Toy博客

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

9月前 作者：小河^^ 分类：Toy博客 阅读(40) 违法举报

这篇具有很好参考价值的文章主要介绍了matlab函数转C++(数字信号处理)。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

前言

近期主要利用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;
}

当运算两比较大时，可以通过动态数组进行改善文章来源地址https://www.toymoban.com/news/detail-740202.html

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;
}

到了这里，关于matlab函数转C++(数字信号处理)的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处： 如若内容造成侵权/违法违规/事实不符，请点击违法举报进行投诉反馈，一经查实，立即删除！

分享到：
领支付宝红包 赞助服务器费用

相关文章

  • 【数字信号处理】带通采样定理及其MATLAB仿真

    【数字信号处理】带通采样定理及其MATLAB仿真

    按照奈奎斯特采样定理(低通采样)，采样频率 f s f_{s} f s ​ 要大于等于信号中最高频率 f m a x f_{max} f ma x ​ 的2倍，才可以保证采样后的数字信号通过DAC转换后，可以无失真的恢复为原信号。然而，如果信号的频率分布在某一有限频带上，并且信号的最高频率 f m a x f_{max} f

    2024年02月16日
    浏览(48)
  • 信号处理之FIR数字滤波器（Matlab仿真）

    信号处理之FIR数字滤波器（Matlab仿真）

            数字滤波器的作用是滤除不感兴趣的信号，留下想要的信号。数字滤波器可分为无限脉冲响应(IIR)数字滤波器、有限脉冲响应(FIR)数字滤波器两种，两者各有优缺点，其中FIR数字滤波器因其具有良好的线性相位特性受到广泛应用，线性相位是指信号中各频率成分的相对

    2024年02月03日
    浏览(50)
  • 数字信号处理实验---LSI系统的分析 Matlab代码

    数字信号处理实验---LSI系统的分析 Matlab代码

    1.试用Matlab计算其幅频特性和相频特性，并绘图。 代码： n = 0:10; %定义采样点n w = [0:1:500]*2*pi/500; % [0,pi]轴被分成1002个点 x1 = power(0.9*exp(1i*pi/3),n); %定义输入序列 x2 = exp(-1i*n); %定义一个系统的冲激响应 x = zeros(1,length(w)); %定义空数组存储系统的频域响应 for i=1:length(x1)     x=x

    2024年01月15日
    浏览(44)

  • matlab数字信号处理实验（5）时域采样与频域采样

    一、实验目的 1、理解时域采样理论与频域采样理论； 2、掌握模拟信号采样前后频谱的变化，以及如何选择采样频率才能使采样后的信号 不丢失信息； 3、掌握频率域采样会引起时域周期化的原因，频率域采样定理及其对频域采样点数 选择的指导作用； 4、对信号在某个表示

    2024年02月07日
    浏览(43)
  • 【数字化处理】仿生假体控制中肌电信号的数字化处理研究（Matlab代码实现）

    【数字化处理】仿生假体控制中肌电信号的数字化处理研究（Matlab代码实现）

     💥💥💞💞 欢迎来到本博客 ❤️❤️💥💥 🏆博主优势： 🌞🌞🌞 博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。 ⛳️ 座右铭： 行百里者，半于九十。 📋📋📋 本文目录如下： 🎁🎁🎁 目录 💥1 概述 📚2 运行结果 🎉3 参考文献 🌈4 Matlab代码及文献 文献

    2024年02月12日
    浏览(42)
  • 一、信号处理 ——impseq函数与stepseq函数(Matlab实现)

    一、信号处理 ——impseq函数与stepseq函数(Matlab实现)

    在脚本中直接运行一次即可，在matlab左侧生成impseq.m文件与stepseq.m文件 1. 单位脉冲函数impseq. 2. 单位阶跃函数stepseq. 仅用于学习记录~

    2024年02月06日
    浏览(38)
  • 数字信号处理-10-并行FIR滤波器MATLAB与FPGA实现

    数字信号处理-10-并行FIR滤波器MATLAB与FPGA实现

    本文介绍了设计滤波器的FPGA实现步骤，并结合杜勇老师的书籍中的并行FIR滤波器部分进行一步步实现硬件设计，对书中的架构做了复现以及解读，并进行了仿真验证。 FIR滤波器的结构形式时，介绍了直接型、级联型、频率取样型和快速卷积型4种。在FPGA实现时，最常用的是最

    2023年04月09日
    浏览(48)
  • 数字信号处理篇之浮点数与定点数的转换（MATLAB）

    数字信号处理篇之浮点数与定点数的转换（MATLAB）

      对于计算机等数字信号处理器件，数字和信号变量都是用二进制进行表示的。在本文中，我们学习了定点数的概念、浮点数与定点数的转换以及在MATLABZ中实现浮点数与定点数的转换。   对于二进制数，大家应该都很熟悉，在学习数电的过程中，我们知道，十进制转二

    2024年02月11日
    浏览(58)
  • 数字信号处理实验---Z变换及系统的零极点分析 Matlab代码

    数字信号处理实验---Z变换及系统的零极点分析 Matlab代码

    一．各种函数的用法 1.tf2zp函数：通常用于将传递函数（Transfer Function）转换为零极增益形式（ZPK form），转换前G(s) = num(s) / den(s)，转换后G(s) = K * (s - z1) * (s - z2) * ... * (s - zn) / (s - p1) * (s - p2) * ... * (s - pn) 2.zp2tf函数：用于将零极增益形式（ZPK form）转换为传递函数（Transfer Fu

    2024年01月23日
    浏览(47)
  • C++ 信号处理

    C++ 信号处理

    信号是由操作系统传给进程的中断，会提早终止一个程序。在 UNIX、LINUX、Mac OS X 或 Windows 系统上，可以通过按 Ctrl+C 产生中断。 有些信号不能被程序捕获，但是下表所列信号可以在程序中捕获，并可以基于信号采取适当的动作。这些信号是定义在 C++ 头文件 csignal 中。 信号

    2024年02月12日
    浏览(44)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包，优惠每天领

二维码1

领取红包

二维码2

领红包