318 lines
9.1 KiB
C++
318 lines
9.1 KiB
C++
|
#include "Calculation.hpp"
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
Calculation::Calculation()
|
|||
|
|
{
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
Calculation::~Calculation()
|
|||
|
|
{
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
/************************************************************************/
|
|||
|
|
/* 一维数据的复数快速傅里叶变换 */
|
|||
|
|
/************************************************************************/
|
|||
|
|
void Calculation::FFT(int n, fftw_complex* in, fftw_complex* out)
|
|||
|
|
{
|
|||
|
|
if (in == NULL || out == NULL) return;
|
|||
|
|
fftw_plan p;
|
|||
|
|
p = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
|
|||
|
|
fftw_execute(p);
|
|||
|
|
fftw_destroy_plan(p);
|
|||
|
|
fftw_cleanup();
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/************************************************************************/
|
|||
|
|
/* 一维数据的实数快速傅里叶变换 */
|
|||
|
|
/************************************************************************/
|
|||
|
|
void Calculation::FFT_R(int n, std::vector<float> & vecData, fftw_complex* out)
|
|||
|
|
{
|
|||
|
|
double in[n];
|
|||
|
|
for (int i = 0; i < n; i++) {
|
|||
|
|
in[i] = vecData[i];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
for (int i = 0; i < n; i++)
|
|||
|
|
{
|
|||
|
|
out[i][0] = (double)vecData[i];
|
|||
|
|
out[i][1] = 0;
|
|||
|
|
}
|
|||
|
|
//create a DFT plan and execute it
|
|||
|
|
fftw_plan plan = fftw_plan_dft_r2c_1d(n, in, out, FFTW_ESTIMATE);
|
|||
|
|
fftw_execute(plan);
|
|||
|
|
//destroy the plan to prevent a memory leak
|
|||
|
|
fftw_destroy_plan(plan);
|
|||
|
|
fftw_cleanup();
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/************************************************************************/
|
|||
|
|
/* 一维数据的快速傅里叶逆变换 */
|
|||
|
|
/************************************************************************/
|
|||
|
|
void Calculation::iFFT(int n, fftw_complex* in, fftw_complex* out)
|
|||
|
|
{
|
|||
|
|
if (in == NULL || out == NULL) return;
|
|||
|
|
fftw_plan p;
|
|||
|
|
p = fftw_plan_dft_1d(n, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
|
|||
|
|
fftw_execute(p);
|
|||
|
|
fftw_destroy_plan(p);
|
|||
|
|
fftw_cleanup();
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
//************************************
|
|||
|
|
// Method: caculateAmp_Pha
|
|||
|
|
// FullName: Calculation::caculateAmp_Pha
|
|||
|
|
// Access: public static
|
|||
|
|
// Returns: void
|
|||
|
|
// Qualifier:
|
|||
|
|
// Parameter: int n
|
|||
|
|
// Parameter: fftw_complex * in
|
|||
|
|
// Parameter: int frequency
|
|||
|
|
// Parameter: double & amplitude
|
|||
|
|
// Parameter: double & phase
|
|||
|
|
// 函数功能是计算特定频率的幅值和相位，原来的讨论中是传入一个特定的频率，然后在给定的频率左右范围内找幅值和相位
|
|||
|
|
// 目前的函数实现是计算FFT变换后特定点的幅值和相位
|
|||
|
|
// 然后还有一个地方需要修改，即给定频率和FFT变换结果序列
|
|||
|
|
//************************************
|
|||
|
|
void Calculation::caculateAmp_Pha(int n, fftw_complex* in, int frequency, double &litude, double &phase)
|
|||
|
|
{
|
|||
|
|
int index = frequency;
|
|||
|
|
amplitude = 2 * sqrt((in[index][0] / n) * (in[index][0] / n) + (in[index][1] / n) * (in[index][1] / n));
|
|||
|
|
phase = 180 * atan(in[index][1] / in[index][0]) / M_PI;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
float Calculation::max(std::vector<float> & vecData)
|
|||
|
|
{
|
|||
|
|
std::vector<float>::iterator it;
|
|||
|
|
it = std::max_element(vecData.begin(), vecData.end());
|
|||
|
|
if (it != vecData.end()) {
|
|||
|
|
return *it;
|
|||
|
|
}
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
float Calculation::min(std::vector<float> & vecData)
|
|||
|
|
{
|
|||
|
|
std::vector<float>::iterator it;
|
|||
|
|
it = std::min_element(vecData.begin(), vecData.end());
|
|||
|
|
if (it != vecData.end()) {
|
|||
|
|
return *it;
|
|||
|
|
}
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
void Calculation::absVec(std::vector<float> & vecAbsData, std::vector<float> & vecData)
|
|||
|
|
{
|
|||
|
|
for (int i = 0; i < vecData.size(); i++) {
|
|||
|
|
vecAbsData.push_back(fabs(vecData[i]));
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
float Calculation::mean(std::vector<float> & vecData)
|
|||
|
|
{
|
|||
|
|
double meanTemp = 0;
|
|||
|
|
for (int i = 0; i < vecData.size(); i++) {
|
|||
|
|
meanTemp = meanTemp += vecData[i];
|
|||
|
|
}
|
|||
|
|
return meanTemp / vecData.size();
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
void Calculation::drop_mean(std::vector<float> & vecDropMeanData, std::vector<float> & vecData)
|
|||
|
|
{
|
|||
|
|
float fMean = mean(vecData);
|
|||
|
|
for (int i = 0; i < vecData.size(); i++) {
|
|||
|
|
vecDropMeanData.push_back(vecData[i] - fMean);
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
float Calculation::srm(std::vector<float> & vecData)
|
|||
|
|
{
|
|||
|
|
double dSrmTemp = 0;
|
|||
|
|
for (int i = 0; i < vecData.size(); i++){
|
|||
|
|
dSrmTemp = dSrmTemp + sqrt(vecData[i]);
|
|||
|
|
}
|
|||
|
|
dSrmTemp = dSrmTemp / vecData.size();
|
|||
|
|
return dSrmTemp * dSrmTemp;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
float Calculation::rms(std::vector<float> & vecData)
|
|||
|
|
{
|
|||
|
|
double rmsTemp = 0;
|
|||
|
|
for (int i = 0; i < vecData.size(); i++) {
|
|||
|
|
rmsTemp = rmsTemp += (vecData[i] * vecData[i]);
|
|||
|
|
}
|
|||
|
|
rmsTemp = rmsTemp / vecData.size();
|
|||
|
|
return sqrt(rmsTemp);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
float Calculation::variance(std::vector<float> & vecDropMeanData)
|
|||
|
|
{
|
|||
|
|
double varianceTemp = 0;
|
|||
|
|
for (int i = 0; i < vecDropMeanData.size(); i++) {
|
|||
|
|
varianceTemp = varianceTemp += (vecDropMeanData[i] * vecDropMeanData[i]);
|
|||
|
|
}
|
|||
|
|
return varianceTemp/vecDropMeanData.size();
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
float Calculation::skew_state(std::vector<float> & vecDropMeanData, float fVariance)
|
|||
|
|
{
|
|||
|
|
double tempSkew = 0;
|
|||
|
|
for (int i = 0; i < vecDropMeanData.size(); i++) {
|
|||
|
|
tempSkew = tempSkew + pow(vecDropMeanData[i], 3);
|
|||
|
|
}
|
|||
|
|
tempSkew = tempSkew / vecDropMeanData.size();
|
|||
|
|
tempSkew = tempSkew / pow(fVariance, 1.5);
|
|||
|
|
return tempSkew;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
float Calculation::kurtosis(std::vector<float> & vecDropMeanData, float fVariance)
|
|||
|
|
{
|
|||
|
|
double tempkurtosis = 0;
|
|||
|
|
for (int i = 0; i < vecDropMeanData.size(); i++) {
|
|||
|
|
tempkurtosis = tempkurtosis + pow(vecDropMeanData[i], 4);
|
|||
|
|
}
|
|||
|
|
tempkurtosis = tempkurtosis / vecDropMeanData.size();
|
|||
|
|
tempkurtosis = tempkurtosis / pow(fVariance, 2);
|
|||
|
|
return tempkurtosis;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void Calculation::hilbert(std::vector<float> & vecData, std::vector<float> & vecHilbertData, int N)
|
|||
|
|
{
|
|||
|
|
|
|||
|
|
double in[N];
|
|||
|
|
for (int i = 0; i < N; i++) {
|
|||
|
|
in[i] = vecData[i];
|
|||
|
|
}
|
|||
|
|
fftw_complex *out;
|
|||
|
|
out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * N);
|
|||
|
|
|
|||
|
|
for (int i = 0; i < N; ++i)
|
|||
|
|
{
|
|||
|
|
out[i][0] = (double)vecData[i];
|
|||
|
|
out[i][1] = 0;
|
|||
|
|
}
|
|||
|
|
//create a DFT plan and execute it
|
|||
|
|
fftw_plan plan = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
|
|||
|
|
fftw_execute(plan);
|
|||
|
|
//destroy the plan to prevent a memory leak
|
|||
|
|
fftw_destroy_plan(plan);
|
|||
|
|
int hN = N >> 1; // half of the length (N /2)
|
|||
|
|
int numRem = hN; // the number of remaining elements
|
|||
|
|
// multiply the appropriate values by 2
|
|||
|
|
// (those that should be multiplied by 1 are left intact because they wouldn't change)
|
|||
|
|
for (int i = 1; i < hN; ++i) // 1,2,...,N/2 - 1 的项乘以2
|
|||
|
|
{
|
|||
|
|
out[i][0] *= 2;
|
|||
|
|
out[i][1] *= 2;
|
|||
|
|
}
|
|||
|
|
// if the length is even, the number of remaining elements decreases by 1
|
|||
|
|
if (N % 2 == 0)
|
|||
|
|
numRem--;
|
|||
|
|
// if it's odd and greater than 1, the middle value must be multiplied by 2
|
|||
|
|
else if (N > 1) // 奇数非空
|
|||
|
|
{
|
|||
|
|
out[hN][0] *= 2;
|
|||
|
|
out[hN][1] *= 2;
|
|||
|
|
}
|
|||
|
|
// set the remaining values to 0
|
|||
|
|
// (multiplying by 0 gives 0, so we don't care about the multiplicands)
|
|||
|
|
memset(&out[hN + 1][0], 0, numRem * sizeof(fftw_complex));
|
|||
|
|
// create an IDFT plan and execute it
|
|||
|
|
plan = fftw_plan_dft_1d(N, out, out, FFTW_BACKWARD, FFTW_ESTIMATE);
|
|||
|
|
fftw_execute(plan);
|
|||
|
|
// do some cleaning
|
|||
|
|
fftw_destroy_plan(plan);
|
|||
|
|
fftw_cleanup();
|
|||
|
|
// scale the IDFT output
|
|||
|
|
for (int i = 0; i < N; ++i)
|
|||
|
|
{
|
|||
|
|
out[i][0] /= N;
|
|||
|
|
out[i][1] /= N;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
for( int n=0; n<N; n++ )//输出
|
|||
|
|
{
|
|||
|
|
// xr[n]=cos(n*pi/6);//原始信号
|
|||
|
|
// y_r[n] = s_i[n];
|
|||
|
|
complex complex_after;
|
|||
|
|
complex_after.real = out[n][1];
|
|||
|
|
complex_after.imag = out[n][0];
|
|||
|
|
float amp = sqrt(complex_after.real * complex_after.real +complex_after.imag * complex_after.imag);
|
|||
|
|
vecHilbertData.push_back(amp);
|
|||
|
|
// printf("%d %f\n",n,vecHilbertData[n]);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
void Calculation::fftShift(fftw_complex* in, int l)
|
|||
|
|
{
|
|||
|
|
double temp;
|
|||
|
|
double temp2;
|
|||
|
|
|
|||
|
|
for (int j = 0;j<l/2;j++) {
|
|||
|
|
temp = in[j+l/2][0];
|
|||
|
|
temp2 = in[j+l/2][1];
|
|||
|
|
in[j+l/2][0] = in[j][0];
|
|||
|
|
in[j+l/2][1] = in[j][1];
|
|||
|
|
in[j][0] = temp;
|
|||
|
|
in[j][1] = temp2;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void Calculation::FFTSpec(std::vector<float> & vecData, std::vector<float> & vecFFTSpecData)
|
|||
|
|
{
|
|||
|
|
fftw_complex *inFFt, *outFFt;
|
|||
|
|
inFFt = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * vecData.size());
|
|||
|
|
outFFt = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * vecData.size());
|
|||
|
|
printf("1---------------------------------------------->%d\n",vecData.size());
|
|||
|
|
|
|||
|
|
for (int j = 0; j < vecData.size(); j++) {
|
|||
|
|
inFFt[j][0] = (double)vecData[j];
|
|||
|
|
inFFt[j][1] = 0;
|
|||
|
|
}
|
|||
|
|
FFT(vecData.size(),inFFt, outFFt);
|
|||
|
|
for(int j = 0; j < vecData.size()/2; j++) {
|
|||
|
|
vecFFTSpecData.push_back(sqrt(outFFt[j][0]*outFFt[j][0] + outFFt[j][1]*outFFt[j][1])*2/vecData.size());
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void Calculation::envSpec(std::vector<float> & vecData, std::vector<float> & vecEnvSpecData)
|
|||
|
|
{
|
|||
|
|
std::vector<float> vecDropMeanData;
|
|||
|
|
drop_mean(vecDropMeanData, vecData);
|
|||
|
|
|
|||
|
|
std::vector<float> vecHilbertData;
|
|||
|
|
hilbert(vecDropMeanData, vecHilbertData, vecDropMeanData.size());
|
|||
|
|
|
|||
|
|
fftw_complex *inHilFFt, *outHilFFt;
|
|||
|
|
inHilFFt = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * vecHilbertData.size());
|
|||
|
|
outHilFFt = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * vecHilbertData.size());
|
|||
|
|
|
|||
|
|
for (int j = 0; j < vecHilbertData.size(); j++) {
|
|||
|
|
inHilFFt[j][0] = (double)vecHilbertData[j];
|
|||
|
|
inHilFFt[j][1] = 0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
FFT(vecHilbertData.size(), inHilFFt, outHilFFt);
|
|||
|
|
fftShift(outHilFFt, vecHilbertData.size());
|
|||
|
|
|
|||
|
|
vecEnvSpecData.push_back(0);
|
|||
|
|
for (int i = vecHilbertData.size() / 2 + 1; i < vecHilbertData.size(); i++) {
|
|||
|
|
float ftemp = 2 * sqrt(outHilFFt[i][0] * outHilFFt[i][0] + outHilFFt[i][1] * outHilFFt[i][1]) / vecHilbertData.size();
|
|||
|
|
vecEnvSpecData.push_back(ftemp);
|
|||
|
|
}
|
|||
|
|
}
|