fftdoc.cpp

数字信号处理的快速傅里叶变换（FFT）代码

// FFTDoc.cpp : implementation of the CFFTDoc class

#include "stdafx.h"
#include "FFT.h"

#include "FFTDoc.h"
#include <math.h>
#include <complex>

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



/////////////////////////////////////////////////////////////////////////////
// CFFTDoc

IMPLEMENT_DYNCREATE(CFFTDoc, CDocument)

BEGIN_MESSAGE_MAP(CFFTDoc, CDocument)
	//{{AFX_MSG_MAP(CFFTDoc)
	ON_COMMAND(ID_INPUT_Xn1, OnINPUTXn1)
	ON_COMMAND(ID_INPUT_Xn2, OnINPUTXn2)
	ON_COMMAND(ID_INPUT_Xn3, OnINPUTXn3)
	ON_COMMAND(ID_INPUT_Xn4, OnINPUTXn4)
	ON_COMMAND(ID_INPUT_Xn5, OnINPUTXn5)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CFFTDoc construction/destruction

CFFTDoc::CFFTDoc()
{
	// TODO: add one-time construction code here
	m_wN = -2.0 * PI / MAXN;
	
	Setxn(64);
	

	FFT();

}

CFFTDoc::~CFFTDoc()
{
}

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

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

	return TRUE;
}



/////////////////////////////////////////////////////////////////////////////
// CFFTDoc serialization

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

/////////////////////////////////////////////////////////////////////////////
// CFFTDoc diagnostics

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

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

/////////////////////////////////////////////////////////////////////////////
// CFFTDoc commands
//序列的快速傅立叶变换
void CFFTDoc::FFT()
{
    UINT i,j,k,N;
	complex<long double> c;

	//清除X(k)中的内容
	for(i=0;i<MAXN;i++)
		ma_FFT[i]=complex<double>(ma_xn[i],0.0);

		////////////////////////////////////////////////////////////////////////
	//时间抽选奇偶分解蝶型算法

	UINT k1,k2;		//存放数组下标的临时变量
	UINT l,			//级数控制、也做群间隔
		 n,			//乘数控制
		 g,			//群数控制
		 group,		//群数
		 mul;		//乘数个数
   double wn;		//第n个乘数的角频率

	for(l=2;l<=MAXN;l*=2)
	{
		group = l/2;   //一级中的群数
		mul   = MAXN/l;     //一级中的乘数个数 
		
		for(n=0;n<mul;n++)//乘数控制
		{
			wn = m_wN * group * n;
			
			for(g=0;g<MAXN;g+=(2*mul))//群数控制
			{
				k1 = n + g;
				k2 = n + g + mul;

				c = complex<long double>(cos(wn),sin(wn));
				ma_FFT[k1] = ma_FFT[k1] +ma_FFT[k2];
				ma_FFT[k2] = (ma_FFT[k1]-ma_FFT[k2]-ma_FFT[k2]) * c;
			}
		}
	}

	/////////////////////////////////////////////////////////////////////////
	//倒位序算法

	
	j=0;				//初始倒位序下标
	N=MAXN/2;			//加快求最高位
	for(i=1;i<MAXN;i++)
	{
		k=N;	//控制下标做右向加法时，最高位
		if(j>=k)
		{
			do
			{
				j=j-k;		//如果当前倒位序下标"不"比当前的最高位 k 小，
								//则当前倒位序下标的最高位应加 k 。
								//由于有进位产生，故以减法完成。
				k=k/2;		//当前的最高位往右移。
			}while(j>=k);	//如果当前倒位序下标"不"比当前的最高位 k 小，
								//继续进行次高位的比较计算。
		}
		j=j+k;		//如果当前倒位序下标比当前的最高位 k 小，
						//说明其最高位是零，故只要直接加 k 即可。
		if(i<j)		//若当前原下标和倒位序下标不同，则交换XK(i)和XK(j)
		{
			c=ma_FFT[i];
			ma_FFT[i]=ma_FFT[j];
			ma_FFT[j]=c;
		}
	}
}

////////////////////////////////////////////
//
//序列快速傅立叶变换的幅频|X|特性`
void CFFTDoc::Aw()
{
	UINT k;
	double dMax=0.0;
	
	if(m_DORF==1)
	{
		for(k=0;k<MAXN;k++)
		{
			ma_Aw[k]=sqrt(ma_FFT[k].real() * ma_FFT[k].real()+
						  ma_FFT[k].imag() * ma_FFT[k].imag());
			if(dMax<=ma_Aw[k])
				dMax = ma_Aw[k];
		}
	}
	//归一化
	for(k=0;k<MAXN;k++)
	{
		ma_Aw[k]=ma_Aw[k]/dMax;
	}
}

////////////////////////////////////////////
//
//序列快速傅立叶变换的相频Arg(X)特性
void CFFTDoc::Pw()
{


	UINT k;
	double dMax=0.0;
	if(m_DORF==1)
	{
		for(k=0;k<MAXN;k++)
		{
			ma_Pw[k] = atan(ma_FFT[k].imag() / ma_FFT[k].real());
			if(dMax<fabs(ma_Pw[k]))
				dMax = fabs(ma_Pw[k]);
		}
	}
	//归一化
	for(k=0;k<MAXN;k++)
	{
		ma_Pw[k] = ma_Pw[k] / dMax;
	}
}

////////////////////////////////////////////
//
//取幅频|X|序列的一个点
double CFFTDoc::GetAw(const UINT k)
{
	return ma_Aw[k];
}

////////////////////////////////////////////
//
//取相频Arg(X)序列的一个点
double CFFTDoc::GetPw(const UINT k)
{
	return ma_Pw[k];
}

////////////////////////////////////////////
//
//取时域序列x(n)的一个点
double CFFTDoc::Getxn(const UINT n)
{
	return ma_xn[n];
}

////////////////////////////////////////////
//
//取时域序列x(n)的当前设置的点数
UINT CFFTDoc::GetN()
{
	return m_N;
}

////////////////////////////////////////////
//
//设置时域序列x(n)的值、取值为"1"的点数
void CFFTDoc::Setxn(const UINT N)
{
	UINT i;

	m_N =  N;
	for(i=0 ; i<m_N ; i++)
	{
		ma_xn[i] = 1.0;
	}
	for(i=m_N ; i<MAXN ; i++)
	{
		ma_xn[i] = 0.0;
	}
}


void CFFTDoc::OnINPUTXn1() 
{
	// TODO: Add your command handler code here
	Setxn(16);
	

	FFT();

	UpdateAllViews(NULL);
}

void CFFTDoc::OnINPUTXn2() 
{
	// TODO: Add your command handler code here
	Setxn(32);
	

	FFT();

	UpdateAllViews(NULL);
	
}

void CFFTDoc::OnINPUTXn3() 
{
	// TODO: Add your command handler code here
	Setxn(64);
	

	FFT();

	UpdateAllViews(NULL);
	
}

void CFFTDoc::OnINPUTXn4() 
{
	// TODO: Add your command handler code here
	Setxn(128);
	

	FFT();

	UpdateAllViews(NULL);
	
}

void CFFTDoc::OnINPUTXn5() 
{
	// TODO: Add your command handler code here
	Setxn(256);
	

	FFT();

	UpdateAllViews(NULL);
	
}