dit-fft.cpp

本人编写的按时间抽取的快速傅立叶变换算法

/*---------------------------------------------------------------------
 Routine fft:to obtain the DFT of Complex Data x(n)
                               By Cooley-Tukey radix-2 DIT Algorithm .
 input parameters:
 x : complex array.input signal is stored in x(0) to x(n-1).
 n : the dimension of x and y.
 isign: if  isign=-1 For Forward Transform
        if  isign=+1 For Inverse Transform.
 output parameters:
 x : complex array. DFT result is stored in x(0) to x(n-1).
 Notes:
     n must be power of 2.
---------------------------------------------------------------------*/

#include<iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct {double real,imag;} complex;

complex cexp(complex a)
{
 complex z;
 z.real=exp(a.real)*cos(a.imag);
 z.imag=exp(a.real)*sin(a.imag);
 return(z);
 }

void mcmpfft(complex x[],int n,int isign)
{
	complex t,z,ce;
	double pisign;
	int mr,m,l,j,i,nn;
	for(i=1;i<=16;i++)    //n must be power of 2
	   {
		nn=(int)pow(2,i);
		if(n==nn) break;
	    }
	if(i>16)
	{
		printf(" N is not a power of 2 ! \n");
		return;
	}
	z.real=0.0;    
	pisign=4*isign*atan(1.);  //pisign的值为+180度或-180度

 //码位倒置
	mr=0;          
	for(m=1;m<n;m++)
	{l=n;
	while(mr+l>=n)   l=l/2;  
	mr=mr%l+l;   
	if(mr<=m)  	continue;
	t.real=x[m].real;t.imag=x[m].imag;
	x[m].real=x[mr].real;x[m].imag=x[mr].imag;
	x[mr].real=t.real;x[mr].imag=t.imag;
	}

	l=1;
	while(1)
	{
	  if(l>=n)     
	  {
		  if(isign==-1)    //isign=-1 For Forward Transform
		  return;
	       for(j=0;j<n;j++)    //  Inverse Transform
		   {
			   x[j].real=x[j].real/n;
			   x[j].imag=x[j].imag/n;
		    }
	       return;
	    }
	  for(m=0;m<l;m++)  //完成当前级所有蝶形运算
	  {
		for(i=m;i<n;i=i+2*l)//完成当前级相同W因子的所有蝶形运算
		 {
		     z.imag=m*pisign/l;  
		     ce=cexp(z);  
		     t.real=x[i+l].real*ce.real-x[i+l].imag*ce.imag;
		     t.imag=x[i+l].real*ce.imag+x[i+l].imag*ce.real;
		     x[i+l].real=x[i].real-t.real;   //原位运算
		     x[i+l].imag=x[i].imag-t.imag;
		     x[i].real=x[i].real+t.real;
		     x[i].imag=x[i].imag+t.imag;
                   }
	       }
	  l=2*l;  //确定下一级蝶形运算中W因子个数
	  }
}

void main()
{
	complex x[32];
	int n,i;
	int isign=0;
	while(isign!=1&&isign!=-1)
	{
		cout<<"Input : '-1'---Forward Transform;'+1'---Inverse Transform\n";
		cin>>isign;
	}
	cout<<"Input the dimension of x:  n=";
	cin>>n;
	cout<<"\nInput complex array:\n";
	for(i=0;i<n;i++)     
	{  
		cin>>x[i].real>>x[i].imag;
		cout<<"x["<<i<<"]="<<x[i].real<<"+j"<<x[i].imag<<endl;
			}
	mcmpfft(x,n,isign);
	cout<<"\nThe result:\n";
	for(i=0;i<n;i++)
	{
		cout<<x[i].real<<"+j"<<x[i].imag<<endl;
	}
}