《2022年数字信号处理实验FFT快速傅里叶变换C语言》由会员分享,可在线阅读,更多相关《2022年数字信号处理实验FFT快速傅里叶变换C语言(11页珍藏版)》请在金锄头文库上搜索。
1、个人资料整理仅限学习使用1 / 11 数字信号处理实验 是按正常顺序排列在存储单元中,即按 X(0,X(1, ,X(7的顺序排列,但是这时输入x(n却不是按自然顺序存储的,而是按x(0,x(4 , ,x(7 的顺序存入存储单元,所以我们要对输入的按正常顺序排列的数据进行变址存储,最后才能得到输出的有序的X(K 。通过观察,可以发现,如果说输出数据是按原位序排列的话,那么输入数据是按倒位序排列的。即如果输入序列的序列号用二进制数,则到位序就为。我们需将输入的数据变为输出的倒位序存储,这里用雷德算法来实现。下面给出雷德算法。DXDiTa9E3d 假如使用 AI存的是顺序位序,而BJ存的是倒位序。例
2、如 N = 8 的时候,倒位序顺序二进制表示倒位序顺序0 0 000 000RTCrpUDGiT 4 1 100 0015PCzVD7HxA 2 2 010 010 jLBHrnAILg 6 3 110 011xHAQX74J0X 1 4 001 100LDAYtRyKfE 5 5 101 101Zzz6ZB2Ltk 3 6 011 110dvzfvkwMI1 7 7 111 111rqyn14ZNXI 由上面的表可以看出,按自然顺序排列的二进制数,其下面一个数总是比其上面一个数大 1,即下面一个数是上面一个数在最低位加1 并向高位进位而得到的。而倒位序二进制数的下面一个数是上面一个数在最高位
3、加1 并由高位向低位进位而得到。 I、J 都是从 0 开始,若已知某个倒位序J,要求下一个倒位序数,则应先判断J 的最高位是否为0,这可与 k=N/2 相比较,因为 N/2 总是等于 100.的。如果 kJ,则 J 的最高位为 0,只要把该位变为1J 与 k=N/2相加即可),就得到下一个倒位序数;如果K=J,则 J 的最高位为 1,可将最高位变为 0J 与 k=N/2 相减即可)。然后还需判断次高位,这可与k=N4 相比较,若次高位为 0,则需将它变为 1加 N4 即可)其他位不变,既得到下一个倒位序数;若次高位是1,则需将它也变为0。然后再判断下一位EmxvxOtOco 。2. 复数运算由
4、于每一个蝶形结构完成的迭代运算为算式中涉及到了复数的运算,而计算机是不能自己实现复数运算的,所以需要我们自己设计进行复数运算的程序。迭代运算式中,= cos2 r/N )- jsin + jI(j, = R(K + jI(K ,SixE2yXPq5 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页,共 11 页个人资料整理仅限学习使用3 / 11 而我们最后希望得到的DFT 结果是复数的模,根据它的模来绘制频谱,所以这里我们定义X(k= 相关程序我们编译为:c.real=a.real*b.real-a.imag*b.imag。 c.imag=
5、a.real*b.imag+a.imag*b.real。根据迭代运算的式子,我们可以将其分解为:R(K+jI(K=R(K+jI(K+ R(j + jI(j* cos2r/N)-jsin= R(K+ R(j cos sin (2r/N I(K=I(K-R(j sincos (2r/N 6ewMyirQFL 同理R(j= R(K- R(j cos sin (2r/N I(j=I(K+R(j sincos (2r/N kavU42VRUs 相关程序编译为:xinip.real=xini.real-t.real。 xinip.imag=xini.imag-t.imag。 xini.real=xini.
6、real+t.real。 xini.imag=xini.imag+t.imag。3. 节点距离运算当输入为倒位序,输出为正常顺序时,第m 级运算,每个蝶形的两节点距离为。4. 旋转因子的计算这里主要解决的是迭代运算的每一项应乘的旋转因为中r 应取多少的问题。第 L 级的 2(L-1 个碟形因子 WPN 中的 P,可表示为 p = j*2(m-L,其中j = 0,1,2,. , 个蝶形因子,第二层循环根据乘数进行控制,保证对于每一个旋转因子第三层循环要执行一次,这样,第三层循环在第二层循环控制下,每一级要进行2(L-1 次循环计算。第三层:由于第L 级共有 N/2L 个蝶形结构,并且同一级内不同
7、蝶型结构的旋转因子分布相同,当第二层循环确定某一旋转因子后,第三层循环要将本级中每个蝶型结构中具有这一旋转引自的蝶形计算一次,即第三层循环每执行完一次要进行N/2L 个蝶形计算。所以,在每一级中,第三层循环完成 N/2L 个蝶形计算;第二层循环使得第三层循环进行 2(L-1 次,因此,第二层循环完成时,共进行2(L-1 *N/2L=N/2个碟形计算。实质是:第二、第三层循环完成了第L 级的计算。0YujCfmUCw 五、程序代码精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页,共 11 页个人资料整理仅限学习使用4 / 11 1.C 语言:#
8、include #include #include #define PI 3.1415926535897932384626433832795028841971 eUts8ZQVRd #define FFT_N 128 /定义傅利叶变换的点数sQsAEJkW5T struct compx double real,imag。 。 /定义一个复数结构GMsIasNXkA struct compx sFFT_N。 /从 S1 开始存放TIrRGchYzg struct compx EE(struct compx a,struct compx b /求复数的模长函数7EqZcWLZNX struct c
9、ompx c。 c.real=a.real*b.real-a.imag*b.imag。 c.imag=a.real*b.imag+a.imag*b.real。 return(c。 void FFT(struct compx *xin /FFT函数 int f,m,nv2,nm1,i,k,l,j=0。 struct compx u,w,t。 nv2=FFT_N/2 。 /变址运算,即把自然顺序变成倒位序,采用雷德算法lzq7IGf02E nm1=FFT_N-1。 for(i=0。i if(i /如果 ij,即进行变址 t=xinj。 xinj=xini。 xini=t。 k=nv2。 /求 j
10、的下一个倒位序 while(k /如果 k!=1。l+ /第一层循环,计算蝶形级数zvpgeqJ1hk 。 for(m=1。m /第二层循环,控制蝶形级数NrpoJac3v1 d=2。 /d蝶形结距离,即第m级蝶形的蝶形结构相距 d 点1nowfTG4KI d1=d/2。 /同一蝶形结中参加运算的两点的距离fjnFLDa5Zo u.real=1.0。 /u为蝶形结构运算系数,初始值为 1tfnNhnE6e5 u.imag=0.0。 w.real=cos(PI/d1。 /w为系数商,即当前系数与 前 一 个 系 数 的 商HbmVN777sL w.imag=-sin(PI/d1。 for(j=0
11、。j /第三层循环,进行旋转因子的计算完成蝶形计算。这里控制 计 算 系 数 不 同 的蝶 形 结V7l4jRB8Hs for(i=j。i /控制计算系数相同蝶形结 ip=i+d1。 /i,ip分别表示参加蝶形运算的两个节点83lcPA59W9 t=EE(xinip,u。 /蝶形运算 xinip.real=xini.real-t.real。 xinip.imag=xini.imag-t.imag。 xini.real=xini.real+t.real。 xini.imag=xini.imag+t.imag。 u=EE(u,w。 /改变系数,进行下一个蝶形运算mZkklkzaaP void ma
12、in( int i。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页,共 11 页个人资料整理仅限学习使用6 / 11 for(i=0。i /给结构体赋值AVktR43bpw si.real=sin(2*3.141592653589793*i/FFT_N。 si.imag=0。 long start=clock(。 FFT(s 。 for(i=0。i si.real=sqrt(si.real*si.real+si.imag*si.imag。ORjBnOwcEd long end=clock(。 printf(%.4fn,si.real。 lo
13、ng t=end-start。printf(%dn,t。 2.MATLAB :n=0:127。b=sin(2*pi*n/128。fft(b 六、实验结果1.C 语言精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页,共 11 页个人资料整理仅限学习使用7 / 11 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页,共 11 页个人资料整理仅限学习使用8 / 11 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页,共 11 页个人资料整理仅限学
14、习使用9 / 11 2.MATLAB Columns 1 through 5 0.0000 -0.0000 -64.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i2MiJTy0dTT Columns 6 through 10 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000igIiSpiue7A Columns 11 through 15 0.0000 + 0.0000i -0.0000 - 0.0000
15、i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000iuEh0U1Yfmh Columns 16 through 20 -0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000iIAg9qLsgBX Columns 21 through 25 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000iWwghWvVhPE
16、Columns 26 through 30 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000iasfpsfpi4k Columns 31 through 35 -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000iooeyYZTjj1 Columns 36 through 40 -0.0000 - 0.0000i 0.0000 + 0.0000i 0
17、.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000iBkeGuInkxI Columns 41 through 45 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页,共 11 页个人资料整理仅限学习使用10 / 11 -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000iPgdO0sRlMo Columns 46 through 50 -0.0000 + 0.0000i -0
18、.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i3cdXwckm15 Columns 51 through 55 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000ih8c52WOngM Columns 56 through 60 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0
19、000iv4bdyGious Columns 61 through 65 -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 0.0000 J0bm4qMpJ9 Columns 66 through 70 -0.0000 -0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000iXVauA9grYP Columns 71 through 75 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i
20、0.0000 + 0.0000i 0.0000 + 0.0000ibR9C6TJscw Columns 76 through 80 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000ipN9LBDdtrd Columns 81 through 85 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000iDJ8T7nHuGT Columns 86 throug
21、h 90 -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000iQF81D7bvUA Columns 91 through 95 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i4B7a9QFw9h Columns 96 through 100 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i -0
22、.0000 + 0.0000i -0.0000 - 0.0000iix6iFA8xoX Columns 101 through 105 -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000iwt6qbkCyDE Columns 106 through 110 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000iKp5zH46zRk Columns 111 th
23、rough 115 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000iYl4HdOAA61 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页,共 11 页个人资料整理仅限学习使用11 / 11 Columns 116 through 120 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000ich4P
24、Jx4BlI Columns 121 through 125 -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000iqd3YfhxCzo Columns 126 through 128 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 +64.0000iE836L11DO5 Elapsed time is 0.000120 seconds. 七、实验分析根据对比可知,我自己设计的C 程序中, FFT 计算时间为70 毫秒,而MATLAB 中,
25、FFT 计算时间为 0.120 毫秒。我认为计算时间其实跟算法设计是有很大关系的,我自己的设计程序算法可能还是不是简便,MATLAB 自带的 FFT函数程序中应该是设计更优化,算法时间更简便的。但两个计算结果都是相同的,并能得到很好的与原信号频谱的拟合效果,说明它们采用的都是FFT 算法。S42ehLvE3M 八、实验心得此次实验由于C 语言基础并不是特别好,而且自从大一后就没有再使用过,所以前期构思时完全属于摸不着头脑,想了好几天。脑袋里有大概的思路,但是并不能形成系统完整的程序。于是到网上搜了很多相关资料和程序进行学习和参考,将这些程序全都弄懂并学习了相关C 语言语法后,才开始进行编程。但
26、编程时仍然存在心里有想法却无法转换成C 语言的问题。认为可能是C 语言不够熟练,对蝶形运算的理解还不够透彻的原因。网上大多数程序的思路是相同的,我想这是由于FFT 这个算法本身已经是固定算法的原因。我从中选出一个我理解的非常透彻的程序进行参考,编出了相关的程序。可能方法比较笨,但已经是我努力过后的结果了。我认为虽然此次大作业比较费心血,而且能力不是特别好,导致程序没有调好之前的几个晚上睡前脑中还是这些程序要怎么调,但是通过这次大作业现在我已经对蝶形运算有了很深的理解,并掌握了整个 FFT 算法的思路,我觉得这就是收获。我认为其实作业完成好坏无关于有没有借鉴,重点是这其中有了收获,并使相关知识水平得到了提升,这就是作业完成了。大作业虽然真的让我很苦恼,而且每次都要憋好久,但是我认为我在这其中学到了很多,这些收获远比平时听课要收获的多,所以我觉得老师的大作业是非常好的促进我们学习的方法。501nNvZFis 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页,共 11 页