效果如下:
在这里插入图片描述在这里插入图片描述
正所谓前任栽树后人乘凉 sp++如何引用不再赘述,由于首次接触信号处理方面的东西(知其然不知其所以然)结果也是试出来的。如有错误感谢指正!!!
频谱图
头文件

#ifndef signal_processing_H
#define signal_processing_H
#include "matrix.h"
#include <QMutex>
#include<qsemaphore.h>
#include <QObject>
#include <vectormath.h>

class signal_processing
{
public:
    /*!
        @Function    : 窗函数结构体
        @Parameter   : name 窗类型 Rectangle Bartlett Hanning Hamming Blackman
        @Parameter   : length 窗长 常用取值128 256 512 1024 2048 4096 不宜超过16384
        @Parameter   : amp
        @Parameter   : alpha
    */
    typedef struct windows_type_t
    {
        QString name="Hamming";
        int length=1024; //数据量太小时窗长不宜过长 窗长*100<数据量最合适
        double amp=1.0;//幅度
        double alpha=1.0;
    }windows_type;

    signal_processing();
    /*!
    @Function    : 获取采样率
    @Description :
    @Parameter   :
    @Output      :
    @Date        : 2023-12-12
*/
    inline double get_sampling_rate()const{return sampling_rate;}
    /*!
        @Function    : 设置采样率
        @Description :
        @Parameter   : value 设置采样率
        @Date        : 2023-12-12
    */
    void set_sampling_rate(const double value);
    /*!
        @Function    : 计算功率谱
        @Description : 开窗类型和开窗长度都会对结果产生影响,窗长决定 fft_size的计算结果
        @Parameter   : splab::Vector<double> &d 容器中存放需要计算的信号
        @Parameter   : QVector<double> &key 频率(MHz) -x轴
        @Parameter   : QVector<double> &value 功率值(dB) -y轴
        @Date        : 2023-12-12
    */
    bool calc_power_spectrum(splab::Vector<double> &d,QVector<double> &key,QVector<double> &value,const windows_type m_window);
    /*!
        @Function    : 计算时频
        @Description : 开窗类型和开窗长度都会对结果产生很大影响,计算结果是一个二位数组(行代表时间 列代表频率)
        @Parameter   : splab::Vector<double> &d 容器中存放需要计算的信号
        @Parameter   : splab::Matrix<double> &value 计算结果 实部real()有效
        @Parameter   : windows_type m_window 开窗信息
        @Date        : 2023-12-12
    */
    bool calc_speech_frequency_addWindowFourier(splab::Vector<double> &d,splab::Matrix< splab::complex<double>> &value,const windows_type m_window);

    ~signal_processing();
//private:
    int cal_numberOfSamples(const int wn_size);
    double cal_sinr(splab::Vector<double> &data); //计算SINR
    void switchToDb(splab::Matrix< splab::complex<double>> &value);
    double sampling_rate=0;//采样率(MHZ)
    double calc_window_enbw(splab::Vector<double> wn);//计算窗口等效噪声带宽
};

#endif // signal_processing_H

实现文件 (.cpp)

#include "signal_processing.h"
#include "wft.h"
#include <QDebug>
#include <iostream>
#include <cstring>
#include <parametricpse.h>
#include <classicalpse.h>
#include <eigenanalysispse.h>
#include <vectormath.h>
#include <dwt.h>
# include <wiener.h>
using namespace std;
using namespace splab;
#define MIN_PONIT_SIZE 128



typedef double  Type;
const   int     Fs = 1000;
const   int     Ls = 10000;
const   int     Lg = 80;
const   int     N  = 40;
const   int     dM = 10;      // over sampling ratio is N/dM
#include <QFile>
Vector<Type> st;

signal_processing::signal_processing()
{
    sampling_rate=1;//单位(MHz)
    Type a = 0;
    Type b = Ls-1;
    Vector<Type> t = linspace( a, b, Ls ) / Type(Fs);
    st = cos( Type(400*PI) * pow(t,Type(2.0)) );

    Type f1 = Type(0.3), f2 = Type(0.4);
    Type amp1 = Type(1.0),amp2 = Type(0.5);
    Vector<Type> tn = linspace(Type(0), Type(100-1),100 );
    Vector<Type> sn = amp1*sin(TWOPI*f1*tn) + amp2*sin(TWOPI*f2*tn);


}

signal_processing::~signal_processing()
{

}

int signal_processing::cal_numberOfSamples(const int wn_size)
{
    int L=128;
    if(L>wn_size)
    {
        while (L<wn_size)
        {
            L=L<<1;
            if(L>wn_size)
            {
                L=L>>1;
                break;
            }
        }
    }

    return L+1;
}

double signal_processing::cal_sinr(splab::Vector<double> &data)
{
    const int N =data.size();
    double Pn = 0.01; // 假设噪声功率
    double Pi = 0.1; // 假设干扰功率
    double Ps = 0; // 信号功率
    double sinr; // SINR
    for (int i = 0; i < N; i++)
    {
        Ps+=pow(data[i],2);
    }
    Ps /= N;
    sinr = 10 * log10(Ps / (Pn + Pi));
    return sinr;

}

void signal_processing::switchToDb(splab::Matrix<splab::complex<double> > &value)
{
    double max_v=-150;
    for(int row=0;row<value.rows();row++)
    {
        for(int column=0;column<value.cols();column++)
        {
            double mag = std::abs(value[row][column]);
            value[row][column] = mag;
            if (mag > max_v)
            {
                max_v = mag;
            }

        }
    }

    for(int row=0;row<value.rows();row++)
    {
        for(int column=0;column<value.cols();column++)
        {
            double val = std::abs(value[row][column]) / max_v;
            value[row][column] = 20 * std::log10(val);
        }
    }
}

double signal_processing::calc_window_enbw(splab::Vector<double> wn)
{
    double sum_1 = sum(wn);
    wn*=wn;
    double sum_2 = sum(wn);
    return  sum_2/(sum_1*sum_1*wn.size());
}


void signal_processing::set_sampling_rate(const double value)
{
    sampling_rate=value;
}


bool signal_processing::calc_power_spectrum(splab::Vector<double> &d, QVector<double> &key, QVector<double> &value, const windows_type m_window)
{
    int d_size=d.size();
    if(d_size<MIN_PONIT_SIZE)return false;

    int wn_length=m_window.length<d_size?m_window.length:(d_size-1);//计算窗长
    int fft_size = wn_length + 1;


    Vector<double> wn =window(m_window.name.toStdString(),wn_length, double(1.0));

    Vector<double> Ps = welchPSE(d,wn,wn_length/2,fft_size); //计算结果为功率密度 (如果要为dB需要转换)

    int result_size=Ps.size()/2;
    value.resize(result_size);
    key.resize(result_size);

    //转换为dB

    double enbw=calc_window_enbw(wn);

    for(int i=0;i<result_size;++i)
    {
        value[i]=10*log10(Ps[i]*enbw);
        key[i]=i/(double)result_size*(sampling_rate/2);
    }

    Ps.destroy();
    return true;
}

#include <dgt.h>
#include <wvd.h>

#include <wiener.h>
//语频图如何计算
bool signal_processing::calc_speech_frequency_addWindowFourier(splab::Vector<double> &d, splab::Matrix<complex<double>> &value,const windows_type m_window)
{
    int d_size=d.size();
    if(d_size<MIN_PONIT_SIZE)return false;

    Vector<double> wn =window(m_window.name.toStdString(),m_window.length, double(1.0));

    value = wft(d,wn,"ppd");

    switchToDb(value);

    return true;
}

原来的wft函数计算量太大,对wft的源码进行修改,添加滑动步长

template <typename Type>
Matrix< complex<Type> > wft( const Vector<Type> &xn, const Vector<Type> &wn,
                             const string &mode )
{
#if 0
    int Lx = xn.size(),
        Lw = wn.size();

    // extends the input signal
    Vector<Type> tmp = wextend( xn, Lw/2, "both", mode );

    Matrix< complex<Type> > coefs( Lw, Lx );
    Vector<Type> sn( Lw );
    Vector< complex<Type> > Sk( Lw );

    for( int i=0; i<Lx; ++i )
    {
        // intercept xn by wn function
        for( int j=0; j<Lw; ++j )
            sn[j] = tmp[i+j] * wn[j];
        Sk = fft(sn);

        // compute the Foureier transform
        coefs.setColumn( Sk, i );
    }

    return coefs;
#else
    int Lx = xn.size(),
        Lw = wn.size();

    // extends the input signal
    Vector<Type> tmp = wextend( xn, Lw/2, "both", mode );

    int step = Lw/2; // 窗口滑动的步长
    int Lcoefs = (Lx - Lw) / step + 1; // 计算矩阵的列数
    Matrix< complex<Type> > coefs( Lw, Lcoefs ); // 矩阵的列数为Lcoefs
    //Matrix< complex<Type> > coefs( Lw, Lx );
    Vector<Type> sn( Lw );
    Vector< complex<Type> > Sk( Lw );

    for( int i=0; i<Lcoefs; ++i )
    {
         int start = i * step; // 计算起始位置
        // intercept xn by wn function
        for( int j=0; j<Lw; ++j )
        {
            sn[j] = tmp[start+j] * wn[j];
        }
        Sk = fft(sn);
        coefs.setColumn( Sk, i );
    }

    return coefs;
#endif
}
Logo

瓜分20万奖金 获得内推名额 丰厚实物奖励 易参与易上手

更多推荐