|
FFT 算法原理 首先，为了能够进行FFT，需要了解DFT。 因为两者之间有本质上相同的东西。 在此之前，先列举离散傅立叶变换对(DFT )。 但是FFT之所以被称为快速傅立叶变换，是利用了以下性质() 先把这三个公式放在这里，接下来用。 这些特性使得可以将长DFT计算分解为几个短序列的DFT计算的组合，从而减少计算量。 为了便于理解，我们使用了按时间提取的快速傅立叶变换(DIT-FFT )方法。 首先，将一个序列x(n )分成两部分。 在以下情况下： 设n为2的整数次方，即N=2^M。 按奇偶校验对x(n )进行分组： DFT也分为两组进行预算。 此时，在我们上面列举的三个性质中，约定性很有用。 让我们回顾一下： 该公式如下所示。 x1(k )和x2 ) k )都是长度为N/2的序列x1 ) k )和x2 ) k )的N/2点的离散傅立叶变换。 其中： 这样，一个n点的DFT被分解为两个N/2的DFT。 但是，x1(k )和x1(k )只有N/2分。 也就是说，公式(1-1)只是DFT的前半部分。 DFT的后半部分被要求能够利用其对称性求出后半部分。 式(1-1)和式(1-2)可以用蝶形信号流程图符号表示。 图： 以N=8为例，可以用下图表示。 那么，通过这样分解，每1个N/2点DFT，只需要(N^2)/4次的复数乘法，就可以明显地节约运算量。 那么，接下来类推，是继续按奇偶校验分解已经出现的x1(k )和x1(k )，还是与上边相同的方法？ 在百度上也能找到的这张照片出山了。 对于这张图像，我想强调的是二次蝶形运算的时候。 之后不应该变成那样吗？ 为什么会变成这样？ 这个问题之前我烦恼了一会儿，但别忘了。 前者的展开是N/2，因为此时n是奇偶校验分离的，是从根据的可约束性得出的，在此不能混淆。 没有必要提到计算效率。 m级计算总共需要： FFT c语言实现 FFT算法c语言的实现，网络方法层出不穷，除了本人懒惰(不擅长看别人的程序)外，还有自给自足的衣食住原则，所以我自己也写了个人容易理解的程序。 而且，为了帮助读者理解，我特意尽量减少了库函数的使用。 一些基本函数是自己写的。 但是，作为FFT算法已经足够了。 目前，该程序只能处理2^N的数据。 理论上，如果2^N不够，可以进行在后面的数列中添加0的操作。 当然，因为简约语句也不能实现，但主要功能如下开源。 左右滑动查看代码 /* 功能：对input内的数据进行快速傅立叶变换并输出 */ #包含 #包含 #define FFT_LENGTH 8 double input [ FFT _ length ]={ 1，1，1，1，1，1，1，1 }； 定义结构完成1 {//多个结构 双现实； //实部 双图像； //虚部 (； 将input的实数结果保存在复数中 结构完成1 result _ dat [8]； /* 虚数乘法 */ 结构完成1 con _ complex (结构完成1 a，结构完成1 b ) { struct complex1 temp;
temp.real=(a.real*b.real)-(a.image*b.image); temp.image=(a.image*b.real)+(a.real*b.image); return temp; } /* 简单的a的b次方 */ int mypow(int a,int b){ int i,sum=a; if(b==0)return 1; for(i=1;i sum*=a; } return sum; } /* 简单的求以2为底的正整数 */ int log2(int n){ unsigned i=1; int sum=1; for(i;;i++){ sum*=2; if(sum>=n)break; } return i; } /* 简单的交换数据的函数 */ void swap(struct complex1 *a,struct complex1 *b){ struct complex1 temp; temp=*a; *a=*b; *b=temp; } /* dat为输入数据的数组 N为抽样次数 也代表周期 必须是2^N次方 */ void fft(struct complex1 dat[],unsigned char N){ /*最终 dat_buf计算出 当前蝶形运算奇数项与W 乘积 dat_org存放上一个偶数项的值 */ struct complex1 dat_buf,dat_org; /* L为几级蝶形运算 也代表了2进制的位数 n为当前级蝶形的需要次数 n最初为N/2 每级蝶形运算后都要/2 i j为倒位时要用到的自增符号 同时 i也用到了L碟级数 j是计算当前碟级的计算次数 re_i i_copy均是倒位时用到的变量 k为当前碟级 cos(2*pi/N*k)的 k 也是e^(-j2*pi/N)*k 的 k */ unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1; //经过观察，发现每级蝶形运算需要N/2次运算，共运算N/2*log2N 次 unsigned char fft_counter=0; //在此要进行补2 N必须是2^n 在此略 //蝶形级数 (L级) L=log2(N); //计算每级蝶形计算的次数(这里只是一个初始值) 之后每次要/2 //n=N/2; //对dat的顺序进行倒位 for(i=1;i i_copy=i; re_i=0; for(j=L-1;j>0;j--){ //判断i的副本最低位的数字 并且移动到最高位 次高位 .. //re_i为交换的数 每次它的数字是不能移动的 并且循环之后要清0 re_i|=((i_copy&0x01)< i_copy>>=1; } swap(&dat,&dat[re_i]); } //进行fft计算 for(i=0;i fft_flag=1; fft_counter=0; for(j=0;j if(fft_counter==mypow(2,i)){ //控制隔几次，运算几次 fft_flag=0; }else if(fft_counter==0){ //休止结束，继续运算 fft_flag=1; } //当不判断这个语句的时候 fft_flag保持 这样就可以持续运算了 if(fft_flag){ dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1))); dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1))); dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf); //计算 当前蝶形运算奇数项与W 乘积 dat_org.real=dat[j].real; dat_org.image=dat[j].image; //暂存 dat[j].real=dat_org.real+dat_buf.real; dat[j].image=dat_org.image+dat_buf.image; //实部加实部 虚部加虚部 dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real; dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image; //实部减实部 虚部减虚部 k++; fft_counter++; }else{ fft_counter--; //运算几次，就休止几次 k=0; } } } } void main{ int i; //先将输入信号转换成复数 for(i=0;i result_dat.image=0; //输入信号是二维的，暂时不存在复数 result_dat.real=input; //result_dat.real=10; //输入信号都为实数 } fft(result_dat,FFT_LENGTH); for(i=0;i input=sqrt(result_dat.real*result_dat.real+result_dat.image*result_dat.image); //取模 printf("%lfn",input); } while(1); } 这个程序中input这个数组是输入信号，在这里只模拟抽样了8次，输出的数据也是input，如果想看其它序列的话，可以改变FFT_LENGTH 的值以及 input里的内容，程序输出的是实部和虚部的模
| 
楼主
标题置顶
标题高亮
