傅里叶变换公式表_快速傅里叶变换推导

(本人为本文原作者,转载请标明出处)

上文请转:

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为

,设满足x[2i]的序列为偶序列,设满足x[2i+1]的序列为奇数列,奇数列和偶数列的个数分别是

个,那么,式5.4可以写成

进一步可以得到

因为

那么,5.7可以写成:

将

用m进行替换,x[2i]用e[i]表示,意为偶序列(even),将x[2i+1]用o[i]表示,意为奇序列(odd),那么式5.9可以写成:

于是一个DFT变换变成了两个

从式5.11中可知,频域复信号的前

可以分解为两个傅里叶变换,下面对后半部分进行推导,设傅里叶变换后半部分的变换写为

通过变形得到

因为

那么,5.13就变为了

其与前半部分的变换公式相比仅仅只是多取了一个负号,那么上式就可以写成

即得到总结公式

其中DFT(e[i])与DFT[o[i]]可以根据上述推导进一步分解,,直到需要变换的点数为2,那么,傅里叶变换的时间复杂度就由

降低到了

,当然,其需要样本的点数满足2的整数次幂,同时也将这种简化方式称为基2快速傅里叶变换(FFT)。

第三节 基2快速傅里叶逆变换

根据式

可以写成

根据正变换的推导过程，同样将复信号序列分为奇数序列与偶数序列，得到式：

简化后得到

即得到式子

以此类推可得后半部分式子表示为：

其中，式5.19与式5.10相差部分仅仅是由

变为了

，那么根据第二节的推导过程，式5.19可以写成：

至此，快速傅里叶变换的逆变换部分推导完成了（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所示

图6-2-1 FFT/IFFT变换结果

从结果来看,FFT/IFFT的运算结果与DFT/IDFT的结果一致,甚至FFT/IFFT的结果比DFT/IDFT的结果在精度损失上更小(这是因为FFT/IFFT在算法上的浮点运算量小于DFT/IDFT因此其因为浮点精度导致的误差更小)，至此快速傅里叶变换的功能代码编写结束。