matlab_ppt教程_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
评价文档：
喜欢此文档的还喜欢
matlab_ppt教程
讲​授​M​A​T​L​A​B​语​言​基​础​入​门​知​识​,​介​绍​M​A​T​L​A​B​产​品​的​体​系​、​M​A​T​L​A​B​桌​面​工​具​的​使​用​方​法​,​重​点​介​绍​M​A​T​L​A​B​的​数​据​可​视​化​、​数​值​计​算​的​基​本​步​骤​以​及​如​何​使​用​M​A​T​L​A​B​语​言​编​写​整​洁​、​高​效​、​规​范​的​程​序​。​并​涉​及​到​一​些​具​体​的​专​业​应​用​工​具​箱​(​如​：​信​号​处​理​工​具​箱​、​图​像​处​理​工​具​箱​等​)​。​
​通​过​本​课​程​的​学​习​,​了​解​、​熟​悉​、​掌​握​ ​M​A​T​L​A​B​的​基​本​编​程​方​法​,​并​具​有​初​步​的​利​用​计​算​机​处​理​、​解​决​实​际​问​题​的​能​力​,​为​进​一​步​学​习​后​续​的​专​业​课​程​做​好​准​备​。
阅读已结束，如果下载本文需要使用
想免费下载本文？
把文档贴到Blog、BBS或个人站等：
普通尺寸(450*500pix)
较大尺寸(630*500pix)
你可能喜欢第九章综合应用1综合应用本章内容要点： ? 9.1 信号处理工具（sptool）的介绍 ? 9.2 语音数字化量化噪声的改善 ? 9.3 系数量化和运算量化的影响 ? 9.4 在双音拨号系统中的应用 ? 9.5 正余弦信号的谱分析 ? 9.6 音乐信号处理 ? 9.7 变采样率数字滤波 ? 9.8 稀疏天线阵列设计 ?
9.9 结束语2综合应用信号处理的概念比较深，其中许多概念不是单纯 用数学推导就能掌握的，往往要通过在实践中 反覆运用，才能真正掌握。开辟这一章的目的 就是通过一些实例来说明如何把理论用到工程 实践中去。同时又加深了对理论的理解。 另外在本章中还要介绍一些MATLAB的工具，本 章先介绍信号处理工具sptool，它是信号处理 工具箱中的一个集成环境。学习和使用这个工 具可以帮助读者把已学的信号处理知识系统化 和集成化，对于在工程中进行信号处理会有很 大的帮助。3综合应用9.2节介绍信号的量化问题，将介绍用压缩扩张 器减小相对误差的非线性处理方法；9.3节介绍 滤波器系数量化和计算量化的影响；9.4节介绍 双频拨号系统，它是一个从双频发送到接收检 测的完整的系统，其中包括goertzel算法开发； 9.5节介绍了频谱分析仪中的数字信号处理技术； 9.6节介绍音响系统中的信号处理，它涉及回声 和混响的生成和处理；9.7节介绍变采样率系统， 包括内插和抽取，这也属于一种非线性处理； 9.8节把时域信号扩展到空域信号处理；4综合应用上述的七节都是以前面学过的八章为基础，但都 把原有的知识拓宽和加深了一步。在本书中， 我们对每一个问题都只作了简要的提示，各个 问题都有很大的展开余地，所以其中任何一个 问题都可以作为实验或课程设计的基础。这要 由教师和读者自行取舍，进行适当的开发。 9.9节是结束语，大略地介绍了一下数字信号处理 器的应用领域，并介绍了MATLAB中与信号处 理有关的各种工具箱函数。59.1信号处理工具介绍(讲本节时，应尽量直接用MATLAB的图形界面) ? 信号处理工具sptool（Signal Processing Tool的 缩写）为信号处理的研究工作提供了一个集成 环境和工具。 ? 信号处理的任务：一是对信号进行分析；二是 滤波器设计。在这两个任务进行过程中，经常 要做第三个任务，那就是要把信号加到所设计 的滤波器中，看它的输出是否满足要求。 ? Sptool把这三个任务集成在一起，加上了适 当的管理功能，配以良好的工作界面，利用它 可以大大提高研究工作的效率。6信号处理工具(sptool)介绍在MATLAB命 令窗中，键 入fdatool, 得到如右图 的界面。各 栏中分别存 入了系统中 原来已保存 的信号、滤 波器和频谱 的名称。图 9.1.1 sptool 的启动界面7信号处理工具(sptool)介绍? 1． 信号和滤波器的导入 先在MATLAB工作空间中放入所需的分析的信号 和滤波器。设滤波器就是例8.8.1中设计并导 入了工作空间的滤波器，再来建立一组信号。 在命令窗中，键入： n=0:200;T=1;s1=sin(5*2*pi*n*T); s2=sin(10*2*pi*n*T); s3= cos(15*2*pi*n*T);s=s1+s2+s3; ? 这时s1,s2,s3和它们的合成信号s都已输入 工作空间。8信号处理工具(sptool)介绍点击左上角的【File】及下拉菜单的【Import】, 出现如下图的新视窗。它分成三栏； 左边一栏【Source】，指定数据来源；如果来源 是工作空间，中间栏的【Workspace Contents】 将显示工作空间中的全部变量；右边一栏指定 导入的目标，即要说明导入的数据是作为信号、 还是滤波器还是频谱。右栏的下部有 【Sampling Frequency】和【Name】两个框需 要用户填写。按照说明的方法，可依次将信号 s和8.6节设计的滤波器输入sptool。9信号处理工具(sptool)介绍10图 9.1.2sptool 的 数 据 导 入 界 面信号处理工具(sptool)介绍2. 信号的时域和频域观测 在 图 9.1.1 上 ， 选 定 sig7, 点 击 【view】( 观 测),就进入【Signal Browser】(信号浏览 器)视窗。如图9.1.3所示。要想求得该信号的频谱，也先选定sig7，然后在 频谱栏下点击【creat】,此时出现频谱观测器 的界面，见图9.1.4。先要在左栏上部选定求 频谱的方法，目前我们只学了一种，即选FFT. 在下面的框中填入点数，例如1024，然后点击 【Apply】,就出现了频谱的曲线。11信号处理工具(sptool)介绍图 9.1.3信号浏览器（Signal Browser）的界面12信号处理工具(sptool)介绍图 9.1.4频谱观测器（Spectrum Browser）的界面13信号处理工具(sptool)介绍? 3． 滤波器的观测和修改 选 定 滤 波 器 的 名 称 filt7, 点 击 其 下 方 的 【View】，就可以观测它的频率响应。如点 击【Edit】，那就不仅能观测，还能修改。 ? 4. 让信号通过滤波器求输出信号 在图9.1.1中，选定信号和滤波器名称，点击滤 波器栏下的按钮【Apply】,这意味着将所选信 号加到所选滤波器中去。此时将产生一个小视 窗，提示用户输入信号和滤波器的名称，并提 供输出信号缺省名称sig8。如无修改，点击OK， 输出信号将自动生成并列写在信号栏中。14信号处理工具(sptool)介绍在本例中滤波器输出sig8基本上是一个较 纯的10Hz正弦波，这是不难想象的。经 过仔细观测和分析，如不满意，可以修 改滤波器后重新实验。 用这个工具，不难在很短的时间内完成多 个滤波器的设计，并分别观测它们的输 出。所有的结果都存储在案，可用于试 验报告。所以这是一个提高效率的有效 工具。159.2语音量化噪声的改善脉冲编码调制(Pulse Code Modulation―PCM)是把 模拟信号量化为二进制数的最简单的方法。以 N个脉冲表示N位二进制数，以脉冲的有无判 断它是0或1。它也是用数字方式传输或存贮信 号的常用方法之一。PCM被广泛应用于电话通 信和利用无线电传输的遥测系统中。 通过电话线传输的语音信号频带限于4kHz以下。 因此其采样频率取8KHz(样本数/秒)，并用N位 二进制序列表示它的值，每个样本量化为2N个 电平之一。所以，传输数字化语音信号所要求 的速率为每秒8000×N位。16语音量化噪声的改善量化处理的数学模型为 xq(n) =x(n)+q(n) (9.2.1) 其中xq(n)表示x(n)的量化值，q(n)表示量化误差， 将其看作一加型噪声。假设采用的是均匀的量 化器，则可用如下均匀概率密度函数p(q)统计 描述量化噪声特性: (9.2.2) 其中，量化步长为△=2CN。量化误差的均方值为: 2 ?2 N ? 2 (9.2.3) 2 E (q ) ? ? 12 12171 ? ? p(q) ? , ? ? q ? ? 2 2语音量化噪声的改善用分贝来度量的噪声均方值为: ? 2 ?2 N ? 10 log E ? 10 log ? ? ? ?6 N ? 10.8 dB(9.2.4) ? 12 ? 可以看出，上述的量化器每增加一位，量化噪声 减小6dB，高质量语音要求每个样本至少量化 为12位，因此传送速率至少为96000位/秒。 最大幅度为±V伏的N位（不含符号位）二进制 A/D变换器的数学模型建立如下。它把电压V分 解为2N-1份，故量化步长为V/(2N-1)，得出二 进制量化子程序bqtize。18语音量化的子程序function y=bqtize(x,N,V) if nargin&3 V=max(abs(x));end % V缺省取最大x ax=abs(x); % 去掉符号 deltax=V/(2^N-1); % 求量化步长 xint=fix(ax./deltax+0.5); % x的量化整数 y=sign(x).*xint.* % 恢复量化原值 这个A/D变换子程序的输入是连续模拟电压x，输 出则是量化了的模拟电压y。均匀量化器在信 号的整个动态范围中的量化步长相同，所以量 化噪声均方值不变。19语音量化噪声的改善均匀量化器在信号的整个动态范围中的量化步长 相同，所以量化噪声均方值不变。然而，语音 信号的特性是小幅度比大幅度出现得频繁。对 小信号而言，量化噪声使信噪比大大下降。解 决的途径之一是用非均匀量化器。不过在技术 上制造非均匀量化器的芯片是困难的。 得到非均匀量化器特性的另一个方法是用压缩― 扩张器。可先使信号通过压缩幅度的非线性器 件，后面再接一均匀量化器，再用逆向扩张幅 度的非线性器件恢复信号。如下图。20语音量化的压缩扩张器图 9.2.1 压缩器、A/D 变换器和扩张器的连接框图压缩器的作用是把小信号放大，大信号缩小，所 以把压缩器后的信号进行量化，小信号的信噪 比就得到提高。而后通过扩张器，信号和量化 噪声同时作非线性变换，信号复原为原来电平， 小信号信噪比则保持较小的水平。21语音量化的压缩扩张器在我国的通信系统中使用的对数压缩器(称为μ律 压缩器)具有如下输入输出幅度特性: V ln(1 ? ? | x / V |) y? sign( x), x ?V, y ?V, ln(1 ? ? )其中，x是归一化输入，y是归一化输出，sign(.) 是符号函数，μ =255,是控制压缩特性的参数。 由量化数字信号恢复模拟信号时，逆关系μ律为:x? V (1 ? ? )| y / V | ? 1?sign( y )22语音量化的压缩扩张器可以用MATLAB把压缩器和扩张器用函数程序 mulawcom和mulawexp来表示，例如把压缩器 函数程序列写如下： function y= mulawcom(x,mu,V) if nargin&3 | V & max(abs(x)) V=max(abs(x)); end y=V/log(1+mu)*log(1+mu/V*abs(x)).*sign(x); 对μ律扩张器和A-律的压扩器也可以按类似的方法 写出相应的子程序。23语音量化压缩器的特性例9.2.1 画出μ 律和A律的压缩器输入输出曲线。 解：写出如下的MATLAB程序hc921 x=0:0.01:1; y=mulawcom(x,255,1); y1=Alawcom(x,87.56,1); plot(x,y,x,y1,':');grid on legend('\mu律','A律') 运行此程序的结果见图9.2.1。请注意二者非常相 似。24语音量化压缩器的特性图 9.2.2μ -律与 A-律非线性压缩器的比较25语音量化噪声的改善例9.2.2 设某A/D变换器把最大输入为5伏的 信号量化为四位二进制（不含符号位），要求 用图形描写其输入输出关系，并画出其绝对误 差和相对误差的曲线。 又若信号像图9.2.1那样经过压缩扩张器，则输 入输出关系有何变化？ 解： A/D变换器的数学模型为绝对量化函数 bqtize，已经在前面得出。利用这个函数，加 上求绝对误差和相对误差的语句，可以方便地 列出以下的MATLAB程序hc922。26语音量化的信噪特性对比x=-5:0.01:5; % 输入自变量数组 xq=bqtize(x,4,5); % 求量化输出 e=x- er=e./abs(x); % 求相对误差er plotyy(x,xq,x,er) % 画输出及误差曲线 plot(x,e,'-.','linewidth',3) 运行此程序所得曲线见图9.2.3左图。 如在第二条语句的前后分别加上压扩语句： x1=mulawcom(x,255,V); % 信号经过压缩器 yq=mulawexp(xq,255,V); % 信号经过扩张器 则所得曲线见图9.2.3右图。27语音量化的信噪特性对比图 9.2.3二进制 A/D 不用压扩器（左）和用压扩器（右）变换输入输出及误差曲线28语音量化的信噪特性对比从左图中可以看出，在输入为-5~5V范围内， 绝对量化误差e呈等幅锯齿波形式，其最 大值恒定，因而其相对误差er在小的输入 幅度x处急遽增大，超过了1。 由右图从中可以看出它的绝对误差随x的增 加而增大，其相对误差则在整个输入范 围内呈等幅锯齿波形式，最大不到0.2。 因此量化误差yq-x1引起的信噪比将比不 用压扩器有显著的提高。299.3 系数量化和运算量化? 9.3.1 数的浮点和定点表示方法 ? 要用计算机处理信号，就不可避免地要遇到量 化问题。这里包括三方面：信号的量化、系统 参数的量化以及信号与系统相互作用，也就是 运算中的量化。关于信号的量化本书已经讨论 过，本节将讨论其他两方面的量化问题。我们 先看一下MATLAB中的量化，作为一种通用计 算机上的算法语言，它用很多位数和浮点方法 来表示数和数的运算。所以它的量化步长很小， 精度很高，完全可以当成模拟量来看待。30双精度浮点表示? 按照IEEE标准，双精度浮点数η用下式表示 ? ? (?1) S ? M ? 2 E ? （9.3.1） ? 其中M是一个小于1大于1/2的二进制分数，称 为尾数，用52个二进制位来表示； ? 指数E是一个带符号的二进制整数，也称为阶 码 。 占 11 个 二 进 制 位 ， 可 以 表 示 从 ? 1023 到 ?1024的数集，它决定了数的动态范围。数的 正负号反映在S上，它只占一位。所以双精度 浮点数总计占1?52?11?64个二进位。31浮点数的相对精度? 浮点数的量化步长可以代表它的相对精度。它 是由M的位数决定的。52位二进制数的量化步 长是2 C52?2.。该数的动态范围则取 决于指数部分。因为2 C7及2 ?。所以MATLAB中数的动态范围为 2.8～1.。在MATLAB命 令窗中，MATLAB中的大部分运算结果的误差 比eps 要大一些。在正常情况下（即系统不给 出警告时），计算结果至少可以有12位十进制 有效数字的可信度。32实际芯片中数的精度? 当要把用MATLAB设计的滤波器付诸实现，特 别是用于实时系统时，要采用嵌入式的数字信 号处理芯片，通常它必须采用较短的字长（8 位、16位或32位），用定点的方法（运算速度 快）来表示数并进行数的运算。因此在工程实 践中，绝不可以盲目轻信MATLAB中计算出的 结果，而必须考虑实际系统的有效字长对结果 进行复核和修正。当然实际系统的位数愈高， 则它将愈好地接近MATLAB计算的结果。33定点数表示法? 如果二进制数码的小数点在整个运算过程固定 不变，即上述浮点数中的E在整个计算系统中 取为常数，则称为定点制。实际上，为了运算 方便，定点制通常把数限制在?1与1之间。这 样就把小数点定位于第一位数码之前，即左边 第一位为符号位，代表数的正负号（0表示正， 1表示负）。数只需用符号位S和尾数M表示。 比如一个8位定点数，若取E?0，其尾数占7位， 则其量化步长为2 ?7，相对精度为0.0078。若取 E?1，其尾数占6位，则该定点数取值在?2与2 之间34定点数的定标? 定点制在整个运算过程中，若取E?0，则所有 运算结果的绝对值都不能超过1。在实际应用 中，当数很大时，应乘以一个比例因子，使整 个运算过程中的数的绝对值都不超过1，运算 完毕后再除以同一比例因子，还原为真值输出。 只有对所有的数和运算都取同样的阶码E才能 这样处理。阶码的确定对于定点系统的设计是 十分重要的任务。 ? 定点数的定标：定标就是把数变为定点表示法 时阶码或小数点位置的确定。35阶码或小数点位置的确定? 如果一个数有整数部分，那就只要把这个数的绝对值以2 为底的对数求出，然后向上取整，即得出其二进制整数位， 例如x?6.54，则log2(x)? 2.7093，E?3。实际上，这个公式 也 可 推 广 到 x 绝 对 值 小 于 一 的 情 况 。 例 如 x?0.035 ， 则 log2(x)? ?4.8365，E? ?4。阶码为负说明此定点数小于1/2， 因为按照（9.3.2），尾数M必须大于1/2。 ? 滤波器系数x通常是一个数组或矩阵，要把x进行量化，可 用9.2.1节中对信号的二进制量化子程序bqtize（x,N,V）， 该子程序的输入输出都用十进制数表示，比较容易读懂。 该子程序的输入变元中，N为二进制尾数的位数，V是量 化后的定点数的最大值，它应该比矩阵x中任何一个元素 的绝对值都大的2的乘幂。V =2^ceil(log2(max(max(abs(x)))))36多个数组用同一阶码情况当系统中有多个滤波器时，滤波器系数x可能包括多个数 组x1,x2,…, 其中有些系数数组（如sos）还是二维的， 量化必须对x中所有的系数按同一阶码进行，且V绝对 值必须大于其中任何一个系数的绝对值，以保证量化 不出现饱和错误。故应取 ? V1 =2^ceil(log2(max(max(abs(x1))))); ? V2 =2^ceil(log2(max(max(abs(x2))))); … ? V = max(V1,V2,…) ? 使一个系统中取一个共同的V，然后再调用子程序 bqtize（x,N,V）。37系数量化对FIR滤波器的影响? 例9.3.1设一个四阶FIR滤波器，其传递函数为：? （9.3.2） ? 算出它的级联结构，比较在两种结构下系数量化对频 率响应的影响。 ? 解：先求出其级联结构参数：0.360530z ?2 ? 0.388726z ?3 ? 0.069735z ?4? ? ? ? b=[0.....069735]; [sos,G]=tf2sos(b,1); 得 sos = 1.5 1.0 0 0 1.8 1.0 0 0 ? G = 0.0697H ( z ) ? 0.069735 ? 0.388726z ?1 ?? 如果把系数b量化为七位二进制数，称为bq，它的值可 由键入下面语句38系数组量化结果? V1=2^ceil(log2(max(abs(b)))), ? bq=bqtize(b,7,V1) ? 得到V1?0.5 ? bq = 0.6 % 求最大量化值0.35940.39060.0703? 如果滤波器由级联结构实现，则对它的参数进行量化? V2=2^ceil(log2(max(max(max(abs(sos))),G))), sosq=bqtize(sos,7,V2) ,Gq=bqtize(G,7,V2)? 得到V2?8 ? sosq = 1.5 1.0 1.0 1.0 ? Gq = 0. 0 039三种情况的幅频特性的比较? 第一种用原来的浮点双精度系数，第二种用对直接结 构量化后的系数，第三种用对级联结构量化后的系数，40系数量化对IIR滤波器的影响? IIR滤波器存在着递归计算的问题，它的 系数误差可能在循环计算中不断扩散， 因此，一般地说，对系数误差的敏感程 度比FIR滤波器严重。特别是当IIR滤波 器的极点靠近z平面上单位圆时，幅频特 性通常要出现很大的峰值。系数的量化 误差造成极点的移动，也很难用解析公 式表示，只能通过实例来计算和观察。41例9.3.2? 已知一IIR滤波器具有如下传递函数：（9.3.3）将其系数经过七位二进制舍入量化，分析其极点的变化。解：此系统在系数为无限精度时的极点可由下列MATLAB程序hc932求得：b=[1,0.4,-0.03,0.232]; % 分子系数向量 a=[1, -2.1,-1.0]; % 分母系数向量 V=2^ceil(log2(max(abs([a,b])))) % 求两系数向量最大量化值 aq=bqtize(a,7,V), bq=bqtize(b,7,V); % 量化后的系数向量 ra=(aq-a)./abs(a), % 量化后的分母系数的相对误差 p=roots(a),pq=roots(aq), % 分母系数量化前、后的极点 z=roots(b);zq=roots(bq); % 分子系数量化前、后的零点42量化引起的极点位移? 程序运行的结果如下：? ? ? ? V=4 aq = bq = ra = 1.8 2.8 0.0 0.3 0..0 -0.4? 这说明量化后的定点系数范围是在?4之间。从ra可看 出分子分母系数的相对误差并不大，下面看看极点位 置的变化。量化前后的四个极点分别为：? ? p?? ? ? ? ? 0.4930 + 0.8643i? ? ? 0.4930 - 0.8643i? ? , pq ? ? ? ? 0.9900 ? ? 0.4999 ? ? ? ? 0.5000 + 0.8660i? 0.5000 - 0.8660i? ? ? 0.9326 ? 0.5361 ? ?43滤波器结构与量化影响? 注意极点的模abs(p)的前两项为0.9950，而abs(pq)的前 两项为1.0000。说明量化使极点从单位园内移到了单位 园上，其后果是使本来稳定的滤波器成为不稳定的系 统。 ? 如果我们以级联方式实现这个系统，对它的级联系数 进行量化。程序为：? [sos,G]=tf2sos(b,a); % 化为二阶级联形式 ? V2=2^ceil(log2(max(max(max(abs(sos))),G))) sosq=bqtize(sos,7,V2), Gq=bqtize(G,7,V2) % 对各系数量化 ? a1q=sosq(1,[4:6]),a2q=sosq(2,[4:6]) % 取出分母系数向量 ? p1q=roots(a1q), p2q=roots(a2q) % 求系数向量的量化后的根? 就会发现，四个极点量化后都在单位圆内部，因此不 会造成滤波器的不稳定性问题44量化引起的极限环问题? 从这个例子可以看到，滤波器的级联结构对系数量化 的敏感性较之直接结构要小，这里的例子还只是四阶。 从多项式理论可以知道，高次多项式系数的微小变化 往往会引起根的显著变动。因此，阶次高的滤波器不 宜采用直接结构。 ? 在某些条件下，一个稳定的滤波器还可能因为量化效 应而在无输入的情况下产生持续的微幅振荡。这是一 种非线性振荡，称为‘极限环’。它的振荡幅度与量 化步长相近。现在用在信号处理中的单片机，位数也 都不小于16位。所以这类问题已经不很突出了，而且 研究这类问题还没有完美的理论，主要靠针对具体问 题的仿真计算来判别量化的影响459.3.4 运算量化和溢出问题? 信号处理的细致过程可以用信号流图来表示。其中主 要的运算就是加法、乘法和时延。加法器和乘法器如 果工作正常，应该能够保持运算的相对精度。由于定 点数的表示范围是一定的，因此在进行运算时，其结 果就有可能出现超过数值表示范围的情况，这种现象 称为溢出。如果结果大于定点数的最大值，称为上溢 （overflow）；如果结果小于定点数的最小值，称为下 溢（underflow）。一般情况下，上溢和下溢统称为溢 出。如不采取特别的措施，上溢将使符号位反号，定 点数突然由正最大值变为负最大值；下溢将使数变为 零。这些都可能意外地造成系统的崩溃。设计系统时 对溢出必须有保护措施，例如上溢时使加法器饱和而 不让它进（符号）位。因此，选择阶码时，不能只考 虑使它大于全部系数的绝对值，最重要的的一个指标 是使溢出的概率降到零或最低限度。 46溢出问题的研究方法? 两个定点数运算的的大体规则如下：相加精度基本不 变，但可能产生溢出；相减不会溢出，但精度降低， 所以要避免两个相近的大数相减；相乘不会上溢，但 字长增加一倍。所以乘法器中的累加器字长至少取16 位，不然会出现截尾误差或舍入误差；运算规则这样 无序，结果又依赖于参与运算的数。所以在运算过程 中在任何一个节点是否会出现溢出是很难靠理论分析 得知的。 ? 量化是一个非线性问题，没有封闭形式的数学解。而 且它和系统的硬件结构形式有关，即使简单地把两个 串联的环节前后颠倒一下结果就不同，不能用我们习 惯的理想数学等效模型去思考处理量化问题。一定要 使仿真用的信号流图以及系统的各个硬件结构真正符 合实际的系统模型，才能得到有意义的结论。47要靠定点仿真来解决? 在Simulink中早已有定点模块库（Fix-Point Blockset） 进行仿真，MATLAB 7.x新推出的定点工具箱（FixPoint Toolbox），提供了解决有限字长问题的工程途径。 Simulink新版本提供了更强的定点仿真的功能，在一个 滤波器系统的输入端加上给定的测试信号，然后记录 滤波器各个节点的信号值，求出其频谱并统计它们出 现溢出的次数和出现的概率。仿真中可以改变阶码， 提高阶码可使得其溢出概率下降，但阶码值取太高将 降低系统的精度，要综合考虑，使系统的特性尽量靠 近理论设计值。解决这样的工程实际问题必须更多地 依靠数学软件才行。48滤波器设计工具中的量化处理? 在8.8节的例8.81中，假如我们设计的是IIR切比雪夫II 型带通滤波器，其理论幅频特性已经显示在图形框中。 然后把左下方的量化按钮激活，界面就变成如图8.9.3 所示。其下半画面的左上角的【Filter arithmatic】是下 拉选项，将它选为Fixpoint，就出现了图9.3.2所示的参 数输入界面。 ? 故意把量化步长加大，以便明显看出定点设计和浮点 设计的差别。把其中的系数字长设为6位，先在其右方 小方框打勾，即由系统自动选择“最佳分数字长”。 再按正下方的【Apply】按钮，则幅频特性画面上将出 现两根曲线，其中虚线所示的是原来用双精度浮点算 出的理论曲线，而实线则表示用6位定点数对滤波器量 化后得到的幅频特性。4950Fdatool中作系数量化? 如果不在“最佳分数字长”的小方框打勾，则分子、 分母和增益（Scale）的小数字长（Fraction length）都 处于激活的可填充状态，用户可自行设定。若设为4位， 即留出一个符号位和一位阶码（也可用其下方的“取 值范围”为?2来定标），则得出的结果就是最佳的量 化定位，读者可试取其他的设定，得出的滤波器响应 都比图9.3.2差。除了系数页面之外，还可从下半平面 的右上方的标签切换到输入/输出和内部结构两个页面， 分别指定输入输出的量化参数和系统内部加法器和乘 法器的量化参数，在本例中取的是系统的默认值。其 输入字长是16位，输出字长是33位，乘法器和加法器 都能保证全精度，没有溢出。51设计结果? 依次按下菜单下方的各个按钮，可以得到滤波器的相 特性、脉冲响应、零极点分布等等图形。按下【系数】 钮，将可以看到原双精度设计所得系数及经过量化后 的 系 数 。 也 可 由 【file】→【Export】→【Coefficient （ASCII file）】将系数存成一个文件。本例所导出的 定点滤波器的三个双二阶环节的参数和增益分别为：?1 -1.875 1 1 -1. ? ?1 -1. -1.25 ? , G ? ?0.0625? sos ? 0.5 ? ? ? ? ? ?1 0 -1 1 -1. ? ? ?? 可以看出，所有的系数都处在?2～2范围之间，这和我 们规定的量化范围相符，也就是说，所取的定点数阶 码位数为1是合适的。529.3.6 定点（Fix-point）工具箱? 在MATLAB 7.0以上的版本中，Mathworks已经开发了 定点（Fix-point）工具箱，可以用它的函数fi来进行 量化，代替本书自编的bqtize函数。例如，在例9.3.2的 程序中，可以改用 ? aq=fi(a,1,8), bq=fi(b,1,8) ? 来进行量化。这里的fi的三个变元中，第一个是待量 化的变量名，第二个变元为‘1’是有符号位，为‘0’ 是无符号位，第三个变元为量化后二进制数的字长。 ? 与aq?bqtize(a,7,V), bq?bqtize(b,7,V)相比， fi中的8位二进制数去掉一个符号位相当于bqtize中 的7位，fi的优点是不需要另外用语句求量化最大值V。 当bqtize中的V值取得正确时，；两者的结果是相同 的。fi的功能很强，这里只给出了一种调用格式，还 有很多种输入变元格式不再列举。 53定点工具箱的其他功能? 定点工具箱中还给出了多种定点运算的函数，包括四 则运算、FFT、到卷积、滤波等信号处理的各种计算。 如果读者的MATLAB版本较高，可以用help命令查到 这些函数的调用方法。 ? 定点的MATLAB程序可以编译成为可执行的mex文件， 因而运行速度可以高达编译后的C代码速度，此外，定 点变量及其运算可以比较方便地与C语言中的定点程序 建立沟通的关系，可以在定点的C程序中调用 MATLAB定点程序来帮助它确定量化的适当参数，在 MATLAB的定点演示程序中提供了相应的例子，有兴 趣的读者可以参阅。因为大部分芯片是用定点格式来 表示数的，在把MATLAB的概念设计程序向数字信号 处理芯片代码实现的自动设计过程中，定点工具箱是 一个不可缺少的软件包。549.4 在双音拨号中的应用双音多频(DTMF―Dual Tone Multi Frequency)是 用按键进行电话拨号的体制。它不单用在电话 中，还可以用于传输十进制数据的其它通信系 统中。DTMF也广泛应用于电子邮件和银行系 统。 第一章1.2节中已经对它的功能进行了描述，4个 低频频率表示四行，3个高频频率表示三列， 两者的组合共可提供4×3＝12个字符。包括0 到9的十个数字，加上字符*和#。本节要讨论 如何解决双音变为数字的问题。55在双音拨号中的应用DTMF通信系统是一个很典型的小型信号处理系 统，它既有模拟信号的生成和传输部分（这要 用到D/A转换）；又有把它转为数字信号（这 要用到A/D转换）并进行数字处理的部分；而 为了提高系统的检测速度和降低成本，还开发 了一种特殊的DFT算法，称为Goertzel算法；这 种算法在国外的几乎每一本数字信号处理的教 材上都要介绍，说明了人们对它的重视。这种 算法既可以用硬件（专用芯片）、也可以用软 件实现，所以DTMF系统的设计问题是理论与 工程相结合的一个很好的典范。56在双音拨号中的应用先研究双音信号的生成问题： DTMF的两个音频可以用计算法或查表法 产生。用计算法得到正弦波形的缺点是 要占一些运算时间；查表法的速度较快， 缺点是要占一定的存贮空间； 两个正弦波的数字样本按比例相加在一起。 因为采样频率是8KHz，硬件必须每125 毫秒输出一个样本。将这个叠合信号送 到D/A变换器变换成模拟音频信号，通过 电话线路传送到交换机。57在双音拨号中的应用在接收端，将收到的模拟音频信号进行A/D变换， 恢复为数字信号，然后检测其中的音频频谱来 确定所发送的数字。检测算法可以用FFT，也 可用一组滤波器来提取所需频率。当要检测的 音频数目比较少时，用滤波器组更节省硬件。 如果用FFT算法实现该DFT计算，计算量(复数乘 法和加法)是N log 2 N。好处是立即得到DFT的 所有N个值，但至少要N个存储器。然而，如 果只希望计算DFT的K个点,而K&& N ，则直接 计算可以节省很多内存。下面介绍的Goertzel 算法就属于后者。58DFT的Goertzel算法Goertzel算法利用相位因子Wnk的周期性，将DFT 运算表示为线性滤波运算，由于WN - kN =1，可 用该因子去乘DFT，则? X (k ) ? ? x(m)WN k ( N ? m ) , k ? 0,1,? , N ? 1 m?0 N ?1它是两个序列的卷积。一是长度为N的输入序列 ? x(n)，另一则是脉冲响应序列 hk (n) ? W N kn ? (n), 可以看作x(n)通过具有该脉冲响应的滤波器。y k ( n) ?N ?1 m ?0?x(m)W N? k ( n ? m)? x ( n) ? W N? k ( n),该滤波器在n＝N点的输出就是频点k处的DFT.59DFT的Goertzel算法单位脉冲响应为hk(n)的滤波器的系统函数 就是(9.4.4)式的z变换，容易求出为1 H k (k ) ? ? 1 ? W N k z ?1这个滤波器只有一个位于单位圆上频率为 ωk=2πk/N处的极点。Goertzel算法的好处 不在于节省时间，而在于节省空间。如 果只需要K个DFT样本，可以只用K个并 行的单极点滤波器来分别计算这K个样 本，，这就可以大大节省硬软件资源。60DFT的Goertzel算法根据这个滤波器的差分方程，因为? y k (n) ? WN k y k (n ? 1) ? x(n),y(?1) ? 0可以用迭代方法计算yk(n)，如下图。要执行该计算， 可以只算一次相 位因子WN - k，将 其存贮起来。迭 代N次后，预期 的输出为： X(k)=yk(N)。图 9.4.1(a)计算单点 DFT 的递推框图61DFT的Goertzel算法上式中包含对 WN? k 的复数运算，想避免它，可 将成对的复共轭极点 WN? k 和WNk ? WN? ( N ?k ) 组 合 在 一起计算，导出双极点滤波器的系统函数:? 1 ? WN k z ?1 H k (k ) ? 1 ? 2 cos(2? k / N ) z ?1 ? z ?2它可以看成用两个差分方程组的联立: 2? kvk (n) ? 2cos N? yk (n) ? vk (n) ? WN k vk (n ? 1),vk (n ? 1) ? vk (n ? 2) ? x(n), (9.4.9)(9.4.10)初始条件为 v k(-1)=vk(-2)=0。运算结构图见后：62DFT的Goertzel算法图 9.4.1(b)计算 DFT 的双极点谐振器实现63DFT的Goertzel算法(9.4.9)中的递推关系对n=0，1，……，N重 复N+1次，每次只需要计算一次实数乘和 两次实数加，而带有复数运算的方程 (9.4.10)仅在n＝N时刻计算一次。所以， 对实数序列x(n)，由于对称性，用这种算 法求出X(k)和X(N-k)的值只需要N次实数 乘法和一次复数乘法运算。 实际上，不必计算复数值X(k)，只要求出 幅度平方值|X(k)|2 ，故计算的最后一步 还可以简化，从而完全避开复数运算。64Goertzel算法子程序按(9.4.7)式编成子程序gfft。它根据输入序列x， 和规定的DFT样本序号k，计算DFT样本X。 function X=gfft(x,k) % 用Goertzel算法计算序号为k的DFT样本 N=length(x);x1=[x,0]; % 递推要N+1次， d1 = 2*cos(2*pi*k/N); %滤波中间项系数v = filter(1,[1,-d1,1],x1);% 用滤波实现卷积W = exp(-i*2*pi*k/N); % 为下一步求W X = v(N+1) - W*v(N); % 求出第k个DFT65在双音拨号中的应用DTMF信号参数选择的考虑因素： (1).为了抗干扰，要检测8个频率及其倍频； (2).为了能分辨清这16个频率，要求采集数据的 长度N&Fs/D，已知Fs=8000，而根据最小相邻 基频间隔（73Hz）大于2D的要求，D&37，得 到N应在200以上。而为了提高检测速度，又希 望N尽量小些。 (3).为了检测尽量精确，希望这16个离散频点靠 近k=整数处，这需要仔细选择N。 最后选定的规范值是N=205。66在双音拨号中的应用本节用一个MATLAB程序hc941来演示DTMF双 频拨号的全过程。包括以下几部分： （1）首先程序要求用户输入一个数字，然后根 据这个数字查出它对应的两个频率，并生成相 应的双频信号x(n)，并产生声音作为标志。 （2）在接收端，程序截取收到的双频信号并进 行傅立叶变换，求出它在八个规定频率上的样 本幅度，画出频谱图。 （3）从频谱中分别取出满足规定电平的高低频 分量的两个下标。再由下标找到它的ASCII码 和数字。67程序hc941的说明此程序分为三段。 第一段是双频信号的生成，从输入一个数 字开始，以产生相应的双频信号结尾； 第二段是接收端信号处理，从收到双频信 号开始，进行截断和量化，并计算八个 基频DFT，以画出DFT幅度分布图结尾； 第三段是根据基频幅度分布图,截取其两个 峰值，查找和显示输入字符。 各段之间用pause命令隔开。68程序hc941的说明d=input(‘键入一位电话号码= ’,‘s’); symbol=abs(d); % 求它的ASCII码 tm=[49,50,51,65;52,53,54,66;55,56,57,67;42,48,35,68]; % 16个键的ASCII码，检测输入码与谁相符 for p=1:4; for q=1:4; if tm(p,q)==abs(d); break,end % 检测列号q end if tm(p,q)==abs(d); break,end % 检测行号p end69程序hc941的说明f1=[697,770,852,941]; % 行频率向量 f2=[77,1633]; % 列频率向量 n=0:2040; % 为了发声，加长序列 % 构成双频信号 x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000); sound(x); % 发出声音 disp('双频信号已经生成并发出'),pause ％ 以上是第一段程序，表示发送端的功能。70程序hc941的说明（第二段）接收检测端的计算8点DFT的程序 N=205; % 截取信号的长度 k = [18 20 22 24 31 34 38 42]; %待求的DFT序号 for m = 1:8; X(m)=gfft(x(1:N),k(m)); % 算八点DFT值 end val = abs(X); stem(k,val,'.'); %画出八点DFT向量 disp('图上显示检测到的八个基频的DFT幅度'); ％ 信号处理工具箱提供的goertzel函数能够取代 带下划线的三行程序。71程序hc941的说明图 9.4.2 接收端双频信号在八个近似基频上的 DFT 幅度72在双音拨号中的应用（第三段程序，选出两个峰值基频，查数字） limit = 80; % 设定门限，这个数与N有关 for s=5:8; if val(s) & limit, break, end % 查找列号 end for r=1:4; if val(r) & limit, break, end % 查找行号 end disp(['接收端检测到的号码为',setstr(tm(r,s-4))]) % 显示接收到的字符73在双音拨号中的应用在信号处理工具箱版本5.1中，也给出了Goertzel 算法的子程序goertzel.m。它的调用格式与gfft 相似，不同在于第二变元k可以是数组。所以 编程时可以省去一个for循环。另外它的k指的 是fft的下标，所以它等于频率的fft编号加一。 因此在程序hc941中，有波纹线的三行for循环 语句可以换成语句 X=goertzel(x(1:N),k+1)。 第三段程序中门限值limit的设定与Goertzel算法所 得的|H(k)|应得值有关，理想值是N/2，实际应 略小些。若取N不同，此门限也必须修改。749.5 正余弦信号的谱分析数字信号处理的一个重要用途是在离散时间域中 确定一个连续时间信号的频谱，通常称为频谱 分析。如果连续时间信号ga(t)是限带的，那末 它的离散时间等效物g(n)应当能给出ga(t)频谱 的近似的估计。 然而，这种近似取决于g(n)在多大程度上与ga(t) 相近，在大多数情况下，ga(t)是在-∞&t&∞范 围内定义的，要估计一个无限长信号的频谱是 不可能的，g(n)必须有有限的数据长度，且定 义在有限范围内。这就要涉及很多因素：如采 样的周期T，截取信号的长度N等。75正余弦信号的谱分析假定表征正余弦信号的基本参数，如振幅、频率 和相位，不随时间故变，则此信号g(n)的付立 叶变换G(ejω)可以用计算它的DTFT得到G ( e j? ) ?n ? ????g (n)e ? j?n实际上，把g(n)首先乘以一个长度为N 的窗函数 w(n) ， 使 它 变 成 一 个 长 为 N 的 有 限 序 列 ， g1(n)=g(n)w(n)，对g1(n)求出的DTFT G1(ejω)， 把它作为原连续模拟信号ga(t)的频谱近似估计。76正余弦信号的谱分析为了快速计算，要用DFT来求DTFT。故在 0≤ω ≤2π 区间等分为R点，求出离散付 立叶变换DFT。为了保证足够的分辨率， DFT的长度R选得比窗长度N大，其方法 是在截断了的序列后面补上R-N个零。计 算采用FFT算法。 我们更详细地考察一下上面的方法。特别 要分析加窗的效果，以及和由DFT样本来 估计DTFT频率采样值的问题。77正余弦信号的谱分析在讨论由G1(k)来估计频谱 G1(ejω) 和 G(ejω) 时，需要 探讨一下这些变换和它们所对应的频率之间的 关系，R点的DFT G1(k)与它的DTFT G1(ejω)的关 系为：G 1 (k) ? G 1 (e j? )? ? 2? k/R0≤k≤R-1(9.5.2)(9.5.3) (9.5.4)78归一化的数字频率ω k和DFT样本序号k的关系为?k ?2? k R模拟角频率和k的关系为(其中T是采样周期)2? k ?k ? RT正余弦信号的谱分析首先考虑单频率正余弦序列的频域分析。设一个 具有数字角频率ω 0的余弦信号为： g (n) ? cos(?0 n ? ? ) (9.5.5) 把这个序列表为指数序列： 查表得知它的DTFT为：G (ej?1 j (?0n?? ) g (n) ? e ? e ? j (?0n?? ) 2) ????(9.5.6)因此， 是一个以2π 为周期的ω的周期信号， 每个周期中在± ω 0处有两个冲击信号。79l ? ?? G(ejω)?e ? ? (? ? ? ?? j? 2? l ) ? e ? j? ? (? ? ? 0 ? 2? l ) 0?正余弦信号的谱分析取g(n)的一个有限长序列g1 (n) ? cos(?0 n ? ? ) 0 ? n ? N ?1假如g1是一个频率为10Hz而采样频率为64Hz的32 点序列，则用计算其DFT的程序如下： N=32; f=10;Fs=64l;n=0:N-1; g1=cos(2*pi*f*n/Fs);k=n; G1=fft(g1,N);stem(k,abs(G1),’.’)% 画出DFT G=fft(g1,1024);plot(N*…) % 画出DTFT曲线 计算结果见图9.5.1上图。它的DFT只有两个点不 等于零，位于k=5和k=27处。80正余弦信号的谱分析81正余弦信号的谱分析在上图中，k=5对应于频率10Hz， k=27对应于频 率54Hz（也就是-10Hz），这样DFT确实正确 地分辨了余弦信号的频率。但是这样理想的结 果是碰巧得到的，因为我们恰好截取了五个完 整的余弦周期(f*N/Fs=5)。如果截取的不是整 数周期，情况就不同了。 例如把频率f改为11Hz，而采样频率仍为64Hz和 窗长度仍为32点，用同样的程序计算此余弦信 号的频谱，则计算结果见图9.5.1下图。频谱图 上k=5和k=27处都有较大的峰值，而其它的点 上幅度不再为零。这样的现象称为频率泄漏， 来源于截断效应。82正余弦信号的谱分析从图9.5.1中可以看到，上下两个子图的DTFT形 状是很相像的。两个DFT样本点所以有那末大 的差别，原因就在于上图的采样点位置正巧都 在频谱的零点。实际工程中，输入信号的参数 是未知的，而且通常含有丰富的各种频谱，这 种理想情况不会出现。 为了理解图中的频谱形状，可以把有限序列g1(n) 看作无限序列g (n)和长度为N的矩形窗序列的 乘积。乘积的频谱应该等于乘子频谱的卷积。 结果，序列g1(n)的DTFT G1(ejω )是窗函数wd(n) 的DTFT WdR(ejω )移频±ω O并加权后的频谱的 的和。83正余弦信号的谱分析? 长度为N的矩形窗序列wd(n)的DTFT可从图 3.2.2推论得到。其DTFT WdR(ejω ) 形状如下图? 所以移频±ω O叠加就得到前面截断正余弦信号 的频谱。84正余弦信号的谱分析对于频率为11Hz的序列，它的数字频率为 ω0=11/64*2π=0.344π，因序列长度N=32，频率 分辨率为Δω=2π/N≈0.0625π。ω0 /Δω=5.5,因此 ω0在k横坐标上是在5和6之间。DFT是由长度 32的矩形窗的频谱，向左右各移动5.5后再相加 并乘1/2而得的。 在数字角频率0~2π(数字下标k=0~31)范围内，存 在着两个尖峰，一个在k=5.5处，表现为k=5和 k=6之间两处的DFT，另一个在k=32-5.5= 26.5处，表现为k=26和k=27处的两个DFT。所 有其它的DFT样本是由窗函数的泄漏引起的 DTFT的旁瓣所造成。85正余弦信号谱分析算例9.5.2当输入信号有一个以上的正余弦分量时，上述问 题变得更加严重。 例：9.5.2 考察DFT的长度对双频率信号频谱分 析的影响。设待分析的信号为 X(n)=0.5 sin(2π f1n)+sin(2π f2n) 0QnQN-1 令两个长度为N=16的正余弦序列的数字频率为 f1=0.22及f2=0.34。 解：现在计算x(n)在取不同R值时的DFT值。为此 采用下列程序，其中R,N,f1,f2都是可设定的。86正余弦信号谱分析算例9.5.2程序hc952 N=input('信号长度'); % R=input('DFT 长度'); % fr=input('[f1,f2]='); % 输入两个频率 n=0:N-1; % 设定自变量向量 x=0.5*sin(2*pi*n*fr(1))+sin(2*pi*n*fr(2)); % 两个正弦信号合成 X=fft(x,R);k=0:R-1; % 求其DFT stem(k,abs(X),'.');grid on % 画图87正余弦信号谱分析算例9.5.2取R为四个不同值16,32,64,128，画出的四个DFT 样本幅度图|X(k)|如图9.5.2上面四个子图。 从第一个子图上，很难看出它有两个峰值。因此 就要提高它的分辨率。把DFT的长度R由16增加 到32，64和128,逐渐可以看出它有两个峰值。 注意这几个图的横坐标为k，可换算为数字频 率f=ω/2π=k/R。这样就能确定峰值的位置大体 在f=0.21和0.35附近。 另外，增加DFT的长度R减小了相邻样本间的频 率间距，提高频谱的视在分辨率，因而可以提 高样本位置的测定精度。88正余弦信号的谱分析图 9.5.2 双频率 fr=[0.22,0.34] 有限长正余弦信号的 DTFT 幅频特性曲线89正余弦信号谱分析算例9.5.2上例中两对冲击函数，本来应该可以分辨的，但 由于和有限长的窗函数的频谱进行卷积，形成 了较宽的频谱图形。所有的窗函数主瓣宽度都 与N成反比。因此增加序列长度N可以有效地 提高频谱的实际分辨率。 虽然矩形窗的主瓣宽度最小，但它的旁瓣幅度太 大，造成严重的频率泄漏，这会使频谱分析的 可靠性和精确性下降。因此，人们宁可选一个 没有旁瓣的窗函数，而靠加大N来提高频谱的 分辨率。下面的例子说明用Hamming窗对改善 分辨率的影响。90正余弦信号谱分析算例9.5.2例9.5.3 在上例中若把两个正弦波的频率取得较 近，令fr=[0.22,0.25]，试问应该怎样选参 数才能在频谱分析中分辨出这两个分量？ 解：要使这两个频率之间隔开一个样本，分辨率 至少应达到Δf=0.03/2=0.015。因为此处的数字 频率是对采样频率Fs进行归一化的，即fr的最 大值为1。因此总的样本数N至少要达到 1/0.015=66。为了分析加窗的影响，只要在程 序hc952a中增加一句x1=x.*hamming(N)’，就构 成hc953的核心语句，然后对x1求DFT即可。91正余弦信号谱分析算例9.5.2程序hc953选择不同N,R和窗函数运行后，可得到 图9.5.3。从图上可以看出，信号的两个峰值是 能够分辨的。 将横坐标放大，可以得到图9.5.4。左边的一列是 加矩形窗的结果，右边的一列是加哈明窗的结 果，可以看出，使用无旁瓣的窗函数得到的频 谱函数比较光滑，便于分辨峰值位置和准确的 数值。所以在频谱分析仪中都要采取加窗的措 施。此外，为了提高实际的分辨率，应该尽量 增加信号长度N及DFT长度R(≥N)，92正余弦信号的谱分析图 9.5.3 双频率 fr=[0.22,0.25]正弦信号的 DFT 和 DTFT 幅频特性曲线93正余弦信号的谱分析图 9.5.4 把图 9.5.3 的横轴放大后的 DFT 和 DTFT 幅频特性曲线94非平稳信号的谱分析在工程实践中，往往不可能任意增加数据长度N。 人们所遇到的大多数信号，例如语音通信、雷 达信号…等等，都不是平稳的。在这种情况下， 增加取样的长度N未必有实际的意义，很多的 时候，要求人们进行非平稳信号的频谱分析， 于是出现了多种多样的分析方法。 这些都要用比较高深的数学工具，这里简单介绍 一种，称为短时间傅立叶变换(STFT－ShortTime Fourier Transfom)方法。 。95非平稳信号的谱分析序列x(n)的STFT定义为X STFT (e j? , n) ?其中，wd(m)是一个适当选择的窗函数序列。 这个式子说明，窗函数wd(m)在x(n)中选择一个 有限长度序列。在这段时间内，近似把非平稳 过程当作平稳过程来进行傅立叶分析。 对于一个非平稳序列，如果窗函数随时间移动， 则得出的频谱也将随时间而变化。对于随时间 而变化的频谱.时间轴和频率轴就占了两维，表 示其幅度的第三维只能用颜色或图形的灰度来 表示。96m?????x(n ? m) wd (m)e ? j?m(9.5.12)非平稳信号谱分析算例9.5.4程序hc954可以给出一个线性调频(Chirp)信号的 STFT频谱图，STFT函数specgram的缺省输入 变元为，FFT长度R=256，序列长度N=R/2，窗 函数为汉宁窗。 t=0:0.001:2; % 1kHz采样频率，持续2秒 x=chirp(t,0,1,150); % 1秒内从直流达到150Hz的线性调频信号 subplot(2,1,1),plot(t,x) % 画出信号曲线 subplot(2,1,2),specgram(x); % 显示频谱图 运行此程序得到图9.5.5。97非平稳信号谱分析算例9.5.4上图为调频信号时 间曲线。下图中横 坐标为时间轴，表 示窗的位置；纵坐 标为频率，向上为 频率增加；灰度表 示幅特性，全白意 味着幅度最大。向 上的白色斜线说明 峰值频率随时间的 增加向高频端移动。图 9.5.5 线性调频信号（上图）和它的 STFT 时频特性图989.6 音乐信号处理由于CD、DVD和家庭视听技术的发展，音乐的 录制和加工已经愈来愈多地采用数字技术。音 乐的录制大致分以下几个步骤： 首先在一个隔音的舞台上把乐队中各个乐器的声 音分别录在一个多磁道磁带的各个独立磁道上； 然后由音响工程师把各个磁道上的信号进行单独 的处理，加入特定的声音效果； 最后在一个混音系统中把这些信号进行最后的合 成，录制在一个立体声的双磁道磁带上。 这些特定的声音效果已愈来愈多地靠数字信号处 理技术来实现，本节将作一简单的介绍。99音乐信号处理时域处理方法 在一个像音乐厅那样的封闭空间中，人们听到的 声音包括直接传播声音、一次反射声音和混响 (reverberation)等几种成分。一次反射声音主要 是由直接传播声音在近距离内的直接回声构成， 而混响则是由密集的回声形成的。 所以在一个隔音的舞台上录制的音乐与音乐厅中 录制的效果就不一样，人们听起来会感到‘不 自然’。这时就需要用数字滤波器来人为地改 变录制的信号，通过增加一些回声，使它接近 于音乐厅中的效果。100音乐信号处理回声可以用延迟单元来生成。直接声音和 它的延迟了R个周期的单个回声可以用如 下的差分方程表示： y ( n) ? x ( n) ? ? x ( n ? R ) ? ? 1 (9.6.1) 其中α&1表示回声的衰减系数。上述差分方 程也可以用传递函数为: H ( z) ? 1 ? ? z ? R (9.6.2) 的FIR滤波器来实现，这实际上是一个梳状 滤波器，其结构如图9.6.1(a)。101音乐信号处理可以用实验来证明这一点。将‘大家好’文件中 的变量x作为输入语声x(n)，让它通过(9.6.2)的 滤波器，把得到的y(n)放音检验，为了得到0.3 秒 的 延 迟 和 30% 的 衰 减 ， 滤 波 器 参 数 R 选 为 0.3Fs=0.3×，α =0.3。 MATLAB程序为hc961如下： load dajiahao x1=[x;zeros(5000,1)]; %把x加长 y=filter([1,zeros(1,],1,x1); ％滤波输出 sound(y,22050) ％ 发声 为了使输出时间足够长，把x加了5000个样本。102音乐信号处理多重回声可以用一个IIR滤波器来实现。其传递函 数为 1 H ( z) ? ? ?1 ?R (9.6.4) 1?? z 在上述程序第二行中，把滤波器的分子分母系数 掉换一下，同时把α 取大一些，又把x加得更 长，写成以下求合成声音y1的语句 x1=[x;zeros(50000,1)] y1=filter(1,[1,zeros(1,],x1); sound(y1) 这样就可以检验多重回声效果。103音乐信号处理(b) 图 9.6.1 音乐回声信号仿真的两种滤波器结构图示为这两种梳状滤波器的结构。不过它们还无 法满足音乐处理的要求。首先是因为这种滤波 器的幅频特性不是常数，对不同频率的声音谐 波响应不均匀；其次是这种回声太单调，每秒 中的回声数目太少会导致声音的颤动。104音乐信号处理频域处理方法 把分别录制的各种乐器或歌手的声音进行混合时， 通常要由音乐工程师修改它们的频率响应。方 法是让信号通过一个均衡器（equalizer），其 作用是使这些声音在混合信号中充分表现。另 外也需要通过‘扩大’或‘削减’在某些频率 范围的信号，以修正低频和高频的信号之间的 关系。 通常用许多一阶和二阶的参数可调的滤波器级联 起来实现这个功能。滤波器的结构选择的主要 要求之一是调整方便，最好是调一个参数只影 响一个应用指标，而且可调参数要少。105音乐信号处理一阶均衡器的可调滤波器结构图如图9.6.1(a)，它 到两个输出端的传递函数分别为H LP ( z ) ? 0.5[1 ? A1 ( z )] H HP ( z ) ? 0.5[1 ? A1 ( z )]? ? z ?1 A1 ( z ) ? 1 ? ? z ?1 ? ?1(9.6.6)其中106音乐信号处理若把低频输出乘以K与高频输出相加，如右图， 就成为低频均衡滤波器；其传递函数为G1 ( z ) ? 0.5 ? K ? [1 ? A1 ( z )] ? 0.5 ? [1 ? A1 ( z )] ? 0.5( K ? 1) ? 0.5(1 ? K ) A1 ( z )(9.6.8)在α 和K给定时,其频率响应可用程序hc962计算:K=2.5; alpha=0.7; b1=[alpha,-1]; a1=[1,-alpha]; b = polyadd(0.5*(1+K)*a1,b1); [H,w]=freqz(b,a1);设定若干个alpha值，得出的曲线如下左图。107音乐信号处理设定的K值小于1时，得出的曲线如下右图。因此，可以根据需要，设置可调的电位器，来 提高低频的增益或高频的增益，并可调节边界 频率。108音乐信号处理二阶均衡器也采用图9.6.1的结构图，不同 点仅仅在于把其中的A1(z)换成A2(z)。? ? ? (1 ? ? ) z ?1 ? z ?2 A2 ( z ) ? 1 ? ? (1 ? ? ) z ?1 ? ?z ?2 ? ?1把A2(z) 替换了A1(z)后，公式(9.6.6)到(9.6.9) 全部适用。因此可以写出 H 2 ( z) ? 0.5 ? K ? [1 ? A2 ( z)] ? 0.5 ? [1 ? A2 ( z)] (9.6.11) 它的频率响应见图9.6.5。109音乐信号处理二阶均衡器的可调参数有三个：K,α和β。可以分 别调节幅度、带宽和中心频率。通常不用更高 阶的均衡器，宁可用多个均衡器级联。1109.7变采样率数字滤波前面所讨论的信号处理的各种方法都是把采样率 Fs视为固定值，即在一个数字中只有一个采样 频率。但在实际系统中，经常会遇到采样率的 转换问题，即要求一个数字系统能工作在“多 采样率”状态。在数字系统愈来愈普及的情况 下，各个数字系统都有自己不同的标准，它们 之间的衔接也不可避免地要遇到这个问题。 近年来，建立在采样率转换基础上的“多采样率 数字信号处理”已成为数字信号处理学科中的 主要内容之一。111变采样率数字滤波直观地想，变换采样率的最简单方法是： 先将以采样率Fs1采集的数字信号进行D／A 转换，变成模拟信号；再按采样率Fs2进 行A／D变换，这就实现了从Fs1到Fs2的采 样率转换。但这样做系统复杂，且易使 信号受到额外的干扰和崎变，所以在实 用的变采样率系统中，改变采样率并不 经过模拟信号，而完全是在数字域实现 的。112变采样率数字滤波采样率转换通常分为“抽取(Decimation)”和“插 值(Interpolation)”。抽取是降低采样率(所以在 英文中也称为Down-sampling)，以去掉多余数 据样本的过程，而插值则是提高采样率(Upsampling)以增加数据样本的过程。 变采样率是一个非线性问题，它涉及相当复杂的 数学推导，通常不在本科阶段讨论它的设计问 题。本书主要介绍抽取和插值的一般概念，一 方面可以加深对采样定理的理解，同时可以帮 助读者了解更多的应用问题，开扩思路。113内插过程的数字滤波信号的整数倍内插 整数倍内插是在已知的相邻两个采样点之间插 入L-1个样本点。由于这些样本点并非已知的 值，所以关键问题是如何求出这L-1个样本值。 前提是要直接在数字域完成这个插值。 用数字方法实现整数内插的步骤如下：在已知采 样序列x(n1T1)的相邻两个样点之间等间隔插入 L-1个零值样本点，得到输出y，则有? x(n / L), y(n) ? ? ?0 n ? 0,? L,?2L, ? n ? 其它值114内插过程的数字滤波示例式中L为大于1的整数，称为插入因子。按上式 y(n)是取值跳动的序列，让它经H(z)进行低通 数字滤波，得到教平滑的内插结果。这种内插 方案的框图如图9.7.1所示。各点的波形和频谱 见图9.7.2。xL↑yH(z)y1图 9.7.1 用数值方法提高采样率的框图115内插中的数字滤波示例9.7.1Nx=10; L=3; fsx=2; T=1/ alpha=0.5; nx=0:Nx-1; % 设定x的自变量向量 x=exp(-alpha*nx*T); % 给出x序列的值 Ny=Nx*L;ny=0:Ny-1; % 设定y的自变量向量 p=ones(1,L*Nx); % 生成p序列 y=zeros(1,L*Nx); % 先y序列初始化（全置零） y(1:L:L*Nx)=x; % 内插给出y在x对应点处的值 b=[0,0.3,0.9,0.2,0. 2,0.9,0]; %用例733的滤波器系数 y1=filter(b,1,y); % 进行滤波116内插中的数字滤波示例9.7.1y2=exp(-alpha*ny*T/L)/L; % 理想插值 X=fftshift(fft(x)); % x频谱，移到对称位置 Y=fftshift(fft(y)); % y频谱，移到对称位置 P=fftshift(fft(p)); % p频谱，移到对称位置 Y1=fftshift(fft(y1)); % y1频谱，移到对称位置 Y2=fftshift(fft(y2)); % y2频谱，移到对称位置 nxm=floor(nx-Nx/2+0.5); % 生成X频率位置序列 nym=floor(ny-Ny/2+0.5); % 生成Y频率位置序列 %以下绘图语句，得出各组曲线如下图。117变采样率数字滤波图 9.7.2 整数倍内插处理过程中的序列时域波形和对应的频谱特性118内插中的数字滤波示例9.7.1现在分析整个整数倍内插过程中的信号变换关系。 在由x变为y时，波形没有变，但采样率提高了， 这相当于把序列x与一个单位样本序列p相乘， 该p的采样频率为Fsy=L×Fsx，生成的y序列采 样频率与p相同，每隔L-1个点有一个非零脉冲， 其值等于该时刻的x(n)。 时域的相乘对应于频域的卷积，因此y的频谱Y相 当于把原来x的频谱X与频谱P卷积，频谱P是把 奈奎斯特频段平分的L个脉冲，时域的相乘对 应于频域的卷积，因此y的频谱Y相当于把原来 x的频谱X与频谱P卷积，频谱P是把x的奈奎斯 特频段扩展L倍的冲激频谱。119抽取过程的数字滤波信号的整数倍抽取 整数倍抽取是在已知序列x中每隔M-1个采样点 取出一个样本点，组成新的序列y。因此y序列 中的所有值都是已知的，可以直接写出下列表 达式 y (n y ) ? x ( M ? n y ) 也就是说，输入序列x中下标为M的整倍数的点 都被取入y序列，而其它的 样本都被舍去不用 了。用MATLAB语言来表达，为： Ny=ceil(Nx/M);nx=0:Nx-1; y=x(0:M:Nx); ny=0:Ny-1;120抽取过程的数字滤波抽取的数学关系看起来较插入简单。但由于采 样率的降低，造成的频率泄漏可能变得相当严 重。要解决好这个问题，在进行抽取之前，必 须对原来的序列进行预先数字滤波，称为‘抗 泄漏滤波’。它的框图可表为图9.7.3xH(z)x1M↓y1图 9.7.3用数值方法降低采样率的框图121抽取中的数字滤波算例9.7.2与框图对应的程序为hc972，先用滤波器滤波得 到x1，再抽取得到y1。故其核心语句为： nx=0:Nx-1; x=exp(-alpha*nx*T); % 给出x序列 b=[0,0.3,0.9,0.2,0. 2,0.9,0]; %例733的滤波器系数 x1=filter(b,1,x); %进行滤波，得出x1 Ny=ceil(Nx/M); ny=0:Ny-1; y1=x1(1:M:Nx); % 实行抽取，得出y1序列 以下分别画出x,x1,y1,y2的频谱见下图。122变采样率数字滤波图 9.7.4整数倍抽取处理过程中的序列时域波形和对应的频谱特 性123抽取中的数字滤波算例9.7.2从图中左边一列看出抽取过程的时域波形，指数 下降的输入信号x经过低通滤波成为波形x1， 然后经过抽取而降低了采样率,成为y1.最下面 的图表示不经过滤波直接抽取的结果。 再从右边的频谱图上可以看出滤波的必要性。直 接抽取波形的频谱abs(Y2)在奈奎斯特频率处 （也就是频谱图的左右边界±π 处）有很高的 幅度，这意味着存在着很大的频率泄漏。而经 过滤波再抽取的信号频谱abs(Y1)在该处的幅度 就小得多，如果用一个好一点的滤波器，结果 会更加理想。124分数倍变采样率数字滤波分数倍变采样率和相应的MATLAB函数 如果变换前后的采样频率的比满足两个整数之比 L/M，则不难想象，可以通过先作L倍内插后 作M倍抽取方法来实现，也可以先抽取后内插。 其滤波器的设计当然有不同的要求和方法。 MATLAB信号处理工具箱提供了一个变采样率函 数resample，用它可以把变采样率的全过程在 一步中完成。这个变采样率函数的典型调用格 式为y=resample(x,L,M)。其中L为采样率提高 的倍数，M为采样率降低的倍数，L和M必须 取整数。125变采样率数字滤波的应用整数倍抽取和内插的应用举例 在数字电话系统中，采样频率取为F＝8kHz，希 望传输尽量接近4kHz的音频带宽。但送话器发 出的信号x（t）的带宽比4kHz大很多。因此， 在A／D变换之前要对其进行模拟预滤波，以 防止采样后发生频率混叠失真。这就要求该滤 波器在4kHz处幅频特性很低，而为了使能传送 的信号频带尽量宽，又要求该滤波器的通带接 近于4kHz(例如3.8kHz)，给过渡带宽留的裕度 只有几百Hz，这样的要求对用模拟低通滤波器 是很难做到的。126变采样率数字滤波的应用为了降低对模拟预滤波器的技术要求，新方案先 用较高的采样率进行采样，比如取采样率Fsx ＝16kHz，这样，模拟预滤波器的过渡带可以 从3.8kHz到8kHz，这就容易设计了。经过A/D 后，再经M=2倍抽取，把采样率降低至8kHz。 此时问题的难点变成了抽取时的滤波器设计技 术。因为这是是在数字域的滤波，用FIR结构 的数字滤波器不难设计成线形相位和过渡带很 窄的特性。这种方案最终可有效地利用通信带 宽，并增加信号数据传送量。1279.8稀疏天线阵列设计本书讨论的信号处理都是以时间作为自变量的。 而在本例中自变量却是一维的空间位置。从本 例可以看到，利用数学模型的相似性，时域的 信号处理理论和方法可以推广到空域的信号处 理中。 可以发现，等间隔分布的线性相位天线阵元所产 生的远场方向图和FIR滤波器的频率响应在数 学模型上非常相似。这就使我们可以根据方向 图的要求来设计天线阵列。在实际中往往需要 去掉其中某些阵元，称为稀疏阵列。128稀疏天线阵列设计图9.8.1所示为N+1个等间隔配置的方向性相同的 阵元。阵元间距为d，因此它们在基线上排列 的坐标为x(n)=n×d，0≤n≤N。在与基线垂直 方向的远场辐射方向图可表为P(u ) ?n ?0? h ( n)eNj[ 2? ( u / ? ) d ]n,其中h(n)第n 个阵元的复数激励源，也称为复数 权函数。而u=sinθ ，θ 为与基线的法线方向的 夹角。如图所示。129稀疏天线阵列设计可以把P(u)看 作h(n)与变量w ? 2? (u / ? ) d构成的离散傅 立叶变换。对 于均匀激励的 阵列，h(n)=常 数a,通常选择 d=λ /2，这时 w的取值范围 为-π ～π 。130稀疏天线阵列设计从(9.8.1)式看出，P(u)的表示式与长度为N+1的 FIR滤波器的频率响应的表示式相同。因此FIR 滤波器的设计方法可以用于设计阵列天线阵。 因为最常见的情况是h(n)=1，它就和一个矩形 窗的频率响应相仿。 如果从中去掉某些阵元就形成了稀疏阵列。这时 某些阵元间的间距将大于λ /2,从而导致旁瓣 幅度的增加和栅格波瓣的出现。为了减小这些 不需要的副瓣，可以适当地调整收发阵元的排 列位置。131稀疏天线阵列设计有效方向图是由发射方向图与接收方向图的卷积 形成的，设计这样的方向图就比较简单些。 如果天线阵是由一个发射阵元和16个非稀疏的接 收阵元组成，那实际上仍然等价于非稀疏阵 列。，这时共需要17个阵元。如果用稀疏接收 阵列（或稀疏发射阵列），就可以节省阵元的 数目。比如，令发射阵列有两个阵元，而接收 阵列有八个阵元，排列如下：hT (n) ? [1 1], hR (n) ? [1 0 1 0 1 0 1 0 1 0 1 0 1 0 1]其中，0表示去除的阵元。132稀疏天线阵列设计则有效方向图函数将由两者的卷积构成，即heff (n) ? hT (n) ? hR (n) ? [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]因此只用了十个阵元就达到了17个阵元的效果。 按照这个思路继续分析，可以找到更经济的方法， 如果用八个阵元稀疏排列如下hT (n) ? [1 1 0 0 0 0 0 0 1 1], hR (n) ? [1 0 1 0 1 0 1]或hT (n) ? [1 1 1 1], hR (n) ? [1 0 0 0 1 0 0 0 1 0 0 0 1]也能够达到同样的效果。133稀疏天线阵列设计算例9.8.1可以用MATLAB计算FIR滤波器频率响应的函数 来计算这个收发阵列的有效方向图。程序 hc981中的核心语句如下 hr1=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1]; ht1=[1,1]; % 给定稀疏收发阵元 he1=conv(ht1,hr1); % 有效阵列等于两者的卷积 [Hr1,w]=freqz(hr1,1,'whole'); % 稀疏阵的方向图 [He1,w]=freqz(he1,1,'whole'); % 有效阵的方向图 % 绘图语句 程序运行的结果见图9.8.2134稀疏天线阵列设计从图上看到，单独的稀疏阵列造成的副瓣较大。 以上方法的本质是用一个稀疏阵列来‘填补’ 另一个稀疏阵列的空位。其结果是产生一个等 效幅度全为1的激励源，相当于一根矩形窗。图 9..8.2稀疏天线收发阵列卷积形成的有效方向图（用 FIR 滤波器频谱算法）135稀疏天线阵列设计矩形窗的旁瓣是比较大的，会造成较大的栅格波 瓣。可以用类似的概念进行插补，使得收发天 线的等效激励源得到更平滑的分布，以减小栅 格波瓣。例如，给出如下的收发阵列：hT (n) ? [1 1 0 0 1 1 0 0 1 1], hR (n) ? [1 0 1 0 1 0 1 0 1]它们卷积的结果是一个三角形的有效激励函数 heff (n) ? [1 1 1 1 2 2 2 2 3 3 2 2 2 2 1 1 1 1]这就像三角窗那样，可以大大减小旁瓣。如果在 阵列最边缘的阵元可以采用小于1的加权函数， 则有效激励函数可以进一步平滑。1369.9结束语数字信号处理技术正以前所未有的速度在发展， 它的发展是多方面的，这里作一个扼要的介绍： 1． 首先是应用领域的拓宽。许多过去用模拟信 号处理的设备，现在都力图用数字信号处理来 代替。包括通信、语音处理、雷达和声纳信号 处理、广播和电视技术、生物和医学信号处理、 地震和地理信息、交通和工业控制等等。出现 了数以万计的用于信号处理的专用集成电路芯 片(ASIC for DSP)，使得数字信号处理的速度 飞速提高，而其成本却不断降低。表9.9.1给出 了DSP芯片应用的一些主要领域。137138139结束语? 2． 数字信号处理应用的发展提出了多种多样 的复杂功能以及愈来愈高技术指标，推动了新 的理论不断出现。本书讨论的只限于一维信号， 要处理图形和图像就涉及二维信号处理，虽然 它的基础是一维信号处理，但究竟更加复杂， 而且带来许多新的更难的问题；比如在三维空 间成像、虚拟现实、数字视频、数字电影等， 使信号处理进入了三维。 在理论方面，它早就涉及了统计理论、非线性 理论等，现在又出现了模糊集合、小波变换、 神经元网络、还有各种谱估计理论等。140结束语3． 信号处理应用的发展也推动了它的研 制手段和工具的发展。除了DSP芯片的 性能日新月异之外，作为一种科学计算 软件，MATLAB为这些数字信号处理系 统提供了快捷的计算和仿真工具，主要 为非实时的顶层设计使用。一般是用来 研究总体方案和算法的可行性。现在除 了信号处理工具箱(Signal)之外，它还有 许多工具箱是与信号处理密切相关的。141结束语例如符号运算工具箱(Symbolic)、系统仿真工具 箱(Simulink)、用于仿真的数字信号处理模块库 (DSP Blockset)、定点运算模块库(Fix-Point Blockset)、统计处理工具箱(Staatistcs)、图像处 理工具箱(Image Processing)，模糊集合工具箱 （Fuzzy Logic）、神经元网络工具箱(Neural Network) 、小波变换工具箱(Wavelet)等等。其 它还有通信系统工具箱（Communication）、 通信模块库（Communication Blocksets）等等， 也都是分析信号处理问题时可能用到的工具。142结束语在本课中，我们只着重于原理和算法的部分，所 以用MATLAB就足以解决问题了。从原理和算 法向DSP芯片过渡是一个相当麻烦的过程。除 了要了解硬件的知识外，在软件上从MATLAB 语言转到芯片的汇编语言也要费很大的功夫。 为了解决这个问题，Mathworks公司提供了与TI 公司和Motorola公司的DSP芯片的接口软件。 其设想是，用MATLAB和Simulink验证过的系 统，通过这个软件编译成为特定规范的C语言， 它能够再编译为特定的DSP芯片的汇编语言。 这两个工具箱分别称为Developer’s Kits for TI DSPs 和Motorola DSP Developer’s Kits。143《数字信号处理》课程结束 谢谢各位的认真学习！144
数字信号处理教程―MATLAB释义与实现第09章―汇集和整理大量word文档,专业文献,应用文书,考试资料,教学教材,办公文档,教程攻略,文档搜索下载下载,拥有海量中文文档库,关注高价值的实用信息,我们一直在努力,争取提供更多下载资源。