FFT-C语言
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-21
一、彻底理解傅里叶变换
快速傅里叶变换(Fast Fourier Transform)是离散傅里叶变换的一种快速算法，简称FFT，通过FFT可以将一个信号从时域变换到频域。
模拟信号经过A/D转换变为数字信号的过程称为采样。为保证采样后信号的频谱形状不失真，采样频率必须大于信号中最高频率成分的2倍，这称之为采样定理。
|
假设采样频率为fs，采样点数为N，那么FFT结果就是一个N点的复数，每一个点就对应着一个频率点，某一点n(n从1开始)表示的频率为：fn=(n-1)*fs/N。 举例说明：用1kHz的采样频率采样128点，则FFT结果的128个数据即对应的频率点分别是0，1k/128，2k/128，3k/128，…，127k/128 Hz。 这个频率点的幅值为：该点复数的模值除以N/2(n=1时是直流分量，其幅值是该点的模值除以N)。 |
二、傅里叶变换的C语言编程
1、码位倒序。
假设一个N点的输入序列，那么它的序号二进制数位数就是t=log2N.
码位倒序要解决两个问题：①将t位二进制数倒序;②将倒序后的两个存储单元进行交换。
如果输入序列的自然顺序号i用二进制数表示，例如若最大序号为15，即用4位就可表示n3n2n1n0，则其倒序后j对应的二进制数就是n0n1n2n3，那么怎样才能实现倒序呢?利用C语言的移位功能!
/*变址计算，将x(n)码位倒置*/ void Reverse(complex *IN_X, int n) { int row = log(n) / log(2); //row为FFT的N个输入对应的最大列数 complex temp; //临时交换变量 unsigned short i = 0, j = 0, k = 0; double t; for (i = 0; i < row; i++) //从第0个序号到第N-1个序号 { k = i; //当前第i个序号 j = 0; //存储倒序后的序号，先初始化为0 t = row; //共移位t次 while ((t--) > 0) //利用按位与以及循环实现码位颠倒 { j = j << 1; j |= (k & 1); //j左移一位然后加上k的最低位 k = k >> 1; //k右移一位，次低位变为最低位 /*j为倒序，k为正序， 从k右向左的值依次移位， j向左依次填入对应的k的右向位*/ } if (j > i) //如果倒序后大于原序数，就将两个存储单元进行交换(判断j>i是为了防止重复交换 { temp = IN_X[i]; IN_X[i] = IN_X[j]; IN_X[j] = temp; } } } //复数加法的定义 complex add(complex a, complex b) { complex c; c.real = a.real + b.real; c.img = a.img + b.img; return c; } //复数乘法的定义 complex mul(complex a, complex b) { complex c; c.real = a.real*b.real - a.img*b.img; c.img = a.real*b.img + a.img*b.real; return c; } //复数减法的定义 complex sub(complex a, complex b) { complex c; c.real = a.real - b.real; c.img = a.img - b.img; return c; }
2、第二个要解决的问题就是蝶形运算
对N个输入行，m表示第m列蝶形。
①第1级(第1列)每个蝶形的两节点“距离”为1，第2级每个蝶形的两节点“距离”为2，第3级每个蝶形的两节点“距离”为4，第4级每个蝶形的两节点“距离”为8。由此推得，
第m级蝶形运算，每个蝶形的两节点“距离”L=2^(m-1)。
②对于16点的FFT，第1级有16组蝶形，每组有1个蝶形;第2级有4组蝶形，每组有2个蝶形;第3级有2组蝶形，每组有4个蝶形;第4级有1组蝶形，每组有8个蝶形。由此可推出，
对于N点的FFT，第m级有N/(2L)组蝶形，每组有L=2^(m-1)个蝶形。
③旋转因子
的确定
以16点FFT为例，第m级第k个旋转因子为
，其中k=0～2^(m-1)-1，即第m级共有2^(m-1)个旋转因子，根据旋转因子的可约性，
，所以第m级第k个旋转因子为
，其中k=0～2^(m-1)-1。
生成一个的序列，代码如下：
/*产生一个周期欧拉形式的Wn的值*/ void InitWn(complex *w,int n) //n为输入的个数,w为权值Wn { int k; for (k = 0; k < n; k++) { w[k].real = cos(2 * PI / n * k); //用欧拉公式计算旋转因子 w[k].img = -1 * sin(2 * PI / n * k); } }
3、算法实现
我们已经知道，N点FFT从左到右共有log2N级蝶形，每级有N/2L组，每组有L个。所以FFT的C语言编程只需用3层循环即可实现：最外层循环完成每一级的蝶形运算(整个FFT共log2N级)，中间层循环完成每一组的蝶形运算(每一级有N/2L组)，最内层循环完成单独1个蝶形运算(每一组有L个)。
/*快速傅里叶变换*/ void fft(complex *IN_X, int n) { int row = log(n) / log(2); //row为FFT的N个输入对应的最大列数 int i = 0, j = 0, k = 0, L = 0; complex top, bottom, product; Reverse(IN_X,n); //调用变址函数 for (i = 0; i < row ; i++) /*第i级蝶形运算 */ { L = 1 << i; //L等于2的i次方,表示每一级蝶形中组间的距离，即一组中蝶形个数 for (j = 0; j < n; j += 2 * L) /*j为组偏移地址，每L个蝶形是一组，每级有N/2L组*/ { for (k = 0; k < L; k++) /*k为一组中蝶形的偏移地址，一个蝶形运算 每个group内的蝶形运算*/ { product = mul(IN_X[j + k + L], Wn[n*k/2/L]); top = add(IN_X[j + k], product); bottom = sub(IN_X[j + k], product); IN_X[j + k] = top; IN_X[j + k + L] = bottom; } } } }
4.全部代码
#include <stdio.h> #include <math.h> #include <stdlib.h> #include <malloc.h> #define N 1000 #define PI atan(1)* 4 /*定义复数类型*/ typedef struct { double real; double img; }complex; complex x[N], *Wn; /*输入序列,复数系数Wn*/ int size_x; /*size_x为采样的信号个数，必须为2的倍数*/ void fft(complex *IN_X, int n); /*对输入IN_X，快速傅里叶变换*/ void InitWn(complex *w, int n); /*生成一个长度为n的Wn欧拉形式序列*/ void Reverse(complex *IN_X, int n); /*对输入IN_X地址*/ complex add(complex, complex); /*复数加法*/ complex mul(complex, complex); /*复数乘法*/ complex sub(complex, complex); /*复数减法*/ void output(); /*输出快速傅里叶变换的结果*/ int main() { int i; /*输出结果*/ system("cls"); //调用系统命令，CLS的功能为清除屏幕输出 printf("输出DIT方法实现的FFT结果\n"); printf("Please input the size of x:\n");//输入序列的大小 scanf_s("%d", &size_x); printf("Please input the data in x[N]:\n");//输入序列的实部和虚部 for (i = 0; i < size_x; i++) { printf("请输入第%d个序列：", i); scanf_s("%lf%lf", &x[i].real, &x[i].img); } printf("输出倒序后的序列\n"); Wn = (complex *)malloc(sizeof(complex) * size_x); //分配变换后的值的内存空间 InitWn(Wn , size_x);//调用变换核 fft(x, size_x);//调用快速傅里叶变换 printf("输出FFT后的结果\n"); output();//调用输出傅里叶变换结果函数 scanf_s("%d", &size_x); //free(Wn); return 0; } /*快速傅里叶变换*/ void fft(complex *IN_X, int n) { int row = log(n) / log(2); //row为FFT的N个输入对应的最大列数 int i = 0, j = 0, k = 0, L = 0; complex top, bottom, product; Reverse(IN_X,n); //调用变址函数 for (i = 0; i < row ; i++) /*第i级蝶形运算 */ { L = 1 << i; //L等于2的i次方,表示每一级蝶形中组间的距离，即一组中蝶形个数 for (j = 0; j < n; j += 2 * L) /*j为组偏移地址，每L个蝶形是一组，每级有N/2L组*/ { for (k = 0; k < L; k++) /*k为一组中蝶形的偏移地址，一个蝶形运算 每个group内的蝶形运算*/ { product = mul(IN_X[j + k + L], Wn[n*k/2/L]); top = add(IN_X[j + k], product); bottom = sub(IN_X[j + k], product); IN_X[j + k] = top; IN_X[j + k + L] = bottom; } } } } /*产生一个周期欧拉形式的Wn的值*/ void InitWn(complex *w,int n) //n为输入的个数,w为权值Wn { int k; for (k = 0; k < n; k++) { w[k].real = cos(2 * PI / n * k); //用欧拉公式计算旋转因子 w[k].img = -1 * sin(2 * PI / n * k); } } /*变址计算，将x(n)码位倒置*/ /* 码位倒序要解决两个问题： ①将t位二进制数倒序;②将倒序后的两个存储单元进行交换。 如果输入序列的自然顺序号i用二进制数表示， 例如若最大序号为15，即用4位就可表示n3n2n1n0， 则其倒序后j对应的二进制数就是n0n1n2n3 */ void Reverse(complex *IN_X, int n) { int row = log(n) / log(2); //row为FFT的N个输入对应的最大列数 complex temp; //临时交换变量 unsigned short i = 0, j = 0, k = 0; double t; for (i = 0; i < row; i++) //从第0个序号到第N-1个序号 { k = i; //当前第i个序号 j = 0; //存储倒序后的序号，先初始化为0 t = row; //共移位t次 while ((t--) > 0) //利用按位与以及循环实现码位颠倒 { j = j << 1; j |= (k & 1); //j左移一位然后加上k的最低位 k = k >> 1; //k右移一位，次低位变为最低位 /*j为倒序，k为正序， 从k右向左的值依次移位， j向左依次填入对应的k的右向位*/ } if (j > i) //如果倒序后大于原序数，就将两个存储单元进行交换(判断j>i是为了防止重复交换 { temp = IN_X[i]; IN_X[i] = IN_X[j]; IN_X[j] = temp; } } } /*输出傅里叶变换的结果*/ void output() { int i; printf("The result are as follows：\n"); for (i = 0; i < size_x; i++) { printf("%.4f", x[i].real); if (x[i].img >= 0.0001)printf("+%.4fj\n", x[i].img); else if (fabs(x[i].img) < 0.0001)printf("\n"); else printf("%.4fj\n", x[i].img); } } //复数加法的定义 complex add(complex a, complex b) { complex c; c.real = a.real + b.real; c.img = a.img + b.img; return c; } //复数乘法的定义 complex mul(complex a, complex b) { complex c; c.real = a.real*b.real - a.img*b.img; c.img = a.real*b.img + a.img*b.real; return c; } //复数减法的定义 complex sub(complex a, complex b) { complex c; c.real = a.real - b.real; c.img = a.img - b.img; return c; }
