mdifftdoc.cpp

对简单的一维信号进行FFT核IFFT处理

// mdiFFTDoc.cpp : implementation of the CMdiFFTDoc class
//

#include "stdafx.h"
#include "mdiFFT.h"
#include <math.h>
//#include <complex>

#include "mdiFFTDoc.h"
#include "ChildFrm.h"

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



/***************************************************/
//以下为fft与ifft函数。FFT的返回值为频谱的最大模值，便于显示时调整大小。

long double FFT(Complex *TD, Complex *FD, int r)
{
    LONG  count;

    int  i,j,k;

    int  bfsize,p;

    double  angle;

    Complex*W,*X1,*X2,*X;

    count=1<<r;

    W=new Complex[count/2];
    X1=new Complex[count];
    X2=new Complex[count];

    for(i=0;i<count/2;i++)
    {
      angle=-i*PI*2/count;
      W[i]=Complex(cos(angle),sin(angle));
    }

    memcpy(X1,TD,sizeof(Complex)*count);

    for(k=0;k<r;k++)
    {
        for(j=0;j<1<<k;j++)
        {
            bfsize=1<<(r-k);
            for(i=0;i<bfsize/2;i++)
            {
               p=j*bfsize;
               X2[i+p]=X1[i+p]+X1[i+p+bfsize/2];
               X2[i+p+bfsize/2]=(X1[i+p]-X1[i+p+bfsize/2])*W[i*(1<<k)];
            }
         }
         X=X1;
         X1=X2;
         X2=X;
    }

	double _max=0;
    for(j=0;j<count;j++)
    {
       p=0;
       for(i=0;i<r;i++)
       {
          if(j&(1<<i))
          { 
             p+=1<<(r-i-1);
          }
       }
       FD[j]=X1[p];
	   if(FD[j].abs()>_max)_max=FD[j].abs();
    }

    delete  W;
    delete  X1;
    delete  X2;
	return _max;
}

 

VOID IFFT(Complex FD[], Complex TD[], int r)
{
    LONG  count;

    int  i;

    Complex*X;

    count=1<<r;

    X=new Complex[count];

    memcpy(X,FD,sizeof(Complex)*count);

    for(i=0;i<count;i++)
    {
      X[i]=Complex(X[i].Re,-X[i].Im);
    }

    FFT(X,TD,r);

    for(i=0;i<count;i++)
    {
      TD[i]=Complex(TD[i].Re/count,-TD[i].Im/count);
    }

    delete X;
}

/*****************************************************************/

/////////////////////////////////////////////////////////////////////////////
// CMdiFFTDoc

IMPLEMENT_DYNCREATE(CMdiFFTDoc, CDocument)

BEGIN_MESSAGE_MAP(CMdiFFTDoc, CDocument)
	//{{AFX_MSG_MAP(CMdiFFTDoc)
	ON_BN_CLICKED(IDC_RADIO5, OnRadio5)
	ON_BN_CLICKED(IDC_RADIO6, OnRadio6)
	ON_BN_CLICKED(IDC_RADIO7, OnRadio7)
	ON_BN_CLICKED(IDC_RADIO8, OnRadio8)
	ON_NOTIFY(UDN_DELTAPOS, IDC_SPIN2, OnDeltaposSpin1)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CMdiFFTDoc construction/destruction

CMdiFFTDoc::CMdiFFTDoc()
{
	// TODO: add one-time construction code here
	m_power=8;
	srand((unsigned)GetTickCount%5000);
	Setsource(Sin);
	CalcValues();
}

CMdiFFTDoc::~CMdiFFTDoc()
{
}

BOOL CMdiFFTDoc::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;

	// TODO: add reinitialization code here
	// (SDI documents will reuse this document)

	return TRUE;
}
//*****************************************************************************//
//此成员函数设置原序列的值，其中Random为随机值
void CMdiFFTDoc::Setsource(MySourceType type)
{
	int i;
	switch(type){
	case Sin:
		{
			//const double PI=acos(-1);
			for(i=0; i<MaxData; i++)
			{
				source[i].Re=/*(DataType)(*/ 100*sin(i*4*PI/MaxData) ;
				source[i].Im=0;
			}
		}break;
	case Traingle:
		{
			double j;
			for(i=0; i<MaxData; i++)
			{
				j=4*i%MaxData*200/MaxData;
				if(j>100)j=200-j;
				if(j<-100)j=-j-200;
				source[i].Re= j ;
				source[i].Im=0;
			}break;
		}
	case Squar:
		for(i=0; i<MaxData; i++)
		{
			source[i].Re= -200*((int)(i*4/MaxData)%2)+100 ;
			source[i].Im=0;
		}break;
	case Random:	///原序列为随机数量（3～10）的随机频率的正弦波随机比例叠加
		{
			int count=rand()*7/RAND_MAX+3;		///正弦波数量
			int j,t;
			double period;
			for(j=0;j<MaxData;j++)
			{
				source[j].Im=0;
				source[j].Re=0;
			}
			for(i=0;i<count;i++)
			{
				period=(rand()*50/RAND_MAX);	///每个分量的周期
				for(j=0;j<MaxData;j++)
				{
					t=(1.0*rand()/RAND_MAX)*(200.0/count)*sin(j*4*PI/period);		///随机的权重
					if(t>100){
						t=rand()*50/RAND_MAX;
					}
					source[j].Re+=t;		///叠加

				}
			}
		}
	}
	CalcValues();
}
void CMdiFFTDoc::CalcValues()			//计算fft及ifft后的值
{
	m_amax=FFT(source,afterfft,m_power);
	IFFT(afterfft,afterifft,m_power);

}
//*******************************************************************************//

/////////////////////////////////////////////////////////////////////////////
// CMdiFFTDoc serialization

void CMdiFFTDoc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO: add storing code here
	}
	else
	{
		// TODO: add loading code here
	}
}

/////////////////////////////////////////////////////////////////////////////
// CMdiFFTDoc diagnostics

#ifdef _DEBUG
void CMdiFFTDoc::AssertValid() const
{
	CDocument::AssertValid();
}

void CMdiFFTDoc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CMdiFFTDoc commands


void CMdiFFTDoc::OnRadio5() 
{
	// TODO: Add your control notification handler code here
	Setsource(Sin);
	UpdateAllViews(NULL);
	
}

void CMdiFFTDoc::OnRadio6() 
{
	// TODO: Add your control notification handler code here
	Setsource(Traingle);
	UpdateAllViews(NULL);
}

void CMdiFFTDoc::OnRadio7() 
{
	// TODO: Add your control notification handler code here
	Setsource(Squar);
	UpdateAllViews(NULL);
	
}

void CMdiFFTDoc::OnRadio8() 
{
	// TODO: Add your control notification handler code here
	Setsource(Random);
	UpdateAllViews(NULL);
	
}

/*void CMdiFFTDoc::OnChangeEdit1() 
{
	// TODO: If this is a RICHEDIT control, the control will not
	// send this notification unless you override the CDocument::OnInitDialog()
	// function and call CRichEditCtrl().SetEventMask()
	// with the ENM_CHANGE flag ORed into the mask.
	
	// TODO: Add your control notification handler code here
//	int power=(CMDIFrameWnd*)(theApp.m_pMainWnd)->MDIGetActive()->GetDlgItemInt(IDC_EDIT1);
}*/

void CMdiFFTDoc::OnDeltaposSpin1(NMHDR* pNMHDR, LRESULT* pResult) 
{
	NM_UPDOWN* pNMUpDown = (NM_UPDOWN*)pNMHDR;
	// TODO: Add your control notification handler code here
	
	/*CString s;*/

	m_power-=pNMUpDown->iDelta;		//由Spin改变m_power
	if(m_power>8)m_power=8;			//范围检查
	else if(m_power<1)m_power=1;

	///显示m_power值（在Edit里）
	POSITION pos=GetFirstViewPosition();
	CString s;
	s.Format("%d",m_power);
	((CChildFrame*)(GetNextView(pos)->GetParentFrame()))->ShowPower(s);
//	AfxMessageBox(s);
//SetDlgItemInt(IDC_EDIT2,);
	
	CalcValues();
	UpdateAllViews(NULL);
	
	*pResult = 0;
}