(本人为本文原作者,转载请标明出处)
上文请转:
DBinary:离散傅里叶变换DFT详解及应用zhuanlan.zhihu.com第一节 DFT的时间复杂度
设:
那么式3.52可以写成
设一个有四个样本的离散序列x[0],x[1],x[2],x[3]那么DFT变换的行列式可以写成:
时间复杂度是一个用于描述算法运行时间的函数,用符号O表示,从DFT的正变换及逆变换公式可以看出,其算法的时间复杂度为
第二节 基2快速傅里叶变换
设对离散序列x[n]的DFT变换写作
其中n的范围为[0,N-1],如果n是一个偶数,那么,就可以将x[n]分为两组,设i为
进一步可以得到
因为
那么,5.7可以写成:
将
于是一个DFT变换变成了两个
从式5.11中可知,频域复信号的前
通过变形得到
因为
那么,5.13就变为了
其与前半部分的变换公式相比仅仅只是多取了一个负号,那么上式就可以写成
即得到总结公式
其中DFT(e[i])与DFT[o[i]]可以根据上述推导进一步分解,,直到需要变换的点数为2,那么,傅里叶变换的时间复杂度就由
第三节 基2快速傅里叶逆变换
根据式
可以写成
根据正变换的推导过程,同样将复信号序列分为奇数序列与偶数序列,得到式:
简化后得到
即得到式子
以此类推可得后半部分式子表示为:
其中,式5.19与式5.10相差部分仅仅是由
至此,快速傅里叶变换的逆变换部分推导完成了(IFFT)。
第六章 用C语言实现FFT、IFFT代码
第一节 用C语言实现FFT代码
由第五章节的推导可以知道,基2 快速傅里叶变换实际上是DFT代码的不断分解计算,直到最终分解为2点DFT。因此,可以用迭代的方式实现FFT代码:
根据公式5.17及欧拉公式2.32,将傅里叶变换写为易于编码的形式
void FFT_Sort(complex x[],int N)
{
//定义一些临时变量,其中er,eb用于排序,i,j,k就是公式里对应的i,j,k
int eb,er,i,j,k;
//定义需要的虚数变量
complex excomplex,Wnk,cx0,cx1;
if (N>>2)//如果样本数大于2
{
/* 将x[]分为奇偶两部分序列并对原序列重新排序,其中将偶数部分放到序列的前端将奇数部分放到序列的末端*/
eb=1;
er=0;
while (er<N)
{
er=eb<<2;
for (i=0;i<N/er;i++)
{
for (j=0;j<eb;j++)
{
excomplex=x[er*i+eb+j];
x[er*i+eb+j]=x[er*i+eb*2+j];
x[er*i+eb*2+j]=excomplex;
}
}
eb<<=1;
}
//迭代计算两部份的DFT
FFT_Sort(x,N>>1);//计算前半部分的偶数列
FFT_Sort(x+(N>>1),N>>1);//计算后半部分奇序列
//根据公式后半部分需要乘以一个wnk,其中__PI是一个宏定义,为3.14159265f
//计算完成后,将两个部分合并起来
for(k=0;k<N>>1;k++)
{
Wnk.re=(float)cos(-2*__PI*k/N);
Wnk.im=(float)sin(-2*__PI*k/N);
cx0=x[k];
cx1=x[k+(N>>1)];
x[k]=complexAdd(cx0,complexMult(Wnk,cx1));
Wnk.re=-Wnk.re;
Wnk.im=-Wnk.im;
x[k+(N>>1)]=complexAdd(cx0,complexMult(Wnk,cx1));
}
}
else
{
//一直分解直到样本点只有2个,那么就是两点DFT了,直接硬编码计算
cx0=x[0];
cx1=x[1];
x[0]=complexAdd(cx0,cx1);
cx1.im=-cx1.im;
cx1.re=-cx1.re;
x[1]=complexAdd(cx0,cx1);
}
}
编写FFT的迭代代码:
其中void FFT_Sort(complex x[],int N)是FFT的迭代函数,其中x[]表示该函数的输入复数信号,同时也作为其输出复信号,N为x[]的样本个数,该函数实现的基本思路是,如果输入的样本数为2,就直接进行2点的DFT变换.否者对x[]进行排序,将偶序列依次排列到该序列的前半部分,同时将奇序列排序到该序列的后半部分,例如如下序列
排序前:1,2,3,4,5,6,7,8,9
排序后:2,4,6,8,1,3,5,7
然后通过迭代调用自身,对前半部分别进行DFT变换,最终就得到了变换的结果为了避免函数对输入原信号的直接修改,编写另一个函数拷贝原信号再进行进一步的傅里叶运算。
void FFT (complex x[],complex X[],int N)
{
int size=1;
complex *cplx;
//计算一个2基数,以用于快速傅里叶变换
while((size<<=1)<N);
//分配新内存用于拷贝原信号
cplx=(complex *)malloc(sizeof(complex)*size);
//拷贝原信号到该内存
memset(cplx,0,sizeof(complex)*N);
memcpy(cplx,x,sizeof(complex)*N);
//进行快速傅里叶变换
FFT_Sort(cplx,size);
//将结果拷贝到输出缓存中
memcpy(X,cplx,sizeof(complex)*N);
//释放分配的内存空间
free(cplx);
}
第二节 用C语言实现IFFT代码
傅里叶逆变换的使用的公式如下
其变换代码可以由第一节的IFFT修改而成:void IFFT_Base2(_IN _OUT complex X[],int N)
{
//定义一些临时变量,其中er,eb用于排序,i,j,k就是公式里对应的i,j,k
int eb,er,i,j,n;
//定义需要的虚数变量
complex excomplex,Wnnk,cx0,cx1;
if (N>>2)
{
/* 将X[]分为奇偶两部分序列并对原序列重新排序,其中将偶数部分放到序列的前端将奇数部分放到序列的末端*/
eb=1;
er=0;
while (er<N)
{
er=eb<<2;
for (i=0;i<N/er;i++)
{
for (j=0;j<eb;j++)
{
excomplex=X[er*i+eb+j];
X[er*i+eb+j]=X[er*i+eb*2+j];
X[er*i+eb*2+j]=excomplex;
}
}
eb<<=1;
}
//迭代计算两部份的IDFT
IFFT_Base2(X,N>>1);
IFFT_Base2(X+(N>>1),N>>1);
//根据公式后半部分需要乘以一个wnk,其中__PI是一个宏定义,为3.14159265f
//计算完成后,将两个部分合并起来
for(n=0;n<N>>1;n++)
{
/*这里与FFT变换的代码有所不同,cos sin值不需要加负号*/
Wnnk.re=(float)cos(2*__PI*n/N);
Wnnk.im=(float)sin(2*__PI*n/N);
cx0=X[n];
cx1=X[n+(N>>1)];
X[n]=complexAdd(cx0,complexMult(Wnnk,cx1));
Wnnk.re=-Wnnk.re;
Wnnk.im=-Wnnk.im;
X[n+(N>>1)]=complexAdd(cx0,complexMult(Wnnk,cx1));
}
}
else
{
//2点IDFT计算代码
cx0=X[0];
cx1=X[1];
X[0]=complexAdd(cx0,cx1);
cx1.im=-cx1.im;
cx1.re=-cx1.re;
X[1]=complexAdd(cx0,cx1);
}
}
同样为了避免函数对输入原信号的直接修改,编写另一个函数拷贝原信号再进行进一步的傅里叶逆变换运算。
void IFFT(_IN complex X[],_OUT complex x[],int N)
{
int size=1,i;
complex *i_px;
//计算一个2基数,以用于快速傅里叶逆变换
while((size<<=1)<N);
//分配一个新的内存空间用于拷贝原信号
i_px=(complex *)malloc(sizeof(complex)*size);
//拷贝原信号到该内存
memset(i_px,0,sizeof(complex)*N);
memcpy(i_px,X,sizeof(complex)*N);
//快速傅里叶逆变换
IFFT_Base2(i_px,size);
//与快速傅里叶变换不同的是,逆变换每个信号需要除以一个N
for (i=0;i<N;i++)
{
i_px[i].re/=N;
i_px[i].im/=N;
}
//将计算结果拷贝到输出缓存中
memcpy(x,i_px,sizeof(complex)*N);
free(i_px);
}
至此,快速傅里叶逆变换的代码也编写完成了。
第三节 验证FFT/IFFT代码
快速傅里叶变换对(FFT/IFFT)的本质仍然是DFT/IDFT的计算方式,只是在算法方面进行了优化让其更加的适合于计算机运算,因此FFT/IFFT结果应该和DFT/IDFT的结果相一致,修改第四章第二三节的代码,验证其结果是否一致以判断代码编写是否正确。
int main()
{
//构造一个长度为8的复信号样本
complex samples[8],_out[8];
int i;
//初始化这个复信号样本序列为1,2,3,4,5,6,7,8
for (i=0;i<8;i++)
{
samples[i].re=i;
samples[i].im=0;
}
//对该信号进行DFT运算
printf("DFT:n");
//输出运算后的结果
DFT(samples,_out,8);
for (i=0;i<8;i++)
{
printf("(%f,%f)n",_out[i].re,_out[i].im);
}
//对该信号进行IDFT逆变换
printf("IDFT:n");
IDFT(_out,samples,8);
//输出结果,应该与原信号相一致
for (i=0;i<8;i++)
{
printf("(%f,%f)n",samples[i].re,samples[i].im);
}
//重新初始化原信号(因为经过逆变换存在精度误差)
for (i=0;i<8;i++)
{
samples[i].re=i;
samples[i].im=0;
}
//进行FFT变换,应该与DFT变换的结果一致
printf("FFT:n");
FFT(samples,_out,8);
//输出结果
for (i=0;i<8;i++)
{
printf("(%f,%f)n",_out[i].re,_out[i].im);
}
//进行IFFT逆变换,应能够还原出原信号
printf("IFFT:n");
IFFT(_out,samples,8);
//输出结果
for (i=0;i<8;i++)
{
printf("(%f,%f)n",samples[i].re,samples[i].im);
}
system("pause");
}
该范例程序(DFT_DEMO)的运行结果显示如图6-2-1所示
从结果来看,FFT/IFFT的运算结果与DFT/IDFT的结果一致,甚至FFT/IFFT的结果比DFT/IDFT的结果在精度损失上更小(这是因为FFT/IFFT在算法上的浮点运算量小于DFT/IDFT因此其因为浮点精度导致的误差更小),至此快速傅里叶变换的功能代码编写结束。