效果如下:
正所谓前任栽树后人乘凉 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
}