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

骑猪看日落 2022-12-06 12:42 88阅读 0赞

这是迭代版的快速傅里叶转换 (FFT)。其实就是翻译了一下https://blog.csdn.net/weixin\_33714884/article/details/94249560?utm\_medium=distribute.pc\_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-4.channel\_param&depth\_1-utm\_source=distribute.pc\_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-4.channel\_param中的代码。  
实际测试发现两个数组的规模都比较小 (<= 50) 的时候，运行是正常的。但是当数据规模到500时，软件就运行出错。不清楚原因。但是用递归版本就没问题。

下面是代码：

#include <stdio.h> 
    #include <stdlib.h>
    #include <math.h>
    #include <time.h>
    
    double PI = 3.1415926;
    int lim = 1;
    int r[1000];
    int l;
    
    //定义复数 
    struct Complex{
        
    	double x; //实部 
    	double y; //虚部 
    };
    
    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;
    }
    
    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;
    }
    
    //迭代版本的FFT，少了一个参数
    void fft(struct Complex * a, int type)  
    {
        
        for (int i = 0; i < lim; i++)
        {
        
        	if (i < r[i]) //交换 
            {
        
            	struct Complex temp = a[i];
    			a[i] = a[r[i]];
    			a[r[i]] = temp;
    		}
    	}
            
        for (int mid = 1; mid < lim; mid *= 2)
        {
        
            struct Complex wn; //把360度等分成n份，这时的复数。
    		wn.x = cos(PI/mid);
    		wn.y = type * sin(PI/mid);
    		
            for (int R = mid * 2, j = 0; j < lim; j += R)
            {
        
                struct Complex w;
    			w.x = 1;
    			w.y = 0;
                
                for (int k = 0; k < mid; k++) //计算arr中的每一个值 
                {
        
                    struct Complex x = a[j + k];
                    struct Complex y = multiply(w, a[j + k + mid]);
    				
                    a[j + k] = add(x, y); //处理前一半 
                    a[j + k + mid] = myMinus(x, y); //处理后一半 
                    w = multiply(w, wn);
                }
            }
        }
    }
    
    int main()
    {
        
    	int n = 50; //数组arr1的长度为50
    	int m = 50; //数组arr2的长度为50
    	time_t t1, t2;
    	
    	t1 = time(0);
    	
        int complexSize = sizeof(struct Complex);
        
        //都要申请很大的空间，不知道为什么 
        struct Complex * a = (struct Complex *)(malloc(1000000 * complexSize));
    	struct Complex * b = (struct Complex *)(malloc(1000000 * complexSize));
    	struct Complex * c = (struct Complex *)(malloc(1000000 * 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;
    	}
            
    	//lim为1, 2, 4, 8, 16, 32....除了1之外，其余都是偶数
        while (lim < n + m - 1)
        {
        
    		lim = lim * 2;
    		l++;
    	}
        
        //可能是蝴蝶操作 
        for (int i = 0; i < lim; i++)
        {
        
        	//r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));  
        	r[i] = (r[i / 2] / 2) | ((i & 1) << (l - 1)); //位运算符 
    	}
    	
        //即两函数的傅里叶变换的乘积等于它们卷积后的傅里叶变换，
    	//能使傅里叶分析中许多问题的处理得到简化。
    	//F(g(x) * f(x)) = F(g(x)) * F(f(x))，其中F表示的是傅里叶变换 
        fft(a, 1); //快速傅里叶转换 
        fft(b, 1); //快速傅里叶转换
        
        for (int i = 0; i < lim; i++)	
        {
        
        	c[i] = multiply(a[i], b[i]); //把a和b中相乘的结果都放到c中 
    	}
            
        fft(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));
    	}
    	printf("\r\n\r\n");
    	
    	t2 = time(0);
    	printf("迭代的快速傅里叶转换FFT: %d\r\n", t2 - t1);
        return 0;
    }

蝴蝶操作那部分代码是不是必须要用位操作呢？

运行结果如下：  
![在这里插入图片描述][20200913110933713.png_pic_center]

[20200913110933713.png_pic_center]: /images/20221123/c071680e88b0460e93a2cfc40e95d2f6.png