matlab-数字信号处理.doc

目录实验一 学习使用MATLAB1一.实验目的1二.实验内容和要求1三.实验主要仪器设备和材料1四.实验方法、步骤及结果测试1五.实验报告要求24六.思考题25实验二 用FFT做谱分析26一.实验目的26二.实验内容和要求26三.实验主要仪器设备和材料27四.实验方法、步骤及结果测试27五.实验报告要求47六.思考题47实验三 用双线性变换法设计IIR数字滤波器48一.实验目的48二.实验内容和要求48三.实验主要仪器设备和材料48四.实验方法、步骤及结果测试48五.实验报告要求53六.思考题53实验四 用窗函数法设计FIR数字滤波器55一.实验目的55二.实验内容和要求55三.实验主要仪器设备和材料56四.实验方法、步骤及结果测试56五.实验报告要求66六.思考题66附录一 MATLAB信号处理工具箱函数67一.表附1-1 波形产生67二.表附1-2 滤波器分析和实现67三.表附1-3 线性系统变换68四.表附1-4 IIR滤波器设计69五.表附1-5 IIR滤波器阶的选择69六.表附1-6 FIR滤波器设计70七.表附1-7 变换70八.表附1-8 统计信号处理71九.表附1-9 窗函数71十.表附1-10 参数化建模72十一.表附1-11 特殊操作72十二.表附1-12 模拟原形滤波器设计73十三.表附1-13 频率变换73十四.表附1-14 滤波器离散化73十五.表附1-15 其他74附录二 实验中用到的一些子程序75实验内容和要求(1) 了解MATLAB 的基本程序设计原则，常量和变量的用法(2) 掌握MATLAB中对矩阵进行输入、运算和比较的方法(3) 了解循环语句的类型，并掌握循环语句的用法(4) 熟悉M文件的作用，并掌握二维图形的绘制一. 实验主要仪器设备和材料计算机，MATLAB6.5或以上版本二. 实验方法、步骤及结果测试关于MATLAB它是由美国的Math Works 公司推出的一个科技应用软件，它的名字是由MATRIX(矩阵）和LABORATORY（实验室）的前三个字母组合而成 MATLAB是一种高性能的、用于工程计算的编程软件，它把科学计算、结果的可视化和编程都集中在一个使用方便的环境中优势在于能很容易求解复数 数值问题，速度快且容易扩展创建新的命令和函数主要组成部分：（1）编程语言：以矩阵和数组为基本单位的编程语言（2）工作环境：包括一系列的应用工具，提供编程和调试程序的环境（3）图形处理：包括绘制二维、三维图形和创建图形用户界面（GUI）等（4）数学库函数：包含大量数学函数，也包括复杂功能（5）应用程序接口：提供接口程序，可使MATLAB与其他语言程序进行交互典型特点：（1）语言简洁紧凑，运算符十分丰富，使用方便灵活（2）既具有结构化的控制语言，又能面向对象编程（3）语法限制不严格，程序设计自由度大，可移植性好（4）具有强大的图形功能（5）包含功能强劲的工具箱（6）最重要、最受欢迎的特点是它的开放性（7）典型应用有两方面： 数值计算和符号计算 建模和动态仿真下面介绍MATLAB的界面、常用命令和使用方法图1-1 MATLAB窗口菜单栏和工具栏：位于窗口顶部，用户可以通过它们来执行某些命令命令窗口：位于右边空白部分，用户的数据输入和结果运算，都在此窗口进行，是 Matlab 极为重要的部分，也是用户使用最频繁的部分工作台和工具箱：位于主窗口左上部分，双击工具箱或前面的”+”号，就能看到工具箱的各项功能工作空间：主窗口的中上部分，可看到 Matlab 的各个工作变量，新打开 Matlab 时，只能看到系统提供的默认输出变量ans历史命令：主窗口的左下部分，主要保存工作过的变量、表达式等，需要时，用户可以直接提取历史命令在命令窗口中使用当前工作目录：主窗口的中下部分，主要保存在当前工作路径下的图形文件和命令文件二、常量和变量MATLAB 中使用的数据有常量和变量作用标量的实数常量，类似于 C 语言中的整形常量和实形常量，常量的表示方法：十进制数字，可以带正负号和小数点，也可采用科学计数法如： 1, 2.5 , 0.0033 , 2 e-7 ， pi , 2+3 i 等变量以其名称在操作语句中第一次合法出现而定义，无需事先定义。如果变量已经存在，那么 matlab 将改变它的内容。变量命名原则：（ 1 ）变量名必须以字母开头，后面可以是任意的字母、数字或下划线（ 2 ）变量名有效的最大长度是 31（ 3 ）变量要区分大小写特殊的常量和变量（ 1 ） ans 缺省变量名，以操作中最近应答作为它的值（ 2 ） pi 圆周率（ 3 ） i 或 j 虚数单位（ 4 ） INF 正无穷大（ 5 ） nan 不定值（或称非数）（ 6 ） version MATLAB的版本三、程序的运算符1 、算术运算符：优先级最高加法：+ 减法： 乘法：* 除法：/ 左除： 乘方：执行时的优先法则按照以下规则：（ 1 ）表达式从左到右执行（ 2 ）乘方运算具有最高的优先级（ 3 ）乘法和除法运算具有相同的次优先级，加法和减法具有相同的最低优先级（ 4 ）括号可以改变优先次序，由最内层的括号向外执行例1.1: A=1 2 3;4 5 6;7 8 9; b=10; A+bans = 11 12 13 14 15 16 17 18 19 B=ones(3); % 输入一个3*3阶的全为1的矩阵 A-Bans = 0 1 2 3 4 5 6 7 8 S=A*B % 矩阵的乘法S = 6 6 6 15 15 15 24 24 24 r=3; U=Ar % 矩阵的乘方U = 468 576 684 1062 1305 1548 1656 2034 2412X=2 3;1 3;Y=1 4;7 9;Z=2 5;6 8;W=2* X+ Y3- Z %按顺序执行W= 311 477 829 1259W=(2* X+Y)3- Z %先计算括号，再计算乘方，最后减去 ZW= 2373 4145 3729 6517W=2*(X+Y)3- Z %先计算括号，再计算乘方，然后乘以2，最后减去ZW= 2068 3425 3914 6472关系运算符：对程序的流程控制非常有用，共 6 个 小于 大于 小于等于 = 大于等于 = 等于 = 不等于 （=）注：所有关系运算符连接的两个操作数必须或者同维、同大小，或者其中一个为标量。有一个操作数是标量时，关系运算符把标量和另一个操作数的每个元素进行比较。对大小相同的两个矩阵，进行关系运算时，是对相应的每一个元素进行比较，如果满足指定关系，则返回 1 ，否则返回 0 。例1.2: A=1 2;3 4 % 输入一个矩阵A = 1 2 3 4 B=1,0;3,5 % 输入一个矩阵B = 1 0 3 5 A A2 %比较矩阵中的每个元素与2的大小关系ans = 0 0 1 1 C=1 2 3;3 2 1;2 1 3; D=2 2 3;3 3 3;1 1 1; C=D %比较两个矩阵中的每个元素是否相等，相同为1，不同为0ans = 0 1 1 1 0 0 0 1 0 C=D %比较两个矩阵中的每个元素是否相等，不同为1，相同为0ans = 1 0 0 0 1 1 1 0 1 C=D ans = 0 1 1 1 0 0 1 1 1 CA=1 0;2 3; %输入一个矩阵B=1 1;2 2; %输入一个矩阵A&B %与运算ans = 1 0 1 1A| B %或运算ans = 1 1 1 1 A %非运算ans = 0 1 0 0XOR(A,B) %逻辑异或ans = 0 1 0 0四、程序控制语句：为编写结构化的程序提供了必不可少的条件，不但可以方便用户编程，还可提高程序的效率和可读性1、 for循环语句：用来执行循环次数已知的情况调用格式： for x=初值：步长：终值 循环体 end其中：变量 x 称为循环变量，初值、终值和步长可以是标量，也可以是表达式当循环语句开始执行时， x 的值被赋为和初值相同的内容，每执行一次循环体的内容， x 的值就会按照步长的大小来改变，如果步长为正数，每执行一次就增加一个步长，否则减小一个步长，一直到变量的值大于或者小于终值，for语句循环结束，继续执行结束语句下面的命令。注：当步长省略不写时，系统默认为 1 。 for循环语句允许嵌套；每一个 “for 关键字必须和一个 “end 配对例1.4: 下面for语句中的循环体将执行 100 次，计算 1 到 100 的和 a=0; for i=1:100a=a+i;end aa = 5050for循环语句可利用嵌套方式使用，如下例计算 1 到 6 的乘法表 for i=1:6for j=1:iA(i,j)=i*j;endend AA = 1 0 0 0 0 0 2 4 0 0 0 0 3 6 9 0 0 0 4 8 12 16 0 0 5 10 15 20 25 0 6 12 18 24 30 362、while循环语句：一般用于事先不能确定循环次数的情况时候调用格式： while 表达式 循环体 end当表达式的值为真时，执行循环体，当循环体执行完毕后，继续判断表达式的值，如果仍为真，继续执行，如此循环，直到表达式的值为假时终止循环注：利用break语句，可在任何时候终止循环（while循环或者for循环）例1.5：求 10 的阶乘 s=1; n=1; while (n ss = 36288003、条件选择语句if根据不同逻辑表达式的值来判断程序在执行过程中需要执行那些语句调用格式： if 逻辑表达式 命令语句体 end程序在执行过程中，首先判断逻辑表达式的值，若为真，则执行命令语句体中所有语句，否则，跳过命令语句体，执行结束语句之后的语句if语句可以通过嵌套的方式使用，但不同的 if 语句之间不能互相交叉调用格式： if 逻辑表达式 1 语句体 1 elseif 逻辑表达式 2 语句体 2 elseif 逻辑表达式 3 语句体 3 else 语句体 end例1.6： a=100;b=345; if a cc = 3.4500例1.7： X=1 2 3; 4 5 6; 7 8 9; i=2; j=3; if i=jX(i,j)=0;elseif abs(i-j)=1X(i,j)=-1;else X(i,j)=10;end XX = 1 2 3 4 5 -1 7 8 9五、M 文件为代替在 MATLAB 提示符下输入的 MATLAB 命令语句，可以把这些命令写入一个文本文件，每当用户输入文件名和它的自变量时，命令语句就由 MATLAB 执行，当文件中最后一个命令被执行后， MATLAB 能再从终端读取命令，这样的文件就叫 M 文件1、M 文件的类型两类：命令（Script）文件和函数（function）文件（ 1 ）、命令文件主要用途：使命令输入更简单化（没有输入参数也没有输出参数）当用户需要重复输入许多相同的命令时，可将它们放在一个命令文件中，每次只要输入文件名，即可得到相同的运行结果实质是将用户在 MATLAB 命令窗口中输入的一串命令用另外一个名称来代替（ 2 ）、函数文件主要用途：扩充 MATLAB 的应用范围和满足用户不同的实际应用需求函数文件可接受输入变量也可返回输出变量函数文件对变量工作的运算是在其本身的工作空间内进行的，它和 MATLAB 的命令窗口的工作空间是相互独立的2、M 文本编辑器：利用 MATLAB 编程语言所写的程序被称作 M 文件后缀均为.m（ 1 ）启动编辑器的方法：（ A ）直接在 MATLAB 命令窗口中运行命令编辑（ B ）单击 MATLAB 命令窗口工具栏上的图标（ C ）从 MATLAB 命令窗口种选择菜单 file “new M-file 。（ 2 ）打开已有 M 文件的操作方法（ A ）直接在命令窗口运行命令edit filename filename 是待打开的文件名。（ B ）单击 MATLAB 命令窗工具栏上的图标，从弹出对话框中选取需打开的文件（ C ）单击命令窗口中的 file open 子菜单，从弹出对话框中选取需打开的文件编辑或修改的文件经保存后可单击 debug 项中的 run 进行运行或者在命令窗中直接输入此 M 文件的名称，就可观察运行结果如果有错误，会在命令窗口中给出相应的提示，根据出错提示找到错误进行修改3 、命令文件最简单的 M 文件，不含输入变量和输出变量，主要用途是使命令输入更简单，调用时 MATLAB 自动执行文件中的全部语句，因此在实现分析问题,解决问题和设计复杂命令等方面十分有用例1.8:建立一个命令文件，输出 500 以内的 Fibonacci 数打开文本编辑器，输入如下命令行，保存文件名为 fibon ，然后再命令窗口中输入命令 fibon,就能立即看到结果 fibonf = Columns 1 through 12 1 1 2 3 5 8 13 21 34 55 89 144 Columns 13 through 14 233 377图1-2 M文件4 、函数文件：可以传递参数，接受输入变量和返回输出变量且对变量的运算是在函数本身的工作空间内进行，和 MATLAB 命令行窗口的工作空间互相独立。函数文件的第一行必须包含关键字 “function函数名的命名规则和变量名的命名规则一样：必须以字母开头，其余部分可以是字母、数字、下划线例1.9: 建立一个函数文件，计算序列的平均值打开编辑器输入如下内容，并保存文件名为 aver.mfunction a=aver (x)% AVERAGE mean value% for vectors,AVERAGE is the mean value% for matrix, AVERAGE is a row vector% Non-vector input results in an errorm,n=size(x);if m=1; m=n;enda=sum(x)/m;然后在命令窗口中输入变量 x。x=26:1000;aver (x)ans = 513函数文件的几点说明：（ 1 ）文件的第 1 行说明了该函数的名称、输入参数、输出参数。（ 2 ）“%”后面的语句用来注释说明，只是为了程序易于理解（ 3 ）文件前面几行由 % 开始的语句构成了 M 文件的帮助信息（ 4 ）变量 m,n 及 y 为局部变量，程序执行完后，变量随之消失（ 5 ）函数体包含所有计算和为输出参数附值的代码（ 6 ）调用函数文件时，传递的参数名称可与函数文件中定义的参数名称不相同六、二维图形的绘制一、一般步骤（ 1 ）数据的准备 A、选定要表现的范围 B、产生自变量的采样向量 C、计算相应的函数值向量（ 2 ）选定图形窗及其子图的位置 典型命令：figure (1) %指定 1 号图形窗 subplot（ 2 ， 2 ， 2 ） %指定 2 号子图（ 3 ）调用绘图命令 典型命令：plot (x, y,-ro) %用红色实线画曲线，数据点类型为 o（ 4 ）设置轴的范围与刻度、坐标分格线 典型命令：axis (0,inf,-1,1) %设置坐标轴的范围 grid on %画坐标分格线（ 5 ）图形注释 典型命令：title(专家系统） %图名 xlabel( x);ylabel( y) %轴名二、plot函数的调用格式二维曲线的绘图命令中，函数plot是最基本,最重要的命令调用格式 1：poly(x,y) 绘制二元组 x,y 的曲线图形调用格式 2：plot (y) 绘制曲线图形 若 y 是实数向量，则绘制 y 与其元素下标所构成的二元组的曲线图 若 y 的元素为复数，则等价于plot (real(y),imag(y) 若 y 是矩阵，则按列绘制曲线图，曲线条数等于 y 矩阵的列数调用格式 3: plot (x1 ， y1 ， x2 ， y2,) 功能：以（ x1,y1)， (x2,y2) 为二元组，绘制多条曲线例1.10：绘制 y=sin (x) ， z=cos(x),w=0.25x-0.5 的曲线，x=0:pi/100:2* pi;y= sin(x)plot (x,y) %绘制 y= sinx 的图形图1-3 绘制 y= sinx 的图形在命令窗口中继续输入以下命令z= cos(x)w=0.25* x-0.5;plot(x , y , x , z , x, w)图1-4 绘制x,y;x,z;x,w的图形三、基本的绘图控制1 、色彩、线型和数据点标志函数 plot 可以接受字符串输入变量，以指定不同线型,数据点标志符号以及每条曲线的颜色调用格式：plot (x ,y,s)其中 s 由 13 个字母组成，用来制定所绘制曲线的色彩,线型和数据点标志。表1-1 PLOT函数参数说明符号代表的线型符号表示的颜色实线（缺省值）c(cyan)青色-虚线m(magenta)品红色：点连线y(yellow)黄色-.点划线r(red)红色none不画线g(green)绿色b(blue)蓝色w(white)白色k(black)黑色标记符号说明+加号标记O圆圈标记*星号标记.点标记X叉型标记S正方形标记（square）D菱形标记（Diamond）朝上的三角形符号V朝下的三角形符号朝右的三角形符号 x=0:pi/100:2*pi; y=sin(x); plot(x,y,-.mo)图1-5 默认的坐标轴刻度生成图形键入命令修改坐标轴的范围 axis(-1,11,-2,2);图1-6 修改坐标轴后的图形键入命令修改坐标轴的刻度内容 set(gca,xtick,-1,3,7,11) %在x轴上的-1,3,7,11处标记刻度图1-7 修改坐标轴刻度内容后的图形3、坐标网络使用”grid”命令可在当前作图中附加网格线调用格式1：grid on 功能：在当前坐标中附加网格线调用格式2：grid off功能：移掉网格线调用格式3：grid功能：切换grid的on与off两种状态例1.12：给上图加上网格线键入命令：grid on 或 grid图1-8 加上网格线后的图形4、图形的标注和图例典型命令有：title xlabel ylabel text gtext调用格式1：title(text)功能：给图形加上标题调用格式2：xlabel(text)功能：对x轴加上标题调用格式3：ylabel(text)功能：对y轴加入标题调用格式4：zlabel(text)功能：对z轴加入标题调用格式5：text(X,Y,string)功能：在指定位置（X,Y）加入一个文本字符串调用格式6：gtext(string)功能：利用鼠标在图形中加入文本字符串调用格式7：legend(sting1，string2,)功能：给图形加入图例，对当前图形建立一个图例说明盒，盒内给出每条不同类型的曲线所代表的不同函数名称下面是演示图形控制的例子： x=linspace(0,2*pi,50); y=sin(x); z=cos(x); w=0.25*x-0.5; plot(x,y,b:,x,z,r-,x,w,-.b*); % 绘图 xlabel(X axis); % x轴的标注 ylabel(function y,z and w); % y轴的标注 title(Three Functions); % 图形的标题 grid on % 添加坐标网格 text(6,-0.5,y=sinx); % 在指定位置添加文本 text(2,-0.5,z=cosx); gtext(w=0.25x-0.5) % 用鼠标选定位置添加文本 legend(y=sinx,z=cosx,w=0.25x-0.5) %制作图例说明盒图1-9 添加图例说明盒五、图形的叠加功能”hold” 命令功能：在已有的图形系统中附加新的图形成分调用格式1：hold on功能：保留当前图形及坐标的全部属性，随后绘制的图形附加到现有图形上调用格式2：hold off功能：返回hold的缺省模式，随后的作图命令“ plot ”将抹去当前已有图形，在绘制新的图形前重新设置坐标轴的属性调用格式3：hold功能：切换hold的on和off两种状态例1.13：演示图形的叠加功能 x=0:pi/50:4*pi; y=sin(x); plot(x,y,r)在做好的图基础上再添加z=sin(x-0.5)和s=sin(x-1)的图形，可继续输入以下命令 hold on z=sin(x-0.5); s=sin(x-1); plot(x,z,bo,x,s,-.m*) hold off 图1-11添加函数z,s后的图形图1-10 y=sin(x) 的图形六、子图利用”subplot”函数，可在同一图形窗口显示多个图形调用格式；subplot(m,n,p)功能：将当前图形窗口分割为mn个绘图子区域，并将第p个区域作为当前的绘图区图形区的号码编排采取行优先原则：所有子图沿第一行从左往右编号，然后是第二行，以此类推例1.14：把一个图形窗口分为四个子图，每个子图绘制不同图形 x=-3*pi:pi/50:3*pi; x11=x.*cos(x); x12=-x.*sin(x)+cos(x); x21=sin(x).*cos(x); x22=(sin(x).2).*(cos(x).2); subplot(2,2,1) plot(x,x11) title(第一个子图) 图1-12 子图的绘制七、枝干图的绘制枝干图将每个离散的数据显示成尾部带有标记符号的线条，称为枝干“stem”命令用于绘制二维枝干图，其线条起点在X轴上例1.15：绘制函数y=e-atcos(bt)的图形命令如下： a=0.02;b=pi/15; t=0:4:120; y=exp(-a*t).*cos(b*t); subplot(2,1,1); plot(t,y); % 绘制用直线连接离散点的二维图形 xlabel(t axis) ylabel(y axis) subplot(2,1,2); stem(t,y); % 二维枝干图 xlabel(t axis) ylabel(y axis)图1-13枝干图的绘制三. 实验报告要求(1) 简述目的以及实验内容(2) 按实验步骤附上实验过程中各例子的执行结果，并与指导书给出的结果进行对比分析(3) 简要回答思考题四. 思考题(1) 例1.14中，程序没有全部给出，余下三个图形的绘制是如何实现的？如果改成3X3的图形又当如何实现？实验二 用FFT做谱分析实验项目名称：用FFT做谱分析实验项目性质：综合性实验所属课程名称：数字信号处理实验计划学时：4一. 实验目的(1) 进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。(2) 熟悉FFT算法原理和FFT子程序的应用。(3) 学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。二. 实验内容和要求(1) 对下列典型信号进行谱分析：，这里给出针对各信号的FFT变换区间N以及对连续信号的采样频率，供实验时参考：(2) 令，用FFT计算8点和16点离散傅立叶变换，并根据DFT的对称性，由求出和（1）中所得结果进行比较。提示：取时，。(3) 令，重复（2）。三. 实验主要仪器设备和材料计算机，MATLAB6.5或以上版本四. 实验方法、步骤及结果测试(1) 复习DFT的定义、性质和用DFT做谱分析的有关内容。(2) 复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图，读懂本实验提供的FFT子程序。(3) 编制信号产生子程序，产生下列典型信号供谱分析用。应当注意，如果给出的是连续信号，则首先要根据其最高频率确定采样速率以及由频率分辨率选择采样点数，然后对其进行软件采样（即计算，），产生对应序列。对信号，频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行谱分析，否则有可能产生较大的分析误差。请实验者根据DFT的隐含周期性思考这个问题。(4) 编写主程序。图2-1给出了主程序框图，供参考。开始读入长度N调用信号产生子程序产生实验信号调用绘图子程序（函数）绘制时间序列波形图调用FFT子程序（函数）计算信号的DFT调用绘图子程序（函数）绘制曲线结束图2-1 主程序框图本实验的主程序比较简单，直接根据上图给出的框图编写主程序即可，编程中的难点是FFT子程序。不过，各种语言的FFT子程序都可以在有关的信号处理的程序库中找到，由于C语言当前最普及，所以为了为实验提供方便，下面给出C语言FFT函数，供参考/* DIT-FFT函数（C语言）*/fft基2DITFFT函数要求：指向复数数组指针X，FFT长度为2m，m为正整数FFT输出结果放在输入复数数组中。/* 计算N点FFT子程序*/* xr：=信号序列实部，xi：=信号序列虚部，N：=FFT变换区间长度N=2M */* 如果信号长度小于N，应该给xr，xi后面补0 */* 计算如果X（K）的实部和虚部分别储存在数组xr和xi中 */Void Fft(double xr, double xi, int N, int M)int L, B, J, P, k, i;double rPartKB, iPartKB;double rCf128, iCf128/* 计算旋转因子 */double PI2= 8.0 * atan(1.0);for(i=0; iN; i+)rCfi=cos(i*PI2/N);iCfi=sin(i*PI2/N);ChangeOrder( xr, xi, N );/* 计算各级蝶形 */for(L=1; L=M; L+)B=(int)(pow(2, (L-1)+0.5);for(J=0; J=B-1; J+)P=J*(int)(pow(2,(M-L)+0.5);for(k=J; k=N-1; k+=(int)(pow(2,L)+0.5)rPartKB = xrk+B*rCfP-xik+B*iCfP;iPartKB = xik+B*rCfP+xrk+B*iCfPxrk+B = xrk- rPartKB;xik+B = xik- iPartKB;xrk = xrk + rPartKB;xik = xik + iPartKB;/* 倒序子程序*/void ChangeOrdor(double xr, double xi, int N )int LH, N1, I, J, K;double T;LH = N/2; J = LH; N1 = N 2;for(I=1;I=N1;I+)if(I=K)J = J-K;K = (int)(K/2+0.5);J=J+K;(5) 按实验内容要求，上机实验并写出实验报告。本实验采用的是MATLAB语言，因此FFT子程序直接调用MATLAB语言中的FFT函数就可以实现。下面给出完整的MATLAB程序%实验二，用FFT做谱分析b=menu(请选择信号x1(n)-x8(n),x1(n),x2(n),x3(n),x4(n),x5(n),x6(n),x7=x4+x5,x8=x4+jx5,Exit);if b=9 b=0;endi=0;close all;while(b)if b=6temp=menu(请选择FFT变换区间长度N,N=16,N=32,N=64);if temp=1 N=16; elseif temp=2 N=32; else N=64; endfs=64;n=0:N-1;x=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);elsetemp=menu(请选择FFT变换区间长度N,N=8,N=16,N=32); if temp=1 N=8; elseif temp=2 N=16; else N=32; endif b=1x=1 1 1 1 0 0 0 0;else if b=2 x=1 2 3 4 4 3 2 1; else if b=3x=4 3 2 1 1 2 3 4;else if b=4 n=0:N-1; x=cos(0.25*pi*n); else if b=5 n=0:N-1;x=sin(pi*n)/8);else if b=7 n=0:N-1; x=cos(n*pi/4)+sin(n*pi/8); else if b=8n=0:N-1;x=cos(n*pi/4)+j*sin(n*pi/8);end endend end end end end end%TO Calculate FFTf=fft(x,N);i=i+1;figure(i);printf(x,abs(f),abs(N),abs(b);if N=16 if b=7 k=conj(f); x4=(f+k)/2;%ReX7(k)=x4(k) figure(i+2); subplot(2,2,1); stem(abs(x4),.); xlabel(k); ylabel(|X4(k)|); title(恢复后的X4(k); x5=(f-k)/2;%jImX7(k)=X5(k) subplot(2,2,3); Stem(abs(x5),.); xlabel(k); ylabel(|X5(k)|); title(恢复后的X5(k); end if b=8 k(1)=conj(f(1); for m=2:N k(m)=conj(f(N-m+2); end fe=(x+k)/2;%求X8(k)的共轭对称分量 fo=(x-k)/2;%求X8(k)的共轭反对称分量 xr=ifft(fe,N);%xr=x4(n) b=4; figure(i+1) printf(xr,abs(fe),abs(N),abs(b); xi=ifft(fo,N)/j;%xi=x5(n) b=5; figure(i+2) printf(xi,abs(f),abs(N),abs(b); endendb=menu(请选择信号x1(n)-x8(n),x1(n),x2(n),x3(n),x4(n),x5(n),x6(n),x7=x4+x5,x8=x4+jx5,Exit);if b=9 b=0;endclose all;end程序运行结果：按照实验内容选择x1(n)x6(n)以及x7(n)=x4(n)+x5(n)、x8(n)=x4(n)+jx5(n)进行谱分析，输出x1(n)x5(n)的波形以及8点DFT和16点DFT如图2-2到图2-11所示，输出x6(n)的16点、32点和64点采样序列及其DFT如图2-12、13、14所示。图2-2 x1(n)的8点DFT图2-3 x1(n)和16点DFT图2-4 x2(n)和8点DFT图2-5 x2(n)和16点DFT图2-6 x3(n)和8点DFT图2-7 x3(n)和16点DFT图2-8 x4(n)和8点DFT图2-9 x4(n)和16点DFT图2-10 x5(n)和8点DFT图2-11 x5(n)和16点DFT图2-12 x6(n)和16点DFT图2-13 x6(n)和32点DFT图2-14 x6(n)和64点DFT选7时，首先计算并图示x7(n)=x4(n)+x5(n)*R8(n)和x7(n)=x4(n)+x5(n)*R16(n)及其DFT，如图2-15、2-16所示；然后程序自动计算并绘图验证DFT的共轭对称性。当N=16时，x4(n)=x4(N-n)，x5(n)=-x5(N-n)。即x4(n)为x7(n)的共轭对称分量，而x5(n)是x7(n)的共轭反对称分量。根据 DFT的共轭对称性，应有如下结果：X7(k)=DFTx7(n)16点=ReX7(k)+j ImX7(k)X4(k)=DFTx4(n)16点=ReX7(k)X5(K)=DFTx5(n)16点=j ImX7(k)图2