2022年数字信号处理实验FFT快速傅里叶变换C语言.docx
精选学习资料 - - - - - - - - - 个人资料整理 仅限学习使用数字信号处理试验一、试验内容<第八章 FFT）针对典型序列,用 的运算结果进行比较；二、试验工具1.VC+6.0 2.matlab2022 三、试验涉及学问图 1 C 语言编程实现基 2-FFT,并与 MATLAB自带的 FFT函数按时间抽选的基 2FFT算法 依据蝶形运算图,如图 1,可以看出,一个点基 2-FFT运算由 n 级蝶形运算构成,每一个蝶形结构完成下属基本迭代运算：b5E2RGbCAP 式中 m 表示第 m 列迭代, k,j 为数据所在行数；上式的蝶形运算如图 2 所示；图 2 1 / 11 名师归纳总结 - - - - - - -第 1 页,共 11 页精选学习资料 - - - - - - - - - 个人资料整理 仅限学习使用四、试验设计思路第一,依据基 2-FFT 的算法特点,可以整理出程序设计的大致思路；基 2-FFT主要运用的就是蝶形运算来降低进行 DFT运算时的运算量；由于是基 2,所以进行 DFT运算的点数 N 必需的；程序设计中主要解决 3 个问题,变址运算,复数运算,节点距离运算,及旋转因子的运算；下面对这三个模块进行说明；p1EanqFDPw 1. 变址运算 由蝶形图我们可以看出,FFT的输出 Xk>是按正常次序排列在储备单元中,即按 X0>,X1>, ,X 7>的次序排列,但是这时输入 xn>却不是按自然次序储备的,而是按 x0>,x4>, ,x7>的次序存入储备单元,所以我们要对输入的按正常次序排列的数据进行变址储备,最终才能得到输出的有序的 XK>；通过观看,可以发觉,假如说输出数据是按原位序排列的话,那么输入数据是按倒位序排列的；即假如输入序列的序列号用二进制数,就到位序就为；我们需将输入的数据变为输出的倒位序储备,这里用雷德算法来实现；下面给出雷德算法；DXDiTa9E3d 假如使用 AI存的是次序位序,而BJ存的是倒位序；例如 N = 8 的时候,倒位序次序 二进制表示 倒位序次序 0 0 000 000 RTCrpUDGiT 4 1 100 001 5PCzVD7HxA 2 2 010 010 jLBHrnAILg 6 3 110 011 xHAQX74J0X 1 4 001 100 LDAYtRyKfE 5 5 101 101 Zzz6ZB2Ltk 3 6 011 110 dvzfvkwMI1 7 7 111 111 rqyn14ZNXI 由上面的表可以看出,按自然次序排列的二进制数,其下面一个数总是比 其上面一个数大 1,即下面一个数是上面一个数在最低位加 1 并向高位进位而 得到的；而倒位序二进制数的下面一个数是上面一个数在最高位加 1 并由高位 向低位进位而得到； I、J 都是从 0 开头,如已知某个倒位序 J,要求下一个倒 位序数,就应先判定 J 的最高位是否为 0,这可与 k=N/2 相比较,由于 N/2 总 是等于 100.的；假如 k>J,就 J 的最高位为 0,只要把该位变为 1<J 与 k=N/2相加即可）,就得到下一个倒位序数；假如K<=J ,就 J 的最高位为 1,可将最高位变为 0<J 与 k=N/2 相减即可）；然后仍需判定次高位,这可与 k=N4 相比 较,如次高位为 0,就需将它变为 1<加 N4 即可）其他位不变,既得到下一个倒位序数；如次高位是1,就需将它也变为0；然后再判定下一位EmxvxOtOco ；2. 复数运算 由于每一个蝶形结构完成的迭代运算为算式中涉及到了复数的运算,而运算机是不能自己实现复数运算的,所以需要 我们自己设计进行复数运算的程序；迭代运算式中,= cos<2 r/N ）- jsin<2 r/N ）我们设 = Rj> + jIj>, = RK> + jIK> ,SixE2yXPq5 2 / 11 名师归纳总结 - - - - - - -第 2 页,共 11 页精选学习资料 - - - - - - - - - 个人资料整理 仅限学习使用而我们最终期望得到的 里我们定义Xk>= 相关程序我们编译为：DFT 结果是复数的模,依据它的模来绘制频谱,所以这c.real=a.real*b.real-a.imag*b.imag； c.imag=a.real*b.imag+a.imag*b.real；依据迭代运算的式子,我们可以将其分解为：RK>+jIK>=RK>+jIK>+ Rj> + jIj>* cos<2 r/N）-jsin<2 r/N）；连续分解得到以下两式 : RK>= RK>+ Rj> cos<2 r/N ）+ Ij> sin 2 r/N> IK>=IK>-Rj> sin<2 r/N ）+Ij>cos 2 r/N> 6ewMyirQFL 同理Rj>= RK>- Rj> cos<2 r/N ）- Ij> sin 2 r/N> Ij>=IK>+Rj> sin<2 r/N ）Ij>cos 2 r/N> kavU42VRUs 相关程序编译为：xinip.real=xini.real-t.real； xinip.imag=xini.imag-t.imag； xini.real=xini.real+t.real； xini.imag=xini.imag+t.imag；3. 节点距离运算当输入为倒位序,输出为正常次序时,第 离为；4. 旋转因子的运算m 级运算,每个蝶形的两节点距这里主要解决的是迭代运算的每一项应乘的旋转由于中 r 应取多少的问题；第 L 级的 2L-1> 个碟形因子 WPN 中的 P,可表示为 p = j*2m-L>,其中j = 0,1,2,. ,<2L-1-1 ）； y6v3ALoS89 将上述方法整理一下,我们可以得到一个完整的程序的思路；主要分为三个模块,复数运算 DFT 模函数, FFT 函数,主函数；其中复数运算 DFT模函数的算法我们已经在前面介绍了；主函数是编译整个程序进行的次序,即我们先设定一个时域函数,从中抽取 N 点的值,然后进行 FFT 运算,然后求模长,最后将对应的频域的值输出；这里比较复杂的是 M2ub6vSTnP FFT 函数；下面进行说明；这里我们共设三个循环来实现 FFT这一函数的功能；第一层：第一我们知道,一个N=2m点 DFT,要进行 m级运算 , 所以第一层循环是对运算的级数进行掌握；其次层：由于第L 级有 2L-1> 个蝶形因子,其次层循环依据乘数进行掌握,保证对于每一个旋转因子第三层循环要执行一次,这样,第三层循环在其次层循环掌握下,每一级要进行 2L-1> 次循环运算；第三层：由于第 L 级共有 N/2L 个蝶形结构,并且同一级内不同蝶型结构的旋转因子分布相同,当其次层循环确定某一旋转因子后,第三层循环要将本级中每个蝶型结构中具有这一旋转引自的蝶形运算一次,即第三层循环每执行完一次要进行 N/2L 个蝶形运算；所以,在每一级中,第三层循环完成 N/2L 个蝶形运算；其次层循环使得第三层循环进行 2L-1> 次,因此,其次层循环完成时,共进行 2L-1> *N/2L=N/2 个碟形运算；实质是：其次、第三层循环完成了第 L 级的运算；0YujCfmUCw 五、程序代码3 / 11 名师归纳总结 - - - - - - -第 3 页,共 11 页精选学习资料 - - - - - - - - - 个人资料整理 仅限学习使用1.C 语言：#include<stdio.h> #include<math.h> #include<time.h> #define PI 3.1415926535897932384626433832795028841971 eUts8ZQVRd #define FFT_N 128 / 点数 sQsAEJkW5T 定义傅利叶变换的struct compx double real,imag； ； /定义一个 从 S1 开头复数结构GMsIasNXkA struct compx sFFT_N； /存放 TIrRGchYzg struct compx EEstruct compx a,struct compx b> /数7EqZcWLZNX struct compx c； c.real=a.real*b.real-a.imag*b.imag c.imag=a.real*b.imag+a.imag*b.real returnc>； void FFTstruct 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 ； fori=0；i<nm1；i+> 假如 i<j,即进行变址 ifi<j> / t=xinj； xinj=xini； xini=t k=nv2； /求 j 的下一个倒位序 whilek<=j> /假如 k<=j, 表示 j 的最高位为 1 ； /把最高位变成 0 j=j-k k=k/2； /比较次高位,依次类推,逐个比较,直到某个位为 0 j=j+k； /把 0 改为 1 4 / 11 名师归纳总结 - - - - - - -第 4 页,共 11 页精选学习资料 - - - - - - - - - ；个人资料整理仅限学习使用 int d,d1,ip f=FFT_N；l+> /第一层循环,运算蝶形 forl=1；f=f/2>.=1级数 zvpgeqJ1hk ； form=1；m<=l； m+> /其次层循环,掌握蝶形级数 NrpoJac3v1 d=2<<m-1>； /d蝶形结距离,即第m 级蝶形的蝶形结 d1=d/2； /构相距 d 点1nowfTG4KI 同一蝶形结中参与运算的两为蝶形结构运算系数,点的距离fjnFLDa5Zo u.real=1.0； /u初始值为 1tfnNhnE6e5 u.imag=0.0； /w为系数商,即当前系数 w.real=cosPI/d1>与 前 一 个 系 数 的 商HbmVN777sL w.imag=-sinPI/d1>； forj=0；j<=d1-1 ；j+> / 第三层循环,进行旋转因子的运算完成蝶形运算；这里控 制 计 算 系 数 不 同 的 蝶 形 结V7l4jRB8Hs fori=j；i<=FFT_N-1；i=i+d> /掌握运算系数相同蝶形结 ip=i+d1； /i；,ip分别表示参与蝶形运算的两个节点83lcPA59W9 蝶形运算 t=EExinip,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=EEu,w>； /转变系数,进行下一个蝶形运算mZkklkzaaP void main> int i；5 / 11 名师归纳总结 - - - - - - -第 5 页,共 11 页精选学习资料 - - - - - - - - - fori=0； i<FFT_N； i+> /个人资料整理仅限学习使用给结构体赋值AVktR43bpw si.real=sin2*3.141592653589793*i/FFT_N>； si.imag=0； long start=clock>； FFTs> ； fori=0；i<FFT_N；i+> ；ORjBnOwcEd si.real=sqrtsi.real*si.real+si.imag*si.imag>long end=clock>； printf"%.4fn",si.real>； long t=end-start；printf"%dn",t>； 2.MATLAB：n=0:127；b=sin2*pi*n/128> fftb> ；六、试验结果1.C 语言6 / 11 名师归纳总结 - - - - - - -第 6 页,共 11 页精选学习资料 - - - - - - - - - 个人资料整理 仅限学习使用7 / 11 名师归纳总结 - - - - - - -第 7 页,共 11 页精选学习资料 - - - - - - - - - 个人资料整理 仅限学习使用8 / 11 名师归纳总结 - - - - - - -第 8 页,共 11 页精选学习资料 - - - - - - - - - 个人资料整理 仅限学习使用2.MATLAB Columns 1 through 5 0.0000 -0.0000 -64.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i 2MiJTy0dTT Columns 6 through 10 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i gIiSpiue7A Columns 11 through 15 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i uEh0U1Yfmh Columns 16 through 20 -0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i IAg9qLsgBX Columns 21 through 25 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i WwghWvVhPE Columns 26 through 30 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i asfpsfpi4k Columns 31 through 35 -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i ooeyYZTjj1 Columns 36 through 40 -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i BkeGuInkxI Columns 41 through 45 9 / 11 名师归纳总结 - - - - - - -第 9 页,共 11 页精选学习资料 - - - - - - - - - 个人资料整理 仅限学习使用 -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i PgdO0sRlMo Columns 46 through 50 -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i 3cdXwckm15 Columns 51 through 55 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i h8c52WOngM Columns 56 through 60 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i v4bdyGious 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.0000i XVauA9grYP Columns 71 through 75 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i bR9C6TJscw Columns 76 through 80 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i pN9LBDdtrd Columns 81 through 85 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i DJ8T7nHuGT Columns 86 through 90 -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i QF81D7bvUA Columns 91 through 95 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 4B7a9QFw9h Columns 96 through 100 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i ix6iFA8xoX Columns 101 through 105 -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i wt6qbkCyDE Columns 106 through 110 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i Kp5zH46zRk Columns 111 through 115 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i Yl4HdOAA61 10 / 11 名师归纳总结 - - - - - - -第 10 页,共 11 页精选学习资料 - - - - - - - - - 个人资料整理 仅限学习使用 Columns 116 through 120 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i ch4PJx4BlI Columns 121 through 125 -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i qd3YfhxCzo Columns 126 through 128 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 +64.0000i E836L11DO5 Elapsed time is 0.000120 seconds. >> 七、试验分析依据对比可知,我自己设计的C 程序中, FFT 运算时间为70 毫秒,而MATLAB中,FFT 运算时间为 0.120 毫秒；我认为运算时间其实跟算法设计是有很大关系的,我自己的设计程序算法可能仍是不是简便,MATLAB自带的 FFT 函数程序中应当是设计更优化,算法时间更简便的；但两个运算结果都是相同的,并能得到很好的与原信号频谱的拟合成效,说明它们采纳的都是 FFT 算法；S42ehLvE3M 八、试验心得此次试验由于 C 语言基础并不是特殊好,而且自从大一后就没有再使用过,所以前期构思时完全属于摸不着头脑,想了好几天；脑袋里有大致的思路,但是并不能形成系统完整的程序；于是到网上搜了许多相关资料和程序进行学习和参考,将这些程序全都弄懂并学习了相关C 语言语法后,才开头进行编程；但编程时仍旧存在心里有想法却无法转换成C 语言的问题；认为可能是C 语言不够娴熟,对蝶形运算的懂得仍不够透彻的缘由；网上大多数程序的思路是相 同的,我想这是由于 FFT 这个算法本身已经是固定算法的缘由；我从中选出一 个我懂得的特别透彻的程序进行参考,编出了相关的程序；可能方法比较笨,但已经是我努力过后的结果了；我认为虽然此次大作业比较费心血,而且才能 不是特殊好,导致程序没有调好之前的几个晚上睡前脑中仍是这些程序要怎么 调,但是通过这次大作业现在我已经对蝶形运算有了很深的懂得,并把握了整 个 FFT 算法的思路,我觉得这就是收成；我认为其实作业完成好坏无关于有没 有借鉴,重点是这其中有了收成,并使相关学问水平得到了提升,这就是作业 完成了；大作业虽然真的让我很苦恼,而且每次都要憋好久,但是我认为我在 这其中学到了许多,这些收成远比平常听课要收成的多,所以我觉得老师的大 作业是特别好的促进我们学习的方法；501nNvZFis 11 / 11 名师归纳总结 - - - - - - -第 11 页,共 11 页
