matlab函数转C++(数字信号处理)-Toy模板网

这篇具有很好参考价值的文章主要介绍了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;
}