递归版的快速傅里叶转换 (FFT)，C语言实现

Dear 丶 2022-12-01 03:57 98阅读 0赞

看了好多个版本的快速傅里叶转换 (FFT)，看懂了理论部分，但是却无法把理论联系到代码。也许是我的水平太差了吧。网上看到的版本中大量使用位操作符<<和>>。开始的时候还真的反应不过来。这里，我不讲FFT的理论，只是实现一下递归版本的FFT实现。有网友说FFT是玄科学，我深表赞同。相同的代码，在A电脑里运行结果是正确的，用B电脑就错误。找了好久才发现只需要把两个数组申请很大的空间即可，但这是为什么呢？反正没人告诉我原因。下面是具体代码：

#include <stdio.h> 
    #include <stdlib.h>
    #include <math.h>
    
    double PI = 3.1415926;
    
    //定义复数 
    struct Complex{
        
    	double x; //实部 
    	double y; //虚部 
    };
    
    //网上大部分都是采用C++实现FFT，使用了操作符重载。
    //但是其它语言没有操作符重载，所以我们还是用函数的方法实现复数的加法、减法和乘法
    struct Complex add(struct Complex a, struct Complex b)//复数的加法 
    {
        
    	struct Complex c;
    	c.x = a.x + b.x;
    	c.y = a.y + b.y;
    	return c;
    }
    
    //复数的减法，不能用minus，好像是与关键字冲突 
    struct Complex myMinus(struct Complex a, struct Complex b) 
    {
        
    	struct Complex c;
    	c.x = a.x - b.x;
    	c.y = a.y - b.y;
    	return c;
    }
    
    struct Complex multiply(struct Complex a, struct Complex b) //复数的乘法
    {
        
    	struct Complex c;
    	c.x = a.x * b.x - a.y * b.y;
    	c.y = a.x * b.y + b.x * a.y;
    	return c;
    }
    
    void fft(int n, Complex * arr, int type)
    {
        
        if(n == 1) //递归的出口 
        {
        
        	return;
    	}
        
    	int m = n / 2; //由于保证了n是足够大的偶数，所以m一定是n的一半 
    	
    	struct Complex *arr1; //储存arr的偶数部分 
    	arr1 = (struct Complex *)(malloc(m * sizeof(struct Complex))); //动态申请内存 
    	
    	struct Complex *arr2; //储存arr的奇数部分
    	arr2 = (struct Complex *)(malloc(m * sizeof(struct Complex))); //动态申请内存 
    	
        for (int i = 0; i < m; i++)
        {
        
        	arr1[i] = arr[i * 2];     //在arr1中储存arr的偶数部分
    		arr2[i] = arr[i * 2 + 1]; //在arr2中储存arr的奇数部分
    	}
            
        fft(m, arr1, type); //对偶数部分进行快速傅里叶变换 
        fft(m, arr2, type); //对奇数部分进行快速傅里叶变换
        
        struct Complex wn; //把360度等分成n份，这时的复数。 
        wn.x = cos(2.0 * PI / n); //复数的x 
        wn.y = type * sin(2.0 * PI / n); //复数的y，注意逆运算的时候，type为-1 
        
        struct Complex w; //x = 1，y = 0的复数 
        w.x = 1;
        w.y = 0;
        
        //计算arr中的每一个值 
        for (int i = 0; i < m; i++) 
        {
        
            arr[i] =         add(arr1[i], multiply(w, arr2[i]));
            arr[i + m] = myMinus(arr1[i], multiply(w, arr2[i]));
            w = multiply(w, wn);
        }
    }
    
    int main()
    {
        
        int n = 9; //数组a的长度为9 
    	int m = 5; //数组b的长度为5
    	int complexSize = sizeof(struct Complex);
        
        //要申请很大的空间，不知道为什么 
        //如果只乘以100或者更小的空间，就会出错。 
        struct Complex * a = (struct Complex *)(malloc(1000 * complexSize));
    	struct Complex * b = (struct Complex *)(malloc(1000 * complexSize));
    	struct Complex * c = (struct Complex *)(malloc(1000 * complexSize));
    	
    	for(int i = 0; i < n; i++) //数组a的初始化 
    	{
        
    		a[i].x = 1;
    		a[i].y = 0;
    	}
    	
    	for(int i = 0; i < m; i++) //数组b的初始化
    	{
        
    		b[i].x = 1;
    		b[i].y = 0;
    	}
            
        int lim = 1; 
    	
    	//lim为1, 2, 4, 8, 16, 32....除了1之外，其余都是偶数
    	//lim为1时，fft直接就递归返回了，所以实际上lim都是偶数 
        while (lim < n + m - 1)
        {
        
    		lim = lim * 2;
    	}
        
        //即两函数的傅里叶变换的乘积等于它们卷积后的傅里叶变换，
    	//能使傅里叶分析中许多问题的处理得到简化。
    	//F(g(x) * f(x)) = F(g(x)) * F(f(x))
    	//其中F表示的是傅里叶变换
        fft(lim, a, 1); //快速傅里叶转换 
        fft(lim, b, 1); //快速傅里叶转换
        
        for (int i = 0; i < lim; i++)	
        {
        
        	c[i] = multiply(a[i], b[i]); //把a和b中相乘的结果都放到c中 
    	}
            
        fft(lim, c, -1); //对c进行快速傅里叶转换的逆转换 
        
        //取实数四舍五入，此时虚数部分应当为0或由于浮点误差接近0
        //理论上计算结果应该是虚部为0，但实际上虚部可能有一点误差
    	//为了抵消这个误差，用下面的方法进行处理
    	//最终显示c的内容 
        for (int i = 0; i < n + m - 1; i++)
        {
        
        	printf("%d, ", (int)(c[i].x / lim + 0.5));
    	}
        return 0;
    }

运行结果如下：

![在这里插入图片描述][20200825212309493.png_pic_center]  
如果main函数中，struct Complex \* a、b、c申请的空间比较小 (100或更小)，则会出错，不知道为什么。

![在这里插入图片描述][20200825212423336.png_pic_center]  
下次再写迭代版本的FFT吧。

[20200825212309493.png_pic_center]: /images/20221123/0aaef4651a4b40c2ad09840f28267fe6.png
[20200825212423336.png_pic_center]: /images/20221123/a0bf80208e3f43ac86ca95ee3bee9000.png