傅里叶变换算法c语言实现（递归与非递归）

1、算法设计原理分析

 

把傅里叶变换公式进行拆分，按奇偶进行分组，使之转外为两个T(n/2)的计算式。

这样就可以找到了递归特性，可以一直拆分下去，直到一组只有一个数值。递归算法实现则可以为在每次递归中对偶数组进行再次划分进行递归，同样对奇数组进行再次划分进行递归。

其次利用傅里叶变化的周期性为N/2和对称性，可以推导出公式X(k)=G(k)+W*H(k),和X(K+n/2)= G(k)-W*H(k) ，所以每次只需要计算一半的计算量。其中系数W可以采用欧拉公式展开计算。

非递归实现则可以参考蝶形算法，由图弄清楚算法具体实现过程。

总共进行了log2(N)次蝶形，作为最外层循环条件，重点在于需要在进行迭代之前，对原数组进行重新排序，得到递归法最底层的排序结果，从而就可以从重排后的数组开始计算，第二层循环条件为次蝶形层中的奇偶组对数，最里层的循环条件为当前层的n值。

2、程序设计

递归实现：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.1415926535
typedef struct Complex               //定义复数结构体
{
    double real;
    double image;
} COMPLEX;

void FastFT(COMPLEX buf[], COMPLEX out[], int N, int step)
{
    int i;
    COMPLEX t;
    double a, b, c, d, w;
    if (step < N)              //每递归一次，step*2,递归层数为log2(N)
    {
        FastFT(out, buf, N, step * 2);
        FastFT(out + step, buf + step, N, step * 2);
        for (i = 0; i < N; i += 2 * step)
        {
            w = -PI * i / N;  //系数的指数值
            a = cos(w);
            b = sin(w);    //欧拉公式转换
            c = out[i + step].real; 
            d = out[i + step].image;  //c,d表示奇数组的第i个数的实数部分和虚数部分
            t.real = a * c - b * d;   //权重系数乘以奇数组数值
            t.image = b * c + a * d;  //t表示奇数组的第i个数乘以系数后的值
            buf[i / 2].real = out[i].real + t.real;
            buf[i / 2].image = out[i].image + t.image; //求出递归返回后d对应位置上的值
            buf[(i + N) / 2].real = out[i].real - t.real;
            buf[(i + N) / 2].image = out[i].image - t.image; //由傅里叶的对称性可以求出后一半对应位置的值
        }
    }
}


int fft(COMPLEX buf[], int n)
{
    int i;
    COMPLEX *out = (COMPLEX *)malloc(n * sizeof(COMPLEX));  //开辟一个和原数组等大的数组，用于存储计算奇数组的值的变化
    if (out == NULL)
    {
        return 0;
    }
    for (i = 0; i < n; i++)
    {
        out[i].real = buf[i].real;
        out[i].image = buf[i].image;
    }
    FastFT(buf, out, n, 1);
    free(out);
    out = NULL;
    return 1;
}


int main()
{
    #define COUNT 8
    COMPLEX buf[COUNT] = { {1.0, 0.0}, {2.0, 0.0}, {3.0, 0.0}, {4.0, 0.0},
                   {5.0, 0.0}, {6.0, 0.0}, {7.0, 0.0} ,{8.0, 0.0} };
    fft(buf, COUNT);
    for(int i=0; i<COUNT;i++)
        printf("(%.4f,%.4fi)\n", buf[i].real, buf[i].image);
    return 0;
}

非递归方式实现（蝶形算法）：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.1415926535

typedef struct Complex               //定义复数结构体 
{
    double real;
    double image;
} COMPLEX;


void RaderReverse(COMPLEX*arr, int N) //定义对数组进行重排的函数 
{     
	int i, j;
	COMPLEX temp;
	//第一个和最后一个数位置不变，故不处理
	for(i=1,j=N/2; i<N/2; i+=2,j+=2) 
	{
		temp.real = arr[j].real;
		temp.image = arr[j].image;
		arr[j].real = arr[i].real;
		arr[j].image = arr[j].image;
		arr[i].real = temp.real;
		arr[i].image = temp.image;
	}
}


void fft(COMPLEX *x,int n) {
	int i=0,j=0,k=0,l=0;
	COMPLEX up,down,product;
	double a,b,c,d,w; //w是系数 
	RaderReverse(x,n); //重排找到蝶形算法最底层 

	for(i=0; i<log(n)/log(2); ++i)  /*log(n)/log(2) 级蝶形运算 stage */  
	{
		l = 1<<i;   //当l前层的n/2值 
		for(j=0;j<n;j+= 2*l)    //当前层总共有n/(2*l)组奇偶对              /*一组蝶形运算 group,每组group的蝶形因子乘数不同*/ 
		{
			for(k=0;k<l;++k)   //k用于每对奇偶对的计算  /*一个蝶形运算 每个group内的蝶形运算的蝶形因子乘数成规律变化*/  
			{
				w=-PI*k/l;    //系数w=-2*PI*k/(2*l)
				a = cos(w);
   				b = sin(w);    //欧拉公式转换 
   				c = x[j+k+l].real;  
   				d = x[j+k+l].image;  //c,d表示奇数组的第i个数的实数部分和虚数部分 
   				product.real=a * c - b * d; 
   				product.image=b * c + a * d;
				up.real= x[j+k].real + product.real;
				up.image= x[j+k].image + product.image;
				down.real= x[j+k].real - product.real;
				down.image= x[j+k].image - product.image;
				
				x[j+k].real= up.real;
				x[j+k].image= up.image;
				x[j+k+l].real= down.real;
				x[j+k+l].image= down.image;

			}
		}
	}
}
int main()
{
   #define COUNT 8
   COMPLEX buf[COUNT] = { {1.0, 0.0}, {2.0, 0.0}, {3.0, 0.0}, {4.0, 0.0},
                   {5.0, 0.0}, {6.0, 0.0}, {7.0, 0.0} ,{8.0, 0.0} }; 
   fft(buf, COUNT);
   for(int i=0; i<COUNT;i++)
		printf("(%.4f,%.4fi)\n", buf[i].real, buf[i].image);
   return 0;
}

3、运行结果

注意：变换数组长度需要为2的整数，这样才可以使分层恰好合适。

 

4、算法分析（复杂度）

可以看出，傅里叶变换约N2次乘法和N2次加法。当N较大时，这个计算量是很大的。利用WN的对称性和周期性，将N点傅里叶变换N／2点的傅里叶变换，／2点DFT总的计算量只是原来的一半，即(N／2)2+(N／2)2=N2／2，这样可以继续分解下去，将N／2再分解为N／4点傅里叶变换等。对于N=偶数点的傅里叶变换都可以分解为2点的傅里叶变换，这样其计算量可以减少为(N／2)log2(N)次乘法和Nlog2(N)次加法。所以递归傅里叶变换的时间复杂度为O(Nlog2(N)),非递归实现的乘法与加法次数也是一样的，所以非递归实现的傅里叶变换的时间复杂度也为O(Nlog2(N))，递归与非递归实现不同之处，在于递归开辟栈来存储递归所需空间，空间复杂度大于非递归实现方式。

感悟：

傅里叶变换算法真的比较难理解，学习这个算法也花了比较就的时间去理解消化，通过傅里叶变换，可以实现时域频域之间的转换。递归方式实现傅里叶变换需明确分组方式是按奇偶分组，计算的核心在于应用傅里叶变换的周期性和对称性，利用推导公式X(k)=G(k)+W*H(k),和X(K+n/2)= G(k)-W*H(k)进行相应的计算，在写递归代码时，开始我认为每次递归都需要重新保存每层的计算结果X（k），尝试了许久后，我发现我的判断、循环条件总有问题，奇偶数组上的数值没对应上，后来我打算不改变数值在原数组中的位置，通过改变部长step来实现奇偶的分组，计算结果保存到相应的位置上。

非递归方式实现，结合蝶形算法的图来理解，就比较容易理解，关键在于如何得到蝶形算法的最低层的数组排序。