“快速傅里叶变换的实现”实验报告

《“快速傅里叶变换的实现”实验报告》由会员分享，可在线阅读，更多相关《“快速傅里叶变换的实现”实验报告（16页珍藏版）》请在金锄头文库上搜索。

1、 广东工业大学 数字图像、数字信号处理及应用数字图像、数字信号处理及应用 实验实验报告报告 题 目 快速傅里叶变换的实现 院、系(部) 信息工程学院 专业及班级 学 号 姓 名 日 期 1 1 实验目的 理解快速傅里叶变换算法的原理； 将 FFT、IFFT 算法进行基于 C 语言的实现； 利用 FFT 算法估计随机过程的功率谱（周期图法） 。 2 实验要求 （1）理解 FFT、IFFT 算法的原理及实现细节，包括复数据和实数据时的实 现区别； （2）熟悉经典功率谱估计法之一周期图法的原理和步骤； （3）完成以下设计实现： a编写 C 语言函数，完成复数据时的 FFT 及 IFFT 算法的浮点实

2、现（按时 间抽取的基 2 算法） ； b. 编写 C 语言函数， 完成单个实函数时的 FFT 及 IFFT 算法的浮点实现 （按 时间抽取的基 2 算法） ； c. 编写程序产生长度为 N 的一段均值为 0、方差为 1 的实高斯白噪声序列 v(n)，产生长度为 N 的三个实正弦序列（归一化频率分别为 0.1、0.25、0.27，信 噪比分别为 30dB、30dB、27dB） ，将三个实正弦序列和实高斯白噪声序列进行 叠加得到观测信号 uN(n)； （该步骤可通过 MATLAB 实现后将数据以文件形式保 存，N 取 32、64、128 三种情况） d. 编写 C 语言主函数测试程序，计算出 uN

3、(n)的功率谱，并用图示显示出功 率谱；（为完成图示显示， 可将功率谱数据结果以文件形式保存， 然后在 MATLAB 中导入并进行图示） e. 利用 MATLAB 检验你所编写的 FFT、IFFT 函数的正确性。 更多的改进方向： f. 编写 C 语言函数，完成 FFT、IFFT 算法的定点实现。 3 实验设备 安装了 VC、MATLAB 软件的 PC 机 4 实验原理 如何在 matlab 上实现与 C 语言的混编？ 由于本次实验中涉及到 C 语言与 matlab 进行数据交互的问题。 如果按照实验要求给出的方案来实现的话，无疑在测试中要进行频繁的文件读写才能 把数据在 matlab 和 C

4、 语言之间传递。这样，测试过程缓慢并且容易出现人 为导致的错误。好在 matlab 为我们考虑到了这点，可以使用在 matlab 使用 mex 命令对 C 语言进行编译，编译完成之后，就可以直接在 matlab 命令行调 用 C 语言程序了，避免了 matlab 与 C 语言交互过程带来的种种不便。准确 来讲使用 mex 命令是把 C 语言编译成 matlab 能够识别的动态链接库 （.mexw64/.mexw32） 。不过使用之前，用户需要解决两个问题： 2 1、如何定义 C 语言与 matlab 的接口函数及接口函数的参数意义； 2、matlab 与 C 语言是如何进行数据的传递的。 下面

5、我们就先解决第一个问题，matlab 与 C 语言的接口函数为： void mexFunction(int nlhs,mxArray *plhs,int nrhs,const mxArray *prhs); 这个函数接口可以与 main 函数接口进行类比，只不过 main 函数是操作系统 对于程序的接口， 而 mexFunction 为 matlab 调用 C 语言的接口函数。 Matlab 与 C 语言的数据交互是通过 mexFunction 的参数进行的（mexFunction 没有 返回值） 。其中， nhls 为输出数据的个数；plhs 为输出数据的指针构成的数组；nrhs 为输入 数

6、据的个数；prhs 为输入数据的指针构成的数组。 例如，使用 a,b=test(c,d,e) 调用 mex 函数 test 时，传给 test 的这三个参数分别是 prhs0=c ，prhs1=d ，prhs2=e。当函数返回时，将会把你放在 plhs0，plhs1里的地址赋给 a 和 b，达到 返回数据的目的。 第二个问题：matlab 与 C 语言是如何进行数据的传递的。 Matlab 传进 C 语言的数据是以数组的形式存在的。对于实数，可以利用 mxGetPr()函数得到 matlab 输入数据的指针。如果想输出 double 数据到 matlab， 可以使用 mxCreateDoubl

7、eMatric()函数来创建一个 double 的矩阵。 具体的输入输出函数可以通过在 matlab 键入”help mex”得到。 什么是 DFT？和 FFT 的由来。 DFT 是离散傅里叶变换的英文缩写，是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其 DTFT 的频域采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作 DFT，也应当将其看作其周期延拓的变换。 在 200 多年前法国数学家、物理学家傅里叶提出后来以他名字命名的傅里叶级数之后，用 DFT 这个工具来分析信号就已经为人们所

8、知。但在很长时间内，这种分析方法并没有引起更多的重视，最主要的原因在于这种方法运算量比较大。 FFT（快速傅里叶变换）是 1965 年由库利和图基共同提出的一种快速计算 DFT 的方法。这种方法充分利用了 DFT 运算中的对称性和周期性，从而将 DFT 运算量从 N2 减少到 N*log2N。当 N 比较小时，FFT 优势并不明显。但当 N 大于 32 开始，点数 越大，FFT 对运算量的改善越明显。比如当 N 为1024 时，FFT 的运算效率比 DFT 提 高了 100 倍。除了运算效率的大幅度提高外，FFT 还大大降低了 DFT 运算带来的累计 量化误差，这点常为人们所忽略。 什么是一次

9、实函数的快速傅里叶变换？ 在 FFT 的实际应用中，常见的情况是，需要求 FFT 的数据有实数样本组成。为了使用 FFT 算法，我们需要把实数样本置入到复数样本的实部中，并3 令复数样本的虚部为 0.计算所得的变换nF满足*FFnn。因为对于0F和2/FN，以及其他（N/2）-1 个独立值来说，这个复值数组都是实值，故它与原来实数据集有相同的自由度（N） 。然而，对实型数据采用上述方法来实现的效率是非常低的。 有两种更好的方法。 第一种是两个实函数同时变换， 俗称 “成批生产法” ，将两个单独的实数据合成输入复数数组，其中，一组实数据作为复数数组的实部，另一组作为复数数据的虚部，对复数数组进行

10、离散傅里叶变换之后，再利用离散傅里叶变换的对称性质，将两组数据的傅里叶变换分离出来。第一种方法的局限性在于必须要求输入的两个序列必须长度相同，第二种方法，一次实函数变换，有效的解决了这个问题。接下来，详细的介绍这个方法的推导过程： 设实序列为)(nx。与二次的实函数一样，先把)(nx分离为奇数序列)0()(1Nnnx和 偶数序列)0()(2Nnnx。再组成复数序列)0()()()(21Nnnjxnxny 对)(ny进行离散傅里叶变换得到： )()()()()()()()()()()()()(2121221121kjYkYkXkXjkXkXkXkjXkjXkXkjXkXkYirriiririr利

11、用离散傅里叶变换的共轭对称性： )()(2j1(k)X)()(21(k)X* 2* 1kNYkYkNYkY最后，把(k)X1、(k)X2合并成X(k) )(W(k)X)()(W(k)XX(k)2k 2N12k 2N1 kXkXkXNk 0 5 软件设计 5.1 总体设计 软件设计主要分为四个文件，四个文件各自完成自己的功能，如下： FFT1_main.cpp （复数据 FFT） 4 用户复数据输入数据顺序颠倒 （remap()函数）逐层蝶式运算数据输出IFFT_main.cpp （复数据 IFFT） 框图与相同，略 FFT2_main.cpp （实数据 FFT） 用户实数据输入实数组奇偶数分离

12、，组成复数序列调用复数FFT数据输出复数序列虚实部分离，恢复出实序 列的离散傅里叶变换IFFT2_main.cpp（实数据 IFFT） 框图与相同，略 5.2 详细设计 因为复数的 FFT 与复数的 IFFT 结构相同，所以这里只贴出 FFT 的详细过程。 数据位序颠数据位序颠倒函数倒函数 5 1（奇数）2（偶数）3（奇数）N-1（奇数）N（偶数）1（前半部分）N/2（前半部分）N/2-1（前半部分）N/2+1（后半部分）N-2（偶数） N/2+2（后半部分）N（后半部分）递归是否结束？否前半部分进行顺序颠倒后半部分进行顺序颠倒数据序列输入数据序列输出是6 int remap(Complex

13、* src, int size_n)/对数据重新排列顺序 if (size_n = 2)/函数递归调用到最后两层时停止 return 0; /*临时序列*/ Complex *temp = (Complex *)malloc(sizeof(Complex)*size_n); for (int i = 0; i= 1;layers+); remap(src, size_n);/对数据重新排列顺序 for (int i = 0; i i) % 2 = 1)/j现在代表第i层的第偶数个数据 7 Complex wn; getWN(k, size_n, /取得蝶式运算的旋转因子 Mul_Complex

14、(/乘 k += 1 #include #include #include “mex.h“ #ifndef M_PI #define M_PI 3.14159265358979323846 #endif #include “FFT.cpp“ /*matlab调用调用C语言的接口语言的接口函数函数*/ void mexFunction(int nlhs, mxArray *plhs, int nrhs, const mxArray *prhs) double *src; double *dst_imag, *dst_real; int size_n, size_n1_2; /异常处理 if(nr

15、hs!=1) mexErrMsgTxt(“USAGE: b=reverse(a)n“); if(!mxIsDouble(prhs0) mexErrMsgTxt(“the Input Matrix must be double!n“); if(mxIsComplex(prhs0) mexErrMsgTxt(“the Input Matrix can not be complex!n“); src = mxGetPr(prhs0); /取得输入指针 size_n = mxGetN(prhs0); /输入的大小 size_n1_2 = size_n/2; int layers = 0; for(int i = size_n; i = 1;layers+); if (size_n != (1 1.real = srci+1;/偶序列 src_comi1.imagin = srci;/奇序列 FFT(src_com, dst1_2N, size_n1_2);/傅里叶变换 dst1_2Nsize_n1_2.imagin = dst1_2N0.imagin;/ ds