waveprocess.cpp

自制软件

// WaveProcess.cpp: implementation of the CWaveProcess class.
//   功能：完成信号处理(FFT)、产生数据文件
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Wave_Fft.h"
#include "WaveProcess.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

#define PI 3.1415
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
/*复数加法*/
COMPLEX Add(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re+c2.re;
	c.im=c1.im+c2.im;
	return c;
}

/*复数减法*/
COMPLEX Sub(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re-c2.re;
	c.im=c1.im-c2.im;
	return c;
}

/*复数乘法*/
COMPLEX Mul(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re*c2.re-c1.im*c2.im;
	c.im=c1.re*c2.im+c2.re*c1.im;
	return c;
}


CWaveProcess::CWaveProcess()
{

}

CWaveProcess::~CWaveProcess()
{

}
/************************************************************************/
/*               产生数据，默认为正弦                                   */
/************************************************************************/
void CWaveProcess::GetData(double *data,int Num,int type,int NoiseType,float Fre)
{
	int i;
	double length=0;
    double index;
	double tempFre=Fre;
	if(type==0)
	{	
		for (i=0;i<Num;i++) 
		{
			data[i]=180*sin(2*Fre*PI*i/Num);
		}
	}
	else
	{
		if (Fre!=0)
		{
			length = (double)(Num-1)/(double)(tempFre);
		}
	    else
		{
			length=2;
		}

		for(i=0;i<Num;i++) 
		{   
           	if(length>1)
			{	
				index = (double)(i)/length;
			}
			else
			{
				index = 0;
			}
	
			while(index>1)
				index -=1;
			if(index<=0.25)
				data[i] = index*800;
			else if(index<=0.75)
				data[i] = (0.5-index)*800;
			else
				data[i] = (index-1)*800;
		}
	}

	if(NoiseType==1)
	{	
		for (i=0;i<Num;i++) 
		{
			data[i]=data[i]+GetRand(-10,10);
		}
	}

}

/**************************************************************************
	功能：随机数
	参数：MIN为下界，MAX为上界
	返回：随机数
/**************************************************************************/

double  CWaveProcess::GetRand(int MIN, int MAX)
{
     int max;
     max=RAND_MAX;
     return (double)(rand()*(MAX-MIN)/max+MIN);
}

/***********************************************************************/
/* 	功能：实现快速傅立叶变换
/* 	参数：TD先存放时域值，后存放变换后的频率值
/* *********************************************************************/
void CWaveProcess::FFT(double Data[],double *FftData, int GapNum, int power)
{
	int i,j,k,bfsize,p;
	double angle;
	COMPLEX *W,*X1,*X2,*X;
    COMPLEX *TD;
    COMPLEX *FD;

	TD=new COMPLEX[GapNum];
	FD=new COMPLEX[GapNum];
	
	/*实数FFT变换*,将数组偶数和奇数分别赋予复数的实部和虚部*/
	for(i=0;i<GapNum/2;i++)
	{
		TD[i].re=Data[2*i];
		TD[i].im=Data[2*i+1];
	}

	GapNum=GapNum/2;
	/*分配运算所需存储器*/
	W=(COMPLEX *)malloc(sizeof(COMPLEX)*GapNum/2);
	X1=(COMPLEX *)malloc(sizeof(COMPLEX)*GapNum);
	X2=(COMPLEX *)malloc(sizeof(COMPLEX)*GapNum);
	
	/*计算加权系数*/
	for(i=0;i<GapNum/2;i++)
	{
		angle=-i*PI*2/GapNum;
		W[i].re=cos(angle);
		W[i].im=sin(angle);
	}
	
	/*将时域点写入存储器*/
	memcpy(X1,TD,sizeof(COMPLEX)*GapNum);
	
	/*蝶形运算*/
	for(k=0;k<power;k++)
	{
		for(j=0;j<1<<k;j++)
		{
			bfsize=1<<(power-k);
			for(i=0;i<bfsize/2;i++)
			{
				p=j*bfsize;
				X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
				X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
			}
		}
		X=X1;
		X1=X2;
		X2=X;
	}
	
	/*重新排序*/
	for(j=0;j<GapNum;j++)
	{
		p=0;
		for(i=0;i<power;i++)
		{
			if (j&(1<<i)) 
				p+=1<<(power-i-1);
		}
		FD[j]=X1[p];
	}
	GapNum=2*GapNum;

	/*实数FFT变换处理*/
    float *x0=new float[GapNum/2];
    float *x1=new float[GapNum/2];
 
    for(i=0;i<GapNum/2;i++)
    {
		x0[i]=FD[i].re;
		x1[i]=FD[i].im;
    }

    for(i=0;i<GapNum/2;i++)
    {
		x0[i]=FD[i].re;
		x1[i]=FD[i].im;
    }
	

    float wr=1.0;
    float wi=0;

    float theta=-8*atan(1.0)/GapNum;
    float ct=cos(theta);
    float st=sin(theta);
    	
    FD[GapNum/2].re=x0[0]-x1[0];
    FD[GapNum/2].im=0.0;
	
    FD[0].re=x0[0]+x1[0];
    FD[0].im=0.0;
	
    for(i=1;i<GapNum/4;i++)
    {
        float x0_1=x0[i];
        float x1_1=x1[i];
        float x0_2=x0[GapNum/2-i];
        float x1_2=x1[GapNum/2-i];
        float x0_1_=0.5*(x0_1+x0_2);
        float x0_2_=0.5*(x1_1-x1_2);
        float x1_1_=0.5*(x1_1+x1_2);
        float x1_2_=-0.5*(x0_1-x0_2);
        float  wr_t=wr;
		
        wr=wr_t*ct-wi*st;
        wi=wr_t*st+wi*ct;
        FD[i].re=x0_1_+wr*x1_1_-wi*x1_2_;
        FD[i].im=x0_2_+wr*x1_2_+wi*x1_1_;
        FD[GapNum/2-i].re=x0_1_-wr*x1_1_ +wi*x1_2_;
        FD[GapNum/2-i].im=-x0_2_+wr*x1_2_+wi*x1_1_;
    }
	
	FD[GapNum/4].re =  x0[GapNum/4];
	FD[GapNum/4].im = -x1[GapNum/4];
	
	for (i = 1; i < GapNum/2; i++)
	{
		FD[i + GapNum/2].re =  FD[GapNum/2 - i].re;
		FD[i + GapNum/2].im = -FD[GapNum/2 - i].im;
	}
	
	/*频谱分析*/
   for(i=0;i<GapNum;i++)
	{
		FftData[i]=pow(pow(FD[i].re,2)+pow(FD[i].im,2),0.5);
	}
	
	/*释放存储器*/
	delete []TD;
	delete []FD;
    delete[] x0;
    delete[] x1;
	free(W);
	free(X1);
	free(X2);

}