离散傅里叶变换是信号处理、图像处理领域常用的信号处理方法。可以有几种程序实现方法,下面是三种常用的算法实现。
直接依据离散傅里叶变换公式实现的方法最简单,但是,其时间复杂度为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下调试通过。关于傅里叶变换的原理,网上较多,不再赘述。