数字信号处理实验6——利用FFT实现快速卷积

本科实验报告实验名称：数字信号处理实验实验1 利用DFT 分析信号频谱一、实验目的1.加深对DFT 原理的理解。

2.应用DFT 分析信号频谱。

3.深刻理解利用DFT 分析信号频谱的原理，分析现实过程现象及解决办法。

二、实验原理1、DFT 和DTFT 的关系有限长序列()x n 的离散时间傅里叶变换()j X e ω在频率区间(02)ωπ≤≤的N个等分点{(0),(1),(),(1)}x x x k x N-……上的N 个取样值可以由下式表示：2120()|()()01(21)N jkn j Nk k X e x n eX k k N πωωπ--====≤≤--∑由上式可知，序列()x n 的N 点DFT ()X k ，实际上就是()x n 序列的DTFT 在N 个等间隔频率点{(0),(1),(),(1)}X X X k X N -……上样本()X k 。

2、利用DFT 求DTFT方法1：由()X k 恢复出()j X eω的方法如图2.1所示：图 2.1.由 N 点DFT 恢复频谱DTFT 的流程由图2.1所示流程图可知：01()()()(22)j j nkn j n N n n k X e x n eX k W e N ωωω∞∞∞---=-∞=-∞=⎡⎤==-⎢⎥⎣⎦∑∑∑ 由式2-2可以得到12()()()(23)Nj k kx e X k N ωπφω==--∑其中()x φ为内插函数12sin()2()(24)sin()2N j N e N ωωφωω--=•-方法2：然而在实际MATLAB 计算中，上诉插值公式不见得是最好的方法。

由于DFT 是DTFT 的取样值，其相邻的两个频率样本点的间距为2Nπ,所以如果我们增加数据的长度N ，使得得到的DFT 谱线就更加精细，其包络就越接近DTFT 的结果，这样可以利用DFT 来近似计算DTFT 。

如果没有更多的数据，可以通过补零来增加数据长度。

洛阳理工学院课程设计报告课程名称_________数字信号处理_____________ 设计题目_ 用FFT实现快速卷积__ _____ 专业_________通信工程_________________ 班级__________x x x x x __________________ 学号__________xx x xx x ________________ 姓名x x xx完成日期______ 2014年6月27日_________课程设计任务书设计题目：_______用FFT实现快速卷积______________设计内容与要求：FFT的出现，使DFT在数字通信、语音信号处理、图像处理、功率谱估计、系统分析与仿真、雷达信号处理、光学、地震及数值分析等各个领域都得到广泛应用。

然而，各种应用一般都以卷积和相关运算为依据。

在实际应用中，为了分析时域离散LTI系统或者序列滤波时，需要计算两个序列的线性卷积。

为了提高运算速度，可以利用FFT来实现。

要求：参考课本上第90页的内容（3.4.1 用DFT计算线性卷积），设计并编写程序来实现重叠相加法计算线性卷积。

课程设计评语成绩：指导教师：_______________年月日目录第1章设计原理 (1)1.1设计思想 (1)1.2计算循环卷积的过程 (1)1.3重叠相加法的图示 (2)第2章设计过程 (3)2.1循环卷积子函数流程图 (3)2.2主函数流程图 (4)2.3循环卷积子函数源程序： (4)2.4主函数源程序： (6)第3章结果与验证 (9)第4章界面设计 (11)第5章分析与总结 (12)参考文献 (13)第1章 设计原理1.1设计思想运用分段处理方法中的重叠相加法计算两个序列的卷积运算。

设一个给定序列是长度为n1的A,另一个导入序列是长度为n2的B,其中B 序列是相对A 序列比较长的，所以可以把B 分为和A 一样长的若干段段，即B 分后每一小段长度为n1。

第十章 上机实验数字信号处理是一门理论和实际密切结合的课程，为深入掌握课程内容，最好在学习理论的同时，做习题和上机实验。

上机实验不仅可以帮助读者深入的理解和消化基本理论，而且能锻炼初学者的独立解决问题的能力。

本章在第二版的基础上编写了六个实验，前五个实验属基础理论实验，第六个属应用综合实验。

实验一 系统响应及系统稳定性。

实验二 时域采样与频域采样。

实验三 用FFT 对信号作频谱分析。

实验四 IIR 数字滤波器设计及软件实现。

实验五 FIR 数字滤波器设计与软件实现实验六 应用实验——数字信号处理在双音多频拨号系统中的应用任课教师根据教学进度，安排学生上机进行实验。

建议自学的读者在学习完第一章后作实验一；在学习完第三、四章后作实验二和实验三；实验四IIR 数字滤波器设计及软件实现在。

学习完第六章进行；实验五在学习完第七章后进行。

实验六综合实验在学习完第七章或者再后些进行；实验六为综合实验，在学习完本课程后再进行。

10.1 实验一: 系统响应及系统稳定性1.实验目的（1）掌握 求系统响应的方法。

（2）掌握时域离散系统的时域特性。

（3）分析、观察及检验系统的稳定性。

2.实验原理与方法在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。

已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应，本实验仅在时域求解。

在计算机上适合用递推法求差分方程的解，最简单的方法是采用MA TLAB 语言的工具箱函数filter 函数。

也可以用MATLAB 语言的工具箱函数conv 函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。

系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。

重点分析实验系统的稳定性，包括观察系统的暂态响应和稳定响应。

系统的稳定性是指对任意有界的输入信号，系统都能得到有界的系统响应。

或者系统的单位脉冲响应满足绝对可和的条件。

洛阳理工学院课程设计报告课程名称数字信号处理课程设计设计题目用FFT实现快速卷积专业通信工程班级学号姓名完成日期2015.06.15课程设计任务书设计题目：用FFT实现快速卷积设计内容与要求：FFT的出现，使DFT在数字通信、语音信号处理、图像处理、功率谱估计、系统分析与仿真、雷达信号处理、光学、地震及数值分析等各个领域都得到广泛应用。

然而，各种应用一般都以卷积和相关运算为依据。

在实际应用中，为了分析时域离散LTI系统或者序列滤波时，需要计算两个序列的线性卷积。

为了提高运算速度，可以利用FFT来实现。

要求：参考课本上第90页的内容（3.4.1 用DFT计算线性卷积），设计并编写程序来实现重叠相加法计算线性卷积。

课程设计评语成绩：指导教师：_______________年月日目录第1章概述 (1)1.1Matlab简介 (1)1.2设计目的 (2)1.3设计原理 (2)1.3.1算法产生背景 (2)1.3.2算法基本思想 (2)第2章程序设计 (5)第3章分析与测试 (7)3.1循环卷积设计 (7)3.2 线性卷积设计 (9)3.3 设计结果 (11)第4章心得体会 (12)参考文献 (13)第1章概述随着信息时代和数字世界的到来，数字信号处理已成为当今一门极其重要的学科和技术领域，数字信号处理在通信、语音、图像、自动控制、医疗和家用电器等众多领域得到了广泛的应用。

任意一个信号都具有时域与频域特性，信号的频谱完全代表了信号，因而研究信号的频谱就等于研究信号本身。

通常从频域角度对信号进行分析与处理，容易从信号的特性获得更加深入的了解。

因此，信号的频谱分析是数字信号处理技术中一种较为重要的工具。

1.1 Matlab简介Matlab语言是当今国际上科学界最具影响力、也是最有活力的软件。

它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。

它是由美国Math Works 公司于1982年推出的软件产品，取名来源于Matrix Laboratory,简称“Matlab”。

数字信号处理中的快速算法实现技巧数字信号处理是现代通信、音频处理等领域中不可或缺的重要技术。

在数字信号处理中，快速算法的实现技巧是非常关键的，它能够显著提高算法的执行效率和计算速度。

本文将介绍几种在数字信号处理中常用的快速算法实现技巧。

一、快速傅里叶变换（FFT）傅里叶变换是数字信号处理中最关键的数学工具之一。

而快速傅里叶变换（FFT）是一种高效计算傅里叶变换的方法。

它通过巧妙地利用对称性和重叠计算的思想，将原本需要O(N^2)的计算复杂度降低到O(NlogN)的复杂度，极大地提高了计算速度。

在实际的快速傅里叶变换实现中，有几个关键的技巧需要注意。

首先是蝶形运算，它是FFT中最基本的计算单元。

通过合理地组织数据的顺序，使得蝶形运算可以被高效地实现。

其次是位翻转技巧，它用于优化数据的存储和访问顺序，提高内存的访问效率。

最后是并行计算技巧，利用现代计算机的多核心特性，同时进行多个蝶形运算，进一步提高计算速度。

二、快速卷积运算卷积运算是数字信号处理中常见的操作，它用于滤波、匹配等算法中。

传统的卷积运算需要进行大量的乘法和累加运算，计算复杂度较高。

而快速卷积运算通过巧妙地利用傅里叶变换的性质，将卷积运算转化为对应频域上的乘法运算，大大降低了计算复杂度。

在实际的快速卷积运算实现中，需要注意几个关键的技巧。

首先是零填补技巧，即在进行傅里叶变换时，将信号和滤波器的长度进行零填补，使得它们的长度相等。

其次是频域乘法技巧，即在频域上进行乘法运算，再将结果进行反变换得到最终的卷积结果。

最后是快速傅里叶变换的优化技巧，在实际计算中，可以通过选择合适的傅里叶变换算法和参数，进一步提高计算速度。

三、快速运算矩阵乘法矩阵乘法是数字信号处理中常用的运算，它用于滤波器系数的更新、信号的变换等操作。

传统的矩阵乘法需要进行大量的乘法和累加运算，计算复杂度较高。

而快速运算矩阵乘法通过优化乘法的顺序和重新组织数据，将计算复杂度降低到更低的水平。

实验1 利用DFT 分析信号频谱一、实验目的1、加深对DFT 原理的理解。

2、应用DFT 分析信号的频谱。

3、深刻理解利用DFT 分析信号频谱的原理，分析实现过程中出现的现象及解决方法。

二、实验设备与环境计算机、MATLAB 软件环境。

三、实验基础理论1.DFT 与DTFT 的关系：有限长序列的离散时间傅里叶变换(e )j X ω 在频率区间(02)ωπ≤≤ 的N 个等间隔分布的点2(0k N 1)kk N πω=≤≤-上的N 个取样值可以有下式表示：2120(e )|(n)e(k)(0k N 1)N jkn j Nkk NX x X πωπω--====≤≤-∑由上式可知，序列(n)x 的N 点DFT (k)X ,实际上就是(n)x 序列的DTFT 在N 个等间隔频率点2(0k N 1)kk N πω=≤≤-上样本(k)X 。

2.利用DFT 求DTFT方法1：由(k)X 恢复出(e )j X ω的方法如下：由流程知：11(e )(n)e[(k)W]e N j j nkn j nNn n k X x X Nωωω∞∞----=-∞=-∞===∑∑∑继续整理可得到：12()(k)()Ni k kx e X N ωπφω==-∑其中(x)φ为内插函数：sin()2()sin()2N N ωφωω=方法2：实际在MATLAB 计算中，上述插值运算不见得是最好的办法。

由于DFT 是DTFT 的取样值，其相邻两个频率样本点的间距为2N π，所以如果我们增加数据的长度N ，使得到的DFT 谱线就更加精细，其包络就越接近DTFT 的结果，这样就可以利用DFT 计算DTFT 。

如果没有更多的数据，可以通过补零来增加数据长度。

3.利用DFT 分析连续信号的频谱采用计算机分析连续时间信号的频谱，第一步就是把连续信号离散化，这里需要进行两个操作：一是采样，二是截断。

对于连续时间非周期信号(t)a x ,按采样间隔T 进行采样，阶段长度M ，那么：1(j )(t)e(nT)e M j tj nTa a a n X x dt T x -∞-Ω-Ω-∞=Ω==∑⎰对(j )a X Ω 进行N 点频域采样，得到：2120(j )|(nT)e(k)M jkn Na a M kn NTX T x TX ππ--Ω==Ω==∑采用上述方法计算信号(t)a x 的频谱需要注意如下三个问题：（1）频谱混叠；（2）栅栏效应和频谱分辨率； （3）频谱泄露。

实验报告格式
图1.1 p=8，q=2 图1.2 p=8,q=4 图1.3 p=8,q=8
图1.4 q=8,p=13
图1.5 q=8,p=14
从上面的图中，根据p和q值相应的变化，可以看出：
，改变q的值：随着q的增大，经过傅里叶变化后的图像显示值的变化比较缓慢，幅度谱变化随着q的变大而变大，但相位的变化较不同。

，随着p的增大，时域信号幅值变换的比较缓慢。

图1.6 a=0.1,f=0.0625 图1.7 a=0.1,f=0.4375
图1.8a=0.1, f=0.5625
a=0.1，f=0.625，检查谱峰出现的位置是否正确，注意频谱的形式，绘制频谱特性曲线。

f=0.4375 f=0.5625，观察在这两种情况下频谱的形状和普峰出现的位置，有无混淆和泄漏现象发生？说明产生现象的发生
图1.9
FFT实现卷积运算
x1=ones(1,10);
x2=8.*sin(0.5.*n.*pi+4);
x3=0.8.*exp(3*n);
图1.10实现卷积
x1=ones(1,10);
x2=8.*sin(0.5.*n.*pi+4);
图1.11
从分别使用FFT和conv来实现卷积运算，在实验结果上可以看出，的结果是一样的。

三．一个综合性例子
图1.12
运行结果：
图1.13 正弦波图1.14正弦波图1.15 方波
图1.16 方波图1.17 三角波图1.18 三角波。

实验1 利用DFT 分析信号频谱一、实验目的1.加深对DFT 原理的理解。

2.应用DFT 分析信号的频谱。

3.深刻理解利用DFT 分析信号频谱的原理，分析实现过程中出现的现象及解决方法。

二、实验设备与环境 计算机、MATLAB 软件环境 三、实验基础理论1.DFT 与DTFT 的关系有限长序列 的离散时间傅里叶变换 在频率区间 的N 个等间隔分布的点 上的N 个取样值可以由下式表示：212/0()|()()01N jkn j Nk N k X e x n eX k k N πωωπ--====≤≤-∑由上式可知，序列 的N 点DFT ,实际上就是 序列的DTFT 在N 个等间隔频率点 上样本 。

2.利用DFT 求DTFT方法1：由恢复出的方法如下：由图2.1所示流程可知：101()()()N j j nkn j nN n n k X e x n eX k W e N ωωω∞∞----=-∞=-∞=⎡⎤==⎢⎥⎣⎦∑∑∑ 由上式可以得到：IDFTDTFT( )12()()()Nj k kX e X k Nωπφω==-∑ 其中为内插函数12sin(/2)()sin(/2)N j N x eN ωωφω--= 方法2：实际在MATLAB 计算中，上述插值运算不见得是最好的办法。

由于DFT 是DTFT 的取样值，其相邻两个频率样本点的间距为2π/N ，所以如果我们增加数据的长度N ，使得到的DFT 谱线就更加精细，其包络就越接近DTFT 的结果，这样就可以利用DFT 计算DTFT 。

如果没有更多的数据，可以通过补零来增加数据长度。

3.利用DFT 分析连续信号的频谱采用计算机分析连续时间信号的频谱，第一步就是把连续信号离散化，这里需要进行两个操作：一是采样，二是截断。

对于连续时间非周期信号,按采样间隔T 进行采样，阶段长度M ，那么：1()()()M j tj nT a a a n X j x t edt T x nT e ∞--Ω-Ω=-∞Ω==∑⎰对进行N 点频域采样，得到2120()|()()M jkn Na a M kn NTX j T x nT eTX k ππ--Ω==Ω==∑因此，可以将利用DFT 分析连续非周期信号频谱的步骤归纳如下： （1）确定时域采样间隔T ，得到离散序列（2）确定截取长度M ，得到M 点离散序列,这里为窗函数。

第十章 上机实验数字信号处理是一门理论和实际密切结合的课程，为深入掌握课程内容，最好在学习理论的同时，做习题和上机实验。

上机实验不仅可以帮助读者深入的理解和消化基本理论，而且能锻炼初学者的独立解决问题的能力。

本章在第二版的基础上编写了六个实验，前五个实验属基础理论实验，第六个属应用综合实验。

实验一 系统响应及系统稳定性。

实验二 时域采样与频域采样。

实验三 用FFT 对信号作频谱分析。

实验四 IIR 数字滤波器设计及软件实现。

实验五 FIR 数字滤波器设计与软件实现实验六 应用实验——数字信号处理在双音多频拨号系统中的应用任课教师根据教学进度，安排学生上机进行实验。

建议自学的读者在学习完第一章后作实验一；在学习完第三、四章后作实验二和实验三；实验四IIR 数字滤波器设计及软件实现在。

学习完第六章进行；实验五在学习完第七章后进行。

实验六综合实验在学习完第七章或者再后些进行；实验六为综合实验，在学习完本课程后再进行。

10.1 实验一: 系统响应及系统稳定性1.实验目的（1）掌握 求系统响应的方法。

（2）掌握时域离散系统的时域特性。

（3）分析、观察及检验系统的稳定性。

2.实验原理与方法在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。

已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应，本实验仅在时域求解。

在计算机上适合用递推法求差分方程的解，最简单的方法是采用MA TLAB 语言的工具箱函数filter 函数。

也可以用MATLAB 语言的工具箱函数conv 函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。

系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。

重点分析实验系统的稳定性，包括观察系统的暂态响应和稳定响应。

系统的稳定性是指对任意有界的输入信号，系统都能得到有界的系统响应。

或者系统的单位脉冲响应满足绝对可和的条件。

第一章1-1 有一个连续信号)2cos()(ψπ+=ft t x a ，式中Hz f 20=，2πψ=，（１） 求出)(t x a 的周期；（２） 用采样间隔s T 02.0=对)(t x a 进行采样，写出采样信号)(ˆt xa 的表达式； （３） 画出对应)(ˆt xa 的时域离散信号（序列）)(n x 的波形，并求出)(n x 的周期。

解：（1）)(t x a 的周期是s fT a 05.01==（2）∑∞-∞=-+=n a nT t fnT t x)()2cos()(ˆδψπ∑∞-∞=-+=n nT t nT )()40cos(δψπ（3）)(n x 的数字频率为πω8.0=，252=ωπ周期5=N 。

)28.0cos()(ππ+=n n x ，画出其波形如题1-1图所示。

题1-1图 1-2 设)sin()(t t x a π=，()()sin()a s s x n x nT nT π==，其中s T 为采样周期。

（1）)(t x a 信号的模拟频率Ω为多少？ （2）Ω和ω的关系是什么？（3）当s T s 5.0=时，)(n x 的数字频率ω为多少？ 解：（1）)(t x a 的模拟频率s rad /π=Ω。

（2）Ω和ω的关系是：s T ⋅Ω=ω。

（3）当s T s 5.0=时，rad πω5.0=。

1-3 判断下面的序列是否是周期的，若是周期的，确定其周期。

（1）)873cos()(ππ-=n A n x ，A 为常数；（2）)81()(π-=n j e n x 。

解： （1）πω73=，3142=ωπ，这是有理数，因此是周期序列，周期是14=T ； （2）81=ω，πωπ162=，这是无理数，因此是非周期序列。

1-4 研究一个线性时不变系统，其单位脉冲响应为指数序列)()(n u a n h n =，10<<a 。

对于矩阵输入序列，1,01()0N n N R n ≤≤-⎧=⎨⎩,其他 求出输出序列，并用MA TLAB 计算，比较其结果。

利用FFT和IFFT计算线性卷积利用FFT和IFFT计算线性卷积/**************************************************************** ********说明:用快速傅立叶变换计算两个有限长序列的线性卷积void dftCOMPLEX *x, COMPLEX *y,int n,int flagx:指向复数变量的指针变量,存放要变换的数据,长度n。

y:指向复数变量的指针变量,存放变换的结果,长度n。

n:整型变量,变换数据的长度。

flag: 整型变量,标志符,flag=1时,做DFT, flag-1时,做IDFT。

***************************************************************** *******/说明:用快速傅立叶变换计算两个有限长序列的线性卷积void convoldouble *x,double *y,int lenx,int leny,int lenx:双精度实型一维数组,长度为1en。

开始时存放实序列xi,程序结束时,存放线性卷积的结果。

y:双精度实型一维数组,长度为n。

存放实序列yi。

lenx:整型变量。

序列xi的长度。

leny:整型变量。

序列yi的长度。

len:整型变量,线性卷积的长度。

1enm+n-1,且必须是2的整数次幂/**************************************************************** ********#include “rfft.c”#include “irfft.c”#in clude “math.h”#include “graphics.h”#include “stdlib.h”#include “string.h”#include “draw.h”void signalxhdouble *x,int ch;void convolutiondouble *x,double *h,double *y,int lenx,int lenhvoid fconvoldouble *x,double *y,int lenx,int leny,int len;void mainint i,choice,select,lenx,leny,flag;double *x;// 输入信号double *h; //单位脉冲响应存放变换的结果double *y;//输出信号printf“\nchoice singal : choice??xn select??hn\n”;printf“1.矩形序列,长度10 2.单位冲击序列\n”;printf“3.单位阶跃序列,长度54.矩形序列,长度6\n”;printf“choice :选择输入信号 select:选择输出信号\n”; printf“choice”;sc anf“%d”,&choice;ifchoice1 lenx10;else ifchoice2 lenx1;else ifchoice3 lenx5;else lenx6;xdouble *calloclenx,sizeofdouble;printf“\nselect”;scanf“%d”,&select;ifselect 1 lenh10;else ifselect 2 lenh1;else ifselect 3 lenh5;else lenh6;leny lenx+lenh-1;hdouble *calloclenh,sizeofdouble;y double *callocleny,sizeofdouble;printf“\nchoice :1.continue 2.over\n”; dosingalxhx,choice;singalxhh,select;convolutionx,h,y,lenx,lenh;printf“Do you want to input length of Convolution:Y/N?”;ifgetchar’Y’||’y’printf“input length:len”;scanf“%d”,len;else len intlogleny/log2+0.9999;printf(”\nDirect Caculation of Linear Convolution\n”);fori0;ileny;i++printf“%10.1lf”,y[i];ifi%43 printf“\n”;printf“\n”;fconvolx,h,lenx,lenh,len;printf(”Fast Caculation of Linear Convolution\n”);fori0;ileny;i++printf“%10.1lf”,x[i];ifi%43 printf“\n”;printf“\n”;scanf“&d”,&flag;while!flag;void singalxhdouble *x,int chint i;double al,a,f,t;printf“1.矩形序列,长度10 2.单位冲击序列\n”;printf“3.1,2.4,-2.4,2.4,长度44.矩形序列,长度6\n”;switchchcase 1: fori0;i10;i++ x[i]1;break;case 2: x[0]1;break;case 3: x[0]1;x[1]2.4;x[2]-2.4;x[3]2.4;break;default: fori0;i5;i++ x[i]1;break;void convolutiondouble *x,double *h,double *y,int lenx,int lenh int len;int i,j,k;lenlenx+lenh-1;fork0;klen;k++y[k]0.0;fori0;ik;i++y[k]y[k]+x[i]*h[k-i];void fconvoldouble *x,double *y,int lenx,int leny,int lenint i,len2;double t;forforilenx;ilen;i++ x[i]0.0;forilenh;ilen;i++ y[i]0.0;rfftx,len;rffty,len;len2len/2;x[0]x[0]*y[0];x[len2]x[len2]*y[len2];fori1;ilen2;i++tx[i]*y[i]-x[len-i]*y[len-i];x[len-i] x[i]* y[len-i]+ x[len-i]* y[i]; x[i]t;irfftx,len;void rfftdouble *x,int nint i,j,k,m,i1,i2,i3,i4,nl,n2,n4;double a,e,cc,ss,xt,t1,t2;forj1,i1;i16;i++m=i;j=2*j;ifjn break;n1=n-1;forj0,i0;i<n;i++=ifi<jxtx[j];x[j]x[i];x[i]=x[j];kn/2;whilekj十1j=j?k;kk/2;j=j+k;fori=0;i<n;i+2=xtx[i];x[i]xt+x[i+1];x[i+1]=xt-x[i+1j]; n21;for k=2;km;k++n4=n2;n2=2*n4;n12*n2;e=6.28318530718/n1;fori=0;in;i+n1xtx[i];x[i]xt+x[i+n2];x[i+n2]=xt-x[i+n2];x[i+n2+n4]-x[i+n2+n4]; ae;forj1;jn4?1;j++i1=i+j;i2i-j+n2;i3=i+j+n2;i4i-j+n1;cccosa;ss=sina;aa+e;t1=cc*x[i3]+ss*x[i4];t2=ss*x[i3]-cc*x[i4];x[i4]=x[i2]-t2;x[i3]=-x[i2]-t2;x[i2]=x[i1]-t1;x[i1]=x[i1]+t;void irfftdouble *x,int nint i,j,k,m,i1,i2,i3,i4,i5,i6,i7,i8,n2,n4,n8,id,is; double a,e,a3,t1,t2,t3,t4,t5,cc1,cc3,ss1,ss3;forj1,i1;i16;i++m=i;j=2*j;ifjn break;n2=2*n;fork1;km;k++is0;id=n2;n2=n2/2;R4n2/4;n8n2/8;e=6.28318530718/n2; doforiis;in;i+idi1=i;i2i1+n4;i3i2+n4;i4i3+n4;tl=x[i1]-x[i3];x[i1]x[i1]+x[i3]; x[i2]2*x[i2];x[i3]t1-2*x[i4];x[i4]t1+2*x[i4];ifn41 continue;i1+n8;i3+n8;i4+n8;t1x[i2]-x[i1]/sqrt2.0; t2x[i4]+x[i3]/sqrt2.0; x[i1]x[i1]+x[i2];x[i2]x[i4]-x[i3];x[i3]2*-t2-t1;x[i4]2*-t2+t1;is2*id-n2;id4*id;whileisn-1;ae;forj=1;jn8;j++a3=3*a;cc1cosa;ss1sina;cc3=cosa3;ss3=sina3;a=j+1*e;id2*n2;doforiis;in-1;ii+idi1i-j;i2=i1+n4;i3=i2+n4;i4=i3+n4;15=i+n4-j;i6=i5+n4;i7=i6+n4;i8=i7+n4;t1=x[i1]-x[i6];x[i1]=x[i1]+x[i6]; t2=x[i5]-x[i2];x[i5]=x[i2]+x[i5]; t3=x[i8]+x[i3];x[i6]=x[i8]-x[i3]; t4=x[i4]+x[i7];x[i2]=x[i4]-x[i7];t5t1-t4;t1t1+t4;t4=t2-t3;t2=t2+t3;x[i3]t5*cc1+t4*ss1; x[i7]-t4*cc1+t5*ss1; x[i4]t1*cc3-t2*ss3; x[i8]t2*cc3+t1*ss3;is2*id-n2;id=4*id;whileisn-1;is=0;id=4;doforiis;in;i+idi1=i+1;t1=x[i];x[i]=t1+x[i1];x[i1]=t1-x[i1];is2*id-2;id4*id;whileisn-1forj=0,i0;in-1;i++ ifijt1=x[j];x[j]=X[i];x[i]=t1;k=n/2;whi1ekj+1j=j-k;k=k/2;j=j十k;fori=0;in;i++ x[i]x[i]/n;。

实验课程名称数字信号处理实验工程名称利用FFT实现快速卷积专业电信年级2008级姓名学号所在学院应用科技学院2009年月曰实验准备实验目地(1＞加深理解FFT在实现数字滤波＜或快速卷积)中地重要作用，更好地利用FFT进行数字信号处理(2＞掌握循环卷积和线性卷积两者之间地关系实验环境＜1)计算机一台＜2) Matlab6.5或以上版本实验内容一•实验原理数字滤波器根据系统地单位脉冲响应h( n＞是有限长还是无限长可分为有限长单位脉冲响应＜Finite Impulse Response)系统＜简记为 FIR系统)和无限长单位脉冲响应＜lnfinite Impulse Respons©系统＜简记为IIR系统).对于FIR滤波器来说,除了可以通过数字网络来实现外，也可以通过FFT地变换来实现.首先我们知道，一个信号序列x(n＞通过FIR滤波器时，其输出应该是x(n＞与h(n＞地卷积：-bey(n) = x(n)*h(n) = ' x(m)h(n - m)m =.::：:当h(n＞是一个有限长序列，即h(n＞是FIR滤波器，且0乞n乞N -1时N 4y(n)二 ' h(m)x(n -m)m £在数字网络类地FIR滤波器中，普遍使用地横截型结构就是按这个卷积公式构成地应用FFT实现数字滤波器实际上就是用FFT来快速计算有限长度序列地线性卷积 .这种方法就是先将输入信号x(n＞通过FFT变换为它地频谱采样值X(k＞,然后再和 FIR滤波器地④完成 X(k> 和 H(k>乘积，丫(k) = X(k)H(k)⑤用FFT 计算Y(k ＞地离散傅里叶反变换得结果,然后用相应地方式把它们结合起来 ，便是总地输出.分段卷积方法主要有两种，即重叠相加法和重 叠保留法.具体内容请参考数字信号处理教材中“快速离散傅里叶变换” 一章中地线性卷积地 法部分,本实验这部分不作重点要求 二、实验内容1 .数字滤波器地脉冲响应为 h(n) =(1/2 $ R N 2(n) ,N 2取8. 输入序列x(n ＞可选下列几种情况2•实验前,预先计算好x(n)* h(n)地值.3 •实验前，预先编制一个应用 FFT 实现数字滤波器地通用程序 .＜或用Matlab 提供地FFT 函数)① x(n) ② x(n) ③ x(n)= R N 1(n),N 1 取 8. 2砧 -cos — N 1 3 R M (n)N 取 8.R N ®) ,N I 取 8. 频响采样值 H(k ＞相乘,H(k ＞可事先存放在存储器中，最后再将乘积H(k ＞X(k ＞通过快速傅里叶变换 ＜简称IFFT )还原为时域序列，即得到输出y(n ＞. 现以FFT 求有限长序列间地卷积及求有限长度序列与较长序列间地卷积为例来讨论FFT 地快速卷积方法.a.序列x(n ＞和h(n ＞地长差不多.设x(n ＞地长为N 1,h(n ＞地长为N 2，要求N J y(n) = x(n) ： y(n) 八 h(m)x(n _ m)m -0用FFT 完成这一卷积地具体步骤如下：①为使两有限长序列地线性卷积可用其循环卷积代替而不发生混叠 ,必须选择循环卷积长度 N _N 1 • N 2 -1,若采用基2-FFT 完成卷积运算 要求N =2m＜m 为整数).②用补零方法使 'x(n)x(n ＞和h(n ＞变成列长为 N 地序列. h(n)③用FFT 计算 x(n ＞和h(n ＞地N 点离散傅里叶变换y(n )=迟［*Y(k)k=0 !_N =柩［丄Y *(k)W 一b.当x(n ＞长度很长时，即N 1 • • N 2,通常不允许等x(n ＞全部采集齐后再进行卷积，否则使输出相 对于输入有较长地延时，另外，若N 1 • N 2 -1太大,h(n ＞要补上太多地零点，很不经济，且FFT 地计算时 间也要很长.为此，采用分段卷积地方法，即把x(n ＞分成长度与h(n ＞相仿地一段段，分别求出每段卷积地_nk N nk N FFT 算请同学们做好实验后，以格式为“实验四+学号+姓名”地形式命名文档＜例如：实验四120602008001蔡和颜.doc).实验报告提交方式：实验报告以电子文档地形式提交给学习委员学习委员统一交给老师.【参考资料一】1.快速傅立叶变换(FFT＞算法长度为N地序列x(n)地离散傅立叶变换X(k)为:N -4X(k)八x(n)W「k,k =0,•…,N -1n =eN点地DFT可以分解为两个 N/2点地DFT，每个N/2点地DFT又可以分解为两个，由N/4点,所以地DFT.依此类推，当N为2地整数次幕时(N = 2M＞,由于每分解一次降低一阶幕次通过M次地分解，最后全部成为一系列 2点DFT运算.以上就是按时间抽取地快速傅立叶变换(FFT>算法.当需要进行变换地序列地长度不是2地整数次方地时候，为了使用以2为基地FFT,可以用末尾补零地方法，使其长度延长至 2地整数次方.b5E2RGbCAP序列X(k)地离散傅立叶反变换为N Jnkx(n) X(k)W N, n =0,•…，N-1N y离散傅立叶反变换与正变换地区别在于W N变为W N1,并多了一个1 N地运算.因为W N和1W N对于推导按时间抽取地快速傅立叶变换算法并无实质性区别，因此可将FFT和快速傅立叶反变换<IFFT )算法合并在同一个程序中.p1EanqFDPw2.利用FFT进行频谱分析若信号本身是有限长地序列，计算序列地频谱就是直接对序列进行FFT运算求得X(k), X (k)就代表了序列在0,2二］之间地频谱值.DXDiTa9E3d幅度谱X(k) r J x R(k) • X2(k)XI (k)相位谱「(k)二arctanX R(k)若信号是模拟信号,用FFT进行谱分析时，首先必须对信号进行采样，使之变成离散信号然后就可按照前面地方法用FFT来对连续信号进行谱分析.按采样定理，采样频率f s应大于2倍信号地最高频率，为了满足采样定理，一般在采样之前要设置一个抗混叠低通滤波器.用FFT对模拟信号进行谱分析地方框图如下所示.< 注意：采样信号地频谱要乘以Ts,其主值区间才代表原来模拟信号地频谱)RTCrpUDGiT【参考资料二】1、 fft功能：一维快速傅里叶变换<FFT)调用格式：y=fft(x,n> :采用n点FFT.当x地长度小于 n时,fft函数在x地尾部补零，以构成n点数据；当x地长度大于n时,fft函数会截断序列 x.当x为矩阵时,fft函数按类似方式处理列长度.5PCzVD7HxA2、 fftshift功能：对 fft 地输出进行重新排列 ,将零频分量移到频谱中心调用格式：y=fftshift(x> ：对 fft 地输出进行重新排列 ,将零频分量移到频谱中心当 x 为向量时 ,fftshift(x> 直接将 x 中地左右两半交换而产生 y. 当 x 为矩阵时 ,fftshift(x> 同时将 x 地左右、上下进行交换而产生 y。

实验三利用FFT计算线性卷积实验记录1.比较当序列长度分别为8,16,32,64,256,512,1024时两种方法计算的时间一、直接线性卷积程序代码for n1=3:10;L=2^n1;x=ones(1,L);h=cos(0.2*pi.*x);ticy=conv(x,h);tocend程序结果：Elapsed time is 0.000104 seconds.Elapsed time is 0.000022 seconds.Elapsed time is 0.000018 seconds.Elapsed time is 0.000021 seconds.Elapsed time is 0.000031 seconds.Elapsed time is 0.000059 seconds.Elapsed time is 0.000152 seconds.Elapsed time is 0.000508 seconds.二、快速卷积程序代码for n0=3:10;L=2^n0;n=0:L;x=heaviside(n)-heaviside(n-L);h=cos(0.2*pi.*n);X=fft(x);H=fft(h);Y=X.*H;ticy=ifft(Y);tocend程序结果Elapsed time is 0.000041 seconds. Elapsed time is 0.000441 seconds. Elapsed time is 0.000075 seconds. Elapsed time is 0.000070 seconds. Elapsed time is 0.001507 seconds. Elapsed time is 0.000544 seconds. Elapsed time is 0.005577 seconds. Elapsed time is 0.002435 seconds.2.考察当L=2048和4096，M=256时两种方法计算时间1）N=2048时一、直接线性卷积程序代码L=2048;M=256;n0=0:L;x=heaviside(n0)-heaviside(n0-L);n1=0:M;h=cos(0.2*pi.*n1);ticy=conv(x,h);toc程序结果Elapsed time is 0.011091 seconds二、快速卷积程序代码L=2048;M=256;n0=0:L;x=heaviside(n0)-heaviside(n0-L);n1=0:M;h=cos(0.2*pi.*n1);ticX=fft(x,2048);H=fft(h,2048);Y=X.*H;y=ifft(Y);toc程序结果Elapsed time is 0.005779 seconds.2）N=4096时一、直接线性卷积程序代码L=4096;M=256;n0=0:L;x=heaviside(n0)-heaviside(n0-L); n1=0:M;h=cos(0.2*pi.*n1);ticy=conv(x,h);toc程序结果Elapsed time is 0.002153 seconds.二、快速卷积程序代码L=4096;M=256;n0=0:L;x=heaviside(n0)-heaviside(n0-L); n1=0:M;h=cos(0.2*pi.*n1);ticX=fft(x,2048);H=fft(h,2048);Y=X.*H;y=ifft(Y);toc程序结果Elapsed time is 0.005149 seconds.3.利用重叠相加法计算线性卷积，考察L=2048和4096时的计算时间1）L=2048程序代码L=2048;M=256;x=ones(1,L);n=0:M-1;h=cos(0.2*pi*n);N=L+M-1;y=zeros(1,N);subLen=M;subNum=L/M;subLenY=2*M-1;ticH=fft(h,2*M-1);for count=1:subNumx_temp=x(((count-1)*subLen+1):(count*subLen));X=fft(x_temp,2*M-1);Y=X.*H;y_temp=ifft(Y,2*M-1);y(((count-1)*subLen+1):(count*subLen+M-1))=y(((count-1)*subLen+1):(count*subLen +M-1))+y_temp;endtoc程序结果Elapsed time is 0.003825 seconds.2）L=4096程序代码L=4096;M=256;x=ones(1,L);n=0:M-1;h=cos(0.2*pi*n);N=L+M-1;y=zeros(1,N);subLen=M;subNum=L/M;subLenY=2*M-1;ticH=fft(h,2*M-1);for count=1:subNumx_temp=x(((count-1)*subLen+1):(count*subLen));X=fft(x_temp,2*M-1);Y=X.*H;y_temp=ifft(Y,2*M-1);y(((count-1)*subLen+1):(count*subLen+M-1))=y(((count-1)*subLen+1):(count*subLen +M-1))+y_temp;endtoc程序结果Elapsed time is 0.006678 seconds.4.利用重叠保留法计算两个序列的线性卷积。

实验1 利用DFT 分析信号频谱一、实验目的1.加深对DFT 原理的理解。

2.应用DFT 分析信号的频谱。

3.深刻理解利用DFT 分析信号频谱的原理，分析实现过程中出现的现象及解决方法。

二、实验设备与环境 计算机、MATLAB 软件环境 三、实验基础理论1.DFT 与DTFT 的关系有限长序列x (n )(0≤n ≤N −1)的离散时间傅里叶变换X (e jω)在频率区间(0≤ω≤2π)的N 个等间隔分布的点kω=2πk /N (0≤k ≤N −1)上的N 个取样值可以由下式表示：212/0()|()()01N jkn j Nk N k X e x n eX k k N πωωπ--====≤≤-∑由上式可知，序列x (n )的N 点DFT X k ,实际上就是x (n )序列的DTFT 在N 个等间隔频率点kω=2πk /N (0≤k ≤N −1)上样本X k 。

2.利用DFT 求DTFT方法1：由恢复出的方法如下：由图2.1所示流程可知：101()()()N j j nkn j nN n n k X e x n eX k W e N ωωω∞∞----=-∞=-∞=⎡⎤==⎢⎥⎣⎦∑∑∑ 由上式可以得到：IDFTDTFTX (ejω)12()()()Nj k kX e X k Nωπφω==-∑ 其中为内插函数12sin(/2)()sin(/2)N j N x eN ωωφω--= 方法2：实际在MATLAB 计算中，上述插值运算不见得是最好的办法。

由于DFT 是DTFT 的取样值，其相邻两个频率样本点的间距为2π/N ，所以如果我们增加数据的长度N ，使得到的DFT 谱线就更加精细，其包络就越接近DTFT 的结果，这样就可以利用DFT 计算DTFT 。

如果没有更多的数据，可以通过补零来增加数据长度。

3.利用DFT 分析连续信号的频谱采用计算机分析连续时间信号的频谱，第一步就是把连续信号离散化，这里需要进行两个操作：一是采样，二是截断。

第十章上机实验数字信号处理是一门理论和实际密切结合的课程，为深入掌握课程内容，最好在学习理论的同时，做习题和上机实验。

上机实验不仅可以帮助读者深入的理解和消化基本理论，而且能锻炼初学者的独立解决问题的能力。

本章在第二版的基础上编写了六个实验，前五个实验属基础理论实验，第六个属应用综合实验。

实验一系统响应及系统稳定性。

实验二时域采样与频域采样。

实验三用FFT对信号作频谱分析。

实验四IIR数字滤波器设计及软件实现。

实验五FIR数字滤波器设计与软件实现实验六应用实验——数字信号处理在双音多频拨号系统中的应用任课教师根据教学进度，安排学生上机进行实验。

建议自学的读者在学习完第一章后作实验一；在学习完第三、四章后作实验二和实验三；实验四IIR数字滤波器设计及软件实现在。

学习完第六章进行；实验五在学习完第七章后进行。

实验六综合实验在学习完第七章或者再后些进行；实验六为综合实验，在学习完本课程后再进行。

10.1 实验一: 系统响应及系统稳定性1.实验目的（1）掌握求系统响应的方法。

（2）掌握时域离散系统的时域特性。

（3）分析、观察及检验系统的稳定性。

2.实验原理与方法在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。

已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应，本实验仅在时域求解。

在计算机上适合用递推法求差分方程的解，最简单的方法是采用MA TLAB语言的工具箱函数filter函数。

也可以用MATLAB语言的工具箱函数conv函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。

系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。

重点分析实验系统的稳定性，包括观察系统的暂态响应和稳定响应。

系统的稳定性是指对任意有界的输入信号，系统都能得到有界的系统响应。

或者系统的单位脉冲响应满足绝对可和的条件。

系统的稳定性由其差分方程的系数决定。

数字信号处理实验——快速傅里叶变换一、概述数字信号处理是一门涉及数字信号采集、处理和分析的学科。

在数字信号处理中，傅里叶变换是一个非常重要的工具，可以将时域信号转换为频域信号，帮助我们理解信号的频谱特性。

而快速傅里叶变换（FFT）是一种高效的计算傅里叶变换的算法，广泛应用于信号处理领域。

本文将介绍数字信号处理实验中快速傅里叶变换的基本原理、算法实现和实验过程。

二、快速傅里叶变换的原理快速傅里叶变换是一种通过递归分治策略减少计算量的傅里叶变换算法。

对于长度为N的离散信号序列，其快速傅里叶变换可以分解为若干个长度为N/2的子序列的快速傅里叶变换的线性组合。

通过这种分治的方法，可以将原始的傅里叶变换计算复杂度从O(N^2)降低到O(NlogN)。

三、快速傅里叶变换的算法实现快速傅里叶变换的算法实现主要分为递归方法和迭代方法两种。

递归方法是一种自顶向下的计算方法，通过将长度为N的信号序列分解为两个长度为N/2的子序列，并利用子序列的快速傅里叶变换计算原始序列的傅里叶变换。

迭代方法则是一种自底向上的计算方法，通过不断合并较短序列的傅里叶变换结果来计算较长序列的傅里叶变换。

在实际应用中，迭代方法通常比递归方法更加高效，特别是对于长度为2的幂次方的信号序列。

四、数字信号处理实验在数字信号处理实验中，快速傅里叶变换通常作为一种重要的信号处理工具使用。

实验过程中，我们首先需要准备一段离散的时域信号序列，然后利用快速傅里叶变换算法将其转换为频域信号序列。

通过对频域信号序列的分析，我们可以获取信号的频谱特性，包括频率成分、谐波分量等信息。

五、实验流程1. 准备一段离散的时域信号序列，可以是从传感器获取的实际物理信号，也可以是由数学函数生成的合成信号。

2. 对时域信号序列进行快速傅里叶变换，得到频域信号序列。

3. 分析频域信号序列，获取信号的频谱特性，如主要频率成分、谐波分量等。

4. 对频域信号序列进行反变换，将其恢复为时域信号序列进行进一步分析。

利用傅里叶变换实现时域的卷积【摘要】本文将探讨利用傅里叶变换实现时域的卷积的原理和应用。

首先介绍傅里叶变换的基本概念和原理，然后详细解释在实现时域卷积中如何利用傅里叶变换来加速计算，并给出具体的实例和应用。

最后总结回顾，强调傅里叶变换在时域卷积中的重要性和优势。

【正文】一、傅里叶变换的基本概念和原理傅里叶变换是将一个函数在频域中表示的数学工具。

它将时域上的信号分解为一系列基本的频率成分，从而能够更好地理解信号的频谱特征。

傅里叶变换可以将一个连续时域函数f(t)转换成一个连续频率域函数F(ω)，其中ω表示频率。

对于离散时间信号或离散频率信号，同样可以利用离散傅里叶变换（DFT）进行频谱分析。

二、利用傅里叶变换实现时域的卷积原理在信号处理领域，时域卷积是一种常用的运算，用于处理信号的滤波、降噪等应用。

然而，在时域中进行卷积运算通常涉及到大量的计算，尤其是对于长时域信号而言，计算复杂度非常高。

利用傅里叶变换可以实现对时域卷积的高效计算。

基本思路是将时域上的卷积运算转换到频域中进行，然后再通过逆傅里叶变换将结果转换回时域。

具体步骤如下：1. 对待卷积的两个时域信号f(t)和g(t)分别进行傅里叶变换，得到它们在频域上的表达F(ω)和G(ω)。

2. 将F(ω)和G(ω)逐元素相乘，即得到结果的频域表示H(ω)。

3. 对H(ω)进行逆傅里叶变换，得到时域上的卷积结果h(t)。

三、利用傅里叶变换实现时域的卷积应用利用傅里叶变换实现时域的卷积广泛应用于信号处理和图像处理领域。

下面以图像处理为例，介绍其具体应用。

图像卷积是一种常见的图像处理操作，用于实现图像的模糊、锐化和滤波等效果。

传统的图像卷积通常需要对图像的每个像素点及其周围像素点进行计算。

对于大尺寸图像而言，计算复杂度非常高。

利用傅里叶变换可以将图像卷积转换成频域上的元素相乘运算，大大加快了计算速度。

具体步骤如下：1. 对待卷积的图像和卷积核分别进行傅里叶变换，得到它们在频域上的表达。