fftdoc.cpp

FFT算法进行ＦＦＴ 、ＩＦＦＴ、功率谱计算

// FFTDoc.cpp : implementation of the CFFTDoc class
//

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

#include "FFTDoc.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#define NUM	1024
#define FRATE 0.4707416
#define HRATE 22.5
/////////////////////////////////////////////////////////////////////////////
// CFFTDoc

IMPLEMENT_DYNCREATE(CFFTDoc, CDocument)

BEGIN_MESSAGE_MAP(CFFTDoc, CDocument)
	//{{AFX_MSG_MAP(CFFTDoc)
	ON_COMMAND(ID_FILE_OPEN, OnFileOpen)
	//}}AFX_MSG_MAP
//	ON_MESSAGE(WM_MSG_WAVE, OnMsgWave)
END_MESSAGE_MAP()

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

CFFTDoc::CFFTDoc()
{
   N = 256;           /* FFT 点数 */
   M = 512;           /* ZFFT 点数 */
     if(N>M) {
     tamp = (double *) calloc(N+1,sizeof(double));
	 famp = (double *) calloc(N+1,sizeof(double));
     samp = (COMPLEX *) calloc(N+1,sizeof(COMPLEX));
     rsamp = (COMPLEX *) calloc(N+1,sizeof(COMPLEX));
  }
  else {
     tamp = (double *) calloc(M+1,sizeof(double));
	 famp = (double *) calloc(M+1,sizeof(double));
     samp = (COMPLEX *) calloc(M+1,sizeof(COMPLEX));
     rsamp = (COMPLEX *) calloc(M+1,sizeof(COMPLEX));
  }
  for (i = 0; i < N; i++)  {
     samp[i].real = 0;
     samp[i].imag = 0; 
	 rsamp[i].real = 0;
     rsamp[i].imag = 0;
	 tamp[i]=0;
	 famp[i]=0;
		}
	tym=0;fym=0;
    m_bPower=false;
	numav=16;
	slice=N;
	ovlap=128;
	wtype=0;
	ntype=0;
	estsize = (slice-ovlap)*(numav-1) + slice;
	sigfloat = (double *) calloc(estsize,sizeof(double));
    mag = (double *) calloc(slice,sizeof(double));
	m_fop=false;
	m_chigh=0;
	m_cfreq=0;
}

CFFTDoc::~CFFTDoc()
{  free(samp);
   free(rsamp);
   free(tamp); free(famp); 
   free(mag);
   free(sigfloat);
}

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
/************************************************************************
fft - 基2 DIF FFT 子程序

输入参数:
	 COMPLEX *x : FFT 输入和输出数据区指针;
	      int m : FFT 长度 ( length = 2^m );
输出参数:
	 输出数据放在 x 所指的输入数据区.
	 无输出参数.

void fft(COMPLEX *x, int m)
*************************************************************************/

void CFFTDoc::fft(COMPLEX *x, int m)
{

  static COMPLEX *w;           /* used to store the w complex array */
    static int mstore = 0;       /* stores m for future reference */
    static int n = 1;            /* length of fft stored for future */

    COMPLEX u,temp,tm;
    COMPLEX *xi,*xip,*xj,*wptr;

    int i,j,k,l,le,windex;

    double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;

    if(m != mstore) {

/* free previously allocated storage and set new m */

	if(mstore != 0) free(w);
	mstore = m;
	if(m == 0) return;       /* if m=0 then done */

/* n = 2**m = fft length */

	n = 1 << m;
	le = n/2;

/* allocate the storage for w */

	w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
	if(!w) {
	    printf("\nUnable to allocate complex W array\n");
	    exit(1);
	}

/* calculate the w values recursively */

	arg = 4.0*atan(1.0)/le;         /* PI/le calculation */
	wrecur_real = w_real = cos(arg);
	wrecur_imag = w_imag = -sin(arg);
	xj = w;
	for (j = 1 ; j < le ; j++) {
	    xj->real = (float)wrecur_real;
	    xj->imag = (float)wrecur_imag;
	    xj++;
	    wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
	    wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
	    wrecur_real = wtemp_real;
	}
    }

/* start fft */

    le = n;
    windex = 1;
    for (l = 0 ; l < m ; l++) {
	le = le/2;

/* first iteration with no multiplies */

	for(i = 0 ; i < n ; i = i + 2*le) {
	    xi = x + i;
	    xip = xi + le;
	    temp.real = xi->real + xip->real;
	    temp.imag = xi->imag + xip->imag;
	    xip->real = xi->real - xip->real;
	    xip->imag = xi->imag - xip->imag;
	    *xi = temp;
	}

/* remaining iterations use stored w */

	wptr = w + windex - 1;
	for (j = 1 ; j < le ; j++) {
	    u = *wptr;
	    for (i = j ; i < n ; i = i + 2*le) {
		xi = x + i;
		xip = xi + le;
		temp.real = xi->real + xip->real;
		temp.imag = xi->imag + xip->imag;
		tm.real = xi->real - xip->real;
		tm.imag = xi->imag - xip->imag;
		xip->real = tm.real*u.real - tm.imag*u.imag;
		xip->imag = tm.real*u.imag + tm.imag*u.real;
		*xi = temp;
	    }
	    wptr = wptr + windex;
	}
	windex = 2*windex;
    }

/* rearrange data by bit reversing */

    j = 0;
    for (i = 1 ; i < (n-1) ; i++) {
	k = n/2;
	while(k <= j) {
	    j = j - k;
	    k = k/2;
	}
	j = j + k;
	if (i < j) {
	    xi = x + i;
	    xj = x + j;
	    temp = *xj;
	    *xj = *xi;
	    *xi = temp;
	}
    }
}
/************************************************************************
ifft - 基2 DIF IFFT 子程序

输入参数:
	 COMPLEX *x : IFFT 输入和输出数据区指针;
	      int m : IFFT 长度 ( length = 2^m );
输出参数:
	 输出数据放在 x 所指的输入数据区.
	 无输出参数.

void ifft(COMPLEX *x, int m)
*************************************************************************/

void CFFTDoc::ifft(COMPLEX *x, int m)
{
    static COMPLEX *w;           /* used to store the w complex array */
    static int mstore = 0;       /* stores m for future reference */
    static int n = 1;            /* length of ifft stored for future */

    COMPLEX u,temp,tm;
    COMPLEX *xi,*xip,*xj,*wptr;

    int i,j,k,l,le,windex;

    double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
    float scale;

    if(m != mstore) {

/* free previously allocated storage and set new m */

	if(mstore != 0) free(w);
	mstore = m;
	if(m == 0) return;       /* if m=0 then done */

/* n = 2**m = inverse fft length */

	n = 1 << m;
	le = n/2;

/* allocate the storage for w */

	w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
	if(!w) {
	    printf("\nUnable to allocate complex W array\n");
	    exit(1);
	}

/* calculate the w values recursively */

	arg = 4.0*atan(1.0)/le;         /* PI/le calculation */
	wrecur_real = w_real = cos(arg);
	wrecur_imag = w_imag = sin(arg);  /* opposite sign from fft */
	xj = w;
	for (j = 1 ; j < le ; j++) {
	    xj->real = (float)wrecur_real;
	    xj->imag = (float)wrecur_imag;
	    xj++;
	    wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
	    wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
	    wrecur_real = wtemp_real;
	}
    }

/* start inverse fft */

    le = n;
    windex = 1;
    for (l = 0 ; l < m ; l++) {
	le = le/2;

/* first iteration with no multiplies */

	for(i = 0 ; i < n ; i = i + 2*le) {
	    xi = x + i;
	    xip = xi + le;
	    temp.real = xi->real + xip->real;
	    temp.imag = xi->imag + xip->imag;
	    xip->real = xi->real - xip->real;
	    xip->imag = xi->imag - xip->imag;
	    *xi = temp;
	}

/* remaining iterations use stored w */

	wptr = w + windex - 1;
	for (j = 1 ; j < le ; j++) {
	    u = *wptr;
	    for (i = j ; i < n ; i = i + 2*le) {
		xi = x + i;
		xip = xi + le;
		temp.real = xi->real + xip->real;
		temp.imag = xi->imag + xip->imag;
		tm.real = xi->real - xip->real;
		tm.imag = xi->imag - xip->imag;
		xip->real = tm.real*u.real - tm.imag*u.imag;
		xip->imag = tm.real*u.imag + tm.imag*u.real;
		*xi = temp;
	    }
	    wptr = wptr + windex;
	}
	windex = 2*windex;
    }

/* rearrange data by bit reversing */

    j = 0;
    for (i = 1 ; i < (n-1) ; i++) {
	k = n/2;
	while(k <= j) {
	    j = j - k;
	    k = k/2;
	}
	j = j + k;
	if (i < j) {
	    xi = x + i;
	    xj = x + j;
	    temp = *xj;
	    *xj = *xi;
	    *xi = temp;
	}
    }

/* scale all results by 1/n */
    scale = (float)(1.0/n);
    for(i = 0 ; i < n ; i++) {
	x->real = scale*x->real;
	x->imag = scale*x->imag;
	x++;
    }
}

void CFFTDoc::SinWave()
{    InitArray();
	for(i=0;i<N;i++)
		{
			samp[i].real=5*sin((double)(i)/N*PI*6);//设置频率
		}
    CallFFT();
    UpdateAllViews(NULL);
//	BOOL  sendresult=PostMessage(WM_MSG_WAVE,0,0);
			
}