离散傅里叶变换的算法实现

离散傅里叶变换是信号处理、图像处理领域常用的信号处理方法。可以有几种程序实现方法，下面是三种常用的算法实现。
直接依据离散傅里叶变换公式实现的方法最简单，但是，其时间复杂度为O(n2)，运算速度较慢，程序如下：
先定义了一种复数类型的结构，及几个基本运算。
typedef struct {
double re;
double im;
} CNumber,*PCNumber;

CNumber Add(CNumber c1,CNumber c2)
{
CNumber c;
c.re=c1.re+c2.re;
c.im=c1.im+c2.im;
return c;
}

CNumber Sub(CNumber c1,CNumber c2)
{
CNumber c;
c.re=c1.re-c2.re;
c.im=c1.im-c2.im;
return c;
}

CNumber Multi(CNumber c1,CNumber c2)
{
CNumber c;
c.re=c1.rec2.re-c1.imc2.im;
c.im=c1.rec2.im+c2.rec1.im;
return c;
}
CNumber Div(CNumber c,int n)
{
CNumber tmp;
tmp.re=c.re/n;tmp.im=c.im/n;
return tmp;
}
//下面是直接依据离散变换公式的程序实现。
void DFT(CNumber t,CNumber f,int r)
{
int i,j;
double angle;
int count=1<<r;
CNumber w;
for (i=0;i<count;i++)
for (j=0;j<count;j++) {
w.re=cos(-2piji/count);
w.im=sin(-2piji/count);
f[i]=Add(f[i],Multi(t[j],w));
}
}
//离散傅里叶反变换的实现。
void IDFT(CNumber t,CNumber f,int r)
{
int i,j;
double angle;
int count=1<<r;
CNumber w;
CNumber tmp;
for (i=0;i<count;i++)
{t[i].im=0;t[i].re=0;}
for (i=0;i<count;i++)
{
for (j=0;j<count;j++)
{
w.re=cos(2piji/count);
w.im=sin(2piji/count);
tmp=Multi(f[j],w);
t[i]=Add(t[i],tmp);
}
t[i]=Div(t[i],count);
}
}
快速傅里叶变换（FFT)，可以极大的提高变换效率。其时间复杂度为O(nlogn)，下面是基于递归方式，代码非常简洁，但是需要较多的内存，实现的C程序代码：
参数说明：t是时域或频域幅值数组，offset是第一个数据偏移，数据长度是2的r次幂，bReverse为真：正变换，假：逆变换。
必须注意的是：逆变换的结果需要除以2的r次幂，才是时域数据。

void FFT_Recursive(CNumber *t,int offset,int r,BOOL bReverse)
{
CNumber a,b,w;
int i,j,k;
if (r==0)
return;
FFT_Recursive(t,offset,r-1,bReverse);
FFT_Recursive(t,offset+(1<<(r-1)),r-1,bReverse);

for (i=0;i<(1<<(r-1));i++)
{   
	if (bReverse)
	{w.re=cos(2*pi/(1<<r)*i);w.im=sin(2*pi/(1<<r)*i);}
	else
	{w.re=cos(-2*pi/(1<<r)*i);w.im=sin(-2*pi/(1<<r)*i);}
	a=Add(t[offset+i],Multi(t[offset+i+(1<<(r-1))],w));
	b=Sub(t[offset+i],Multi(t[offset+i+(1<<(r-1))],w));
	t[offset+i]=a;t[offset+i+(1<<(r-1))]=b;
}

}
注意在变换前，应该将数据数组，"逆序"排序，其实现代码：
void Inverse(PCNumber t,int r)
{
int i,j,k;
PCNumber pTmp;
size_t s=sizeof(CNumber)(1<<r);
pTmp=(PCNumber)malloc(s);
for (j=0;j<(1<<r);j++)
{
k=0;
for (i=0;i<r;i++)
if (j&(1<<i)) k+=1<<(r-i-1);
pTmp[j]=t[k];
}
memcpy(t,pTmp,s);
delete pTmp;
}
最后，我实现了一种非递归的FFT，代码如下，可以使用同一个函数同时完成傅里叶逆变换，时间复杂度为O(nlogn)。
void FFT(CNumber t[],int r,BOOL bInverse)
{
int count;
int i,j,k,m,n,e;
CNumber a,b;
PCNumber w;
double angle;
count=1<<r;
PCNumber pTmp;
pTmp=(PCNumber)malloc(sizeof(CNumber)count);
for (j=0;j<count;j++)
{
k=0;
for (i=0;i<r;i++)
if (j&(1<<i)) k+=1<<(r-i-1);
pTmp[j]=t[k];
}
memcpy(t,pTmp,sizeof(CNumber)count);
free(pTmp);
w=(PCNumber)malloc(sizeof(CNumber)count/2);
for (i=0;i<count/2;i++)
{
if (!bInverse)
angle=-ipi2/count;
else
angle=ipi2/count;
w[i].re=cos(angle);
w[i].im=sin(angle);
}

for(k=0;k<r;k++)
   for(j=0;j<(1<<(r-k))/2;j++)
	   for (i=0;i<(1<<(k+1))/2;i++)
	   {
		  m=j*(2<<k)+i;
		  n=m+(1<<k);
		  e=i*(1<<(r-k-1));
		  a=Add(t[m],Multi(t[n],w[e]));
		  b=Sub(t[m],Multi(t[n],w[e]));
		  t[m]=a;
		  t[n]=b;
	   }
if (bInverse)
	for (i=0;i<count;i++)
		t[i]=Div(t[i],1<<r);
free(w);		

}
花了两天时间，终于完成了上述代码，均在vs 2012下调试通过。关于傅里叶变换的原理，网上较多，不再赘述。