递归版的快速傅里叶转换 (FFT),C语言实现
看了好多个版本的快速傅里叶转换 (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;
}
运行结果如下:
如果main函数中,struct Complex * a、b、c申请的空间比较小 (100或更小),则会出错,不知道为什么。
下次再写迭代版本的FFT吧。
还没有评论,来说两句吧...