《数字信号处理教程-MATLAB释义与实现》第九章

第九章

综合应用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中与信号处理有关的各种工具箱函数。59.1信号处理工具介绍

(讲本节时，应尽量直接用MATLAB的图形界面)信号处理工具sptool（SignalProcessingTool的缩写）为信号处理的研究工作提供了一个集成环境和工具。信号处理的任务：一是对信号进行分析；二是滤波器设计。在这两个任务进行过程中，经常要做第三个任务，那就是要把信号加到所设计的滤波器中，看它的输出是否满足要求。Sptool把这三个任务集成在一起，加上了适当的管理功能，配以良好的工作界面，利用它可以大大提高研究工作的效率。

6信号处理工具(sptool)介绍 在MATLAB命令窗中，键入fdatool,得到如右图的界面。各栏中分别存入了系统中原来已保存的信号、滤波器和频谱的名称。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】，指定数据来源；如果来源是工作空间，中间栏的【WorkspaceContents】将显示工作空间中的全部变量；右边一栏指定导入的目标，即要说明导入的数据是作为信号、还是滤波器还是频谱。右栏的下部有【SamplingFrequency】和【Name】两个框需要用户填写。按照说明的方法，可依次将信号s和8.6节设计的滤波器输入sptool。

9信号处理工具(sptool)介绍10信号处理工具(sptool)介绍2.信号的时域和频域观测在图9.1.1上，选定sig7,点击【view】(观测),就进入【SignalBrowser】(信号浏览器)视窗。如图9.1.3所示。要想求得该信号的频谱，也先选定sig7，然后在频谱栏下点击【creat】,此时出现频谱观测器的界面，见图9.1.4。先要在左栏上部选定求频谱的方法，目前我们只学了一种，即选FFT.在下面的框中填入点数，例如1024，然后点击【Apply】,就出现了频谱的曲线。11信号处理工具(sptool)介绍12信号处理工具(sptool)介绍13信号处理工具(sptool)介绍3．

滤波器的观测和修改选定滤波器的名称filt7,点击其下方的【View】，就可以观测它的频率响应。如点击【Edit】，那就不仅能观测，还能修改。4.

让信号通过滤波器求输出信号在图9.1.1中，选定信号和滤波器名称，点击滤波器栏下的按钮【Apply】,这意味着将所选信号加到所选滤波器中去。此时将产生一个小视窗，提示用户输入信号和滤波器的名称，并提供输出信号缺省名称sig8。如无修改，点击OK，输出信号将自动生成并列写在信号栏中。14信号处理工具(sptool)介绍在本例中滤波器输出sig8基本上是一个较纯的10Hz正弦波，这是不难想象的。经过仔细观测和分析，如不满意，可以修改滤波器后重新实验。用这个工具，不难在很短的时间内完成多个滤波器的设计，并分别观测它们的输出。所有的结果都存储在案，可用于试验报告。所以这是一个提高效率的有效工具。159.2语音量化噪声的改善脉冲编码调制(PulseCodeModulation—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)其中，量化步长为△=2–N。量化误差的均方值为: (9.2.3) 17语音量化噪声的改善用分贝来度量的噪声均方值为: (9.2.4)可以看出，上述的量化器每增加一位，量化噪声减小6dB，高质量语音要求每个样本至少量化为12位，因此传送速率至少为96000位/秒。最大幅度为±V伏的N位（不含符号位）二进制A/D变换器的数学模型建立如下。它把电压V分解为2N-1份，故量化步长为V/(2N-1)，得出二进制量化子程序bqtize。

18语音量化的子程序functiony=bqtize(x,N,V)ifnargin<3V=max(abs(x));end %V缺省取最大xax=abs(x);%去掉符号deltax=V/(2^N-1);%求量化步长xint=fix(ax./deltax+0.5); %x的量化整数y=sign(x).*xint.*deltax;%恢复量化原值这个A/D变换子程序的输入是连续模拟电压x，输出则是量化了的模拟电压y。均匀量化器在信号的整个动态范围中的量化步长相同，所以量化噪声均方值不变。

19语音量化噪声的改善均匀量化器在信号的整个动态范围中的量化步长相同，所以量化噪声均方值不变。然而，语音信号的特性是小幅度比大幅度出现得频繁。对小信号而言，量化噪声使信噪比大大下降。解决的途径之一是用非均匀量化器。不过在技术上制造非均匀量化器的芯片是困难的。得到非均匀量化器特性的另一个方法是用压缩—扩张器。可先使信号通过压缩幅度的非线性器件，后面再接一均匀量化器，再用逆向扩张幅度的非线性器件恢复信号。如下图。

20语音量化的压缩扩张器压缩器的作用是把小信号放大，大信号缩小，所以把压缩器后的信号进行量化，小信号的信噪比就得到提高。而后通过扩张器，信号和量化噪声同时作非线性变换，信号复原为原来电平，小信号信噪比则保持较小的水平。21语音量化的压缩扩张器在我国的通信系统中使用的对数压缩器(称为μ律压缩器)具有如下输入输出幅度特性:其中，x是归一化输入，y是归一化输出，sign(.)是符号函数，μ=255,是控制压缩特性的参数。由量化数字信号恢复模拟信号时，逆关系μ律为: 22语音量化的压缩扩张器可以用MATLAB把压缩器和扩张器用函数程序mulawcom和mulawexp来表示，例如把压缩器函数程序列写如下：functiony=mulawcom(x,mu,V)ifnargin<3|V<max(abs(x))V=max(abs(x));endy=V/log(1+mu)*log(1+mu/V*abs(x)).*sign(x);对μ律扩张器和A-律的压扩器也可以按类似的方法写出相应的子程序。23语音量化压缩器的特性例9.2.1画出μ律和A律的压缩器输入输出曲线。解：写出如下的MATLAB程序hc921x=0:0.01:1;y=mulawcom(x,255,1);y1=Alawcom(x,87.56,1);plot(x,y,x,y1,':');gridonlegend('\mu律','A律')运行此程序的结果见图9.2.1。请注意二者非常相似。24语音量化压缩器的特性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-xq;er=e./abs(x);%求相对误差erplotyy(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语音量化的信噪特性对比28语音量化的信噪特性对比从左图中可以看出，在输入为-5~5V范围内，绝对量化误差e呈等幅锯齿波形式，其最大值恒定，因而其相对误差er在小的输入幅度x处急遽增大，超过了1。

由右图从中可以看出它的绝对误差随x的增加而增大，其相对误差则在整个输入范围内呈等幅锯齿波形式，最大不到0.2。因此量化误差yq-x1引起的信噪比将比不用压扩器有显著的提高。299.3系数和运算量化的影响要用计算机处理信号，就不可避免地要遇到量化问题。这里包括三方面:(1)信号的量化；(2)系统参数的量化；以及(3)信号与系统相互运算中的量化。在MATLAB中，参数和运算都有十六位十进制（64位二进制浮点数）的精度，量化误差通常可以忽略不计。当要把用MATLAB设计的滤波器付诸实现，特别是用于实时系统时，要采用嵌入式的数字信号处理芯片。它们通常只有较短的（8位、16位、32位）二进制定点数精度，必须考虑量化的影响。绝不可盲目轻信MATLAB中计算出的结果。30系数和运算量化的影响MATLAB按二进制双精度格式来表示数，表达一个数需要8个字节，也就是64个二进制位。按照IEEE标准，双精度浮点数η用下式表示其中M是一个小于一大于1/2的二进制分数，称为尾数，占用52个二进制位，近似16位十进制；而指数E是一个带符号的二进制整数。占11个二进制位，总共可表示2048个整数，即可以表示从-1023到+1024的数集，它决定了数的动态范围。数的正负号反映在S上，它只占一位。所以用MATLAB计算通常不必考虑量化误差。31系数和运算量化的影响系数量化：系数和信号的不同在于它的值并非来自A/D变换，而是由设计者在计算机内设定的，各个系数没有必要用同样的量化步长。通常应该尽量利用系统的硬件字长，把各个系数设置得相对量化误差为最小。本节将先给出系数相对量化的算式，然后考虑滤波器系数的量化误差对滤波器性能的影响。因为量化是一个非线性问题，很难有封闭形式的解，我们只能用一些数字例子说明在进行滤波器结构选择时应考虑的原则。

32系数和运算量化的影响系数相对量化的数学算式和子程序：十进制的相对量化，就是截取有效位的长度。若xd=[3.14159,0.38455,54.3753,-0.134273]，键入xdr2=chop(xd,2)和xdr3=chop(xd,3)，得到的结果分别为xdr2=3.100.3854.0-0.13Xdr3=3.140.38554.4-0.134xdr2与xd之间的最大误差为十进制第三位的±5。所以其最大相对误差的范围为±(0.5~5)%。

33系数和运算量化的影响要实现二进制数的相对量化，应把原量化子程序bqtize中的步长deltax取成输入量乘以二进制负幂，故它成为随各个输入量不同的数组。得到子程序brqtize.m

functiony=brqtize(x,N)ax=abs(x);%去掉符号m=ceil(log(ax)/log(2));%ax=2^mdeltax=2.^(m-N);%求出量化步长数组xint=round(ax./deltax);%量化整数y=sign(x).*xint.*deltax;%量化输出

34系数和运算量化的影响程序中采用了

(9.3.1)来确定数x的最高二进制位数m，其中表示向上取整。得到m，就可以知道该数的二进制上限，从而确定最小量化单位Δ=2(m-N)，参看式

(1.1.2)就可以理解此子程序中后两条语句。

35系数和运算量化的影响这个子程序的优点一是适用于x为数组的情况，不用一个一个系数去量化，特别方便。二是对所有的系数，其相对量化误差都以各个系数本身的绝对值来计算，最大相对误差都是2^-N。因此它可以对滤波器系数向量实现统一的量化。仍取前面给出的xd，将它量化至4位二进制。键入xdq=bqtize(xd,4)得到

xdq=3.250.375056.0-0.1406求其相对误差(xdq-xd)./abs(xd)，得到ans=0.0345-0.02480.0299-0.0473相对误差比较均匀，都在2^(-4)=0.0625以下。

36系数和运算量化的影响例9.3.2设一个四阶FIR滤波器，其传递函数为：算出它的级联结构，比较在两种结构下系数量化为四位二进制对频率响应的影响。解：先用[sos,G]=tf2sos(b,1)求其级联结构参数，分别求b,sos,G的量化值bq=brqtize(b,4);sosq=brqtize(sos,4),Gq=brqtize(G,4)得出的bq=0.07030.37500.3750…然后用这几个参数计算滤波器的幅特性，并与量化前参数产生的幅特性进行比较，得出如下曲线。37系数和运算量化的影响以理想情况作为标准，可以看出，在本例中，用级联结构时，系数量化对系统幅频特性的影响较直接结构为小。

38系数和运算量化的影响系数量化对IIR滤波器的影响举例IIR滤波器分母系数存在着递归计算的问题，它的系数误差可能在循环计算中不断扩散，因此，一般地说，对系数误差的敏感程度比FIR滤波器严重。特别是当IIR滤波器的极点靠近z平面上单位圆时，幅频特性通常要出现很大的峰值。系数的量化误差造成的影响很难用解析公式表示，只能通过实例来计算和观察。39系数和运算量化的算例9.3.2一IIR滤波器具有如下传递函数将系数量化为七位二进制，分析其极点的变化。解：系数量化前、后的极点由下列核心语句求得：aq=brqtize(a,7)%量化后分母系数向量

p=roots(a),%量化前的极点pq=roots(aq),%量化后的极点其量化前、后零极点的位置画在z平面上如下图。键入abs(pq),得[1.0079,1.0079,0.9613,0.5000]40系数和运算量化的影响 量化前的零、极点用‘o’,’x’表示，而量化后的零、极点用‘+’,’*’表示。从计算的数据看出，量化后极点的模大于1，出了单位圆，系统由稳定变为不稳定。41系数和运算量化的影响从这个例子看到，滤波器的级联结构对系数量化的敏感性较之直接结构要小。高次多项式系数的微小变化往往会引起根的显著变动。因此，阶次高的滤波器不宜采用直接结构。系数量化还会造成其它一些问题，如引起极限环振荡等，这些问题一般出现在量化位数较小的单片式嵌入数字系统中。当位数达到16位时，这类问题已经不很突出了。研究系数量化问题还没有完美的理论，主要靠针对具体问题的仿真计算来判别量化的影响。

42系数和运算量化的影响运算量化的影响运算精度的规律如下：两个数相加，输出的有效数位与输入大数的有效数位同；两数相减，输出的相对误差将增大，增大的倍数等于输入大数与差之比；两数相乘，输出的有效数位等于输入两数的最低有效数位。所以只要避免两个数量相近的大数的相减，运算器的量化误差与输入数据的量化误差数量级相同。信号运算的细致过程可以用信号流图来表示。加法器和乘法器如果工作正常（没有溢出），应该能够保持运算的量化精度。43系数和运算量化的影响加法和乘法器的问题一般不在于量化，而在于发生溢出。溢出后如果不加特殊处理，加法器的输出会从正的最大值跳变到负的最大值，或反之。这种情况有时会引起大范围振荡，彻底破坏滤波器的工作，那时，上面的量化模型全都无效了，所以这是绝不容许发生的。要避免这种现象发生，除了对输入量的动态范围进行控制之外，应该把加法器设置成溢出后饱和的工作状态，不让它在溢出后出现正负号的突跳。449.4在双音拨号中的应用双音多频(DTMF—DualToneMultiFrequency)是用按键进行电话拨号的体制。它不单用在电话中，还可以用于传输十进制数据的其它通信系统中。DTMF也广泛应用于电子邮件和银行系统。第一章1.2节中已经对它的功能进行了描述，4个低频频率表示四行，3个高频频率表示三列，两者的组合共可提供4×3＝12个字符。包括0到9的十个数字，加上字符*和#。本节要讨论如何解决双音变为数字的问题。45在双音拨号中的应用

DTMF通信系统是一个很典型的小型信号处理系统，它既有模拟信号的生成和传输部分（这要用到D/A转换）；又有把它转为数字信号（这要用到A/D转换）并进行数字处理的部分；而为了提高系统的检测速度和降低成本，还开发了一种特殊的DFT算法，称为Goertzel算法；这种算法在国外的几乎每一本数字信号处理的教材上都要介绍，说明了人们对它的重视。这种算法既可以用硬件（专用芯片）、也可以用软件实现，所以DTMF系统的设计问题是理论与工程相结合的一个很好的典范。46在双音拨号中的应用先研究双音信号的生成问题：DTMF的两个音频可以用计算法或查表法产生。用计算法得到正弦波形的缺点是要占一些运算时间；查表法的速度较快，缺点是要占一定的存贮空间；两个正弦波的数字样本按比例相加在一起。因为采样频率是8KHz，硬件必须每125毫秒输出一个样本。将这个叠合信号送到D/A变换器变换成模拟音频信号，通过电话线路传送到交换机。

47在双音拨号中的应用在接收端，将收到的模拟音频信号进行A/D变换，恢复为数字信号，然后检测其中的音频频谱来确定所发送的数字。检测算法可以用FFT，也可用一组滤波器来提取所需频率。当要检测的音频数目比较少时，用滤波器组更节省硬件。如果用FFT算法实现该DFT计算，计算量(复数乘法和加法)是Nlog2N。好处是立即得到DFT的所有N个值，但至少要N个存储器。然而，如果只希望计算DFT的K个点,而K<<N，则直接计算可以节省很多内存。下面介绍的Goertzel算法就属于后者。48DFT的Goertzel算法Goertzel算法利用相位因子Wnk的周期性，将DFT运算表示为线性滤波运算，由于WN-kN=1，可用该因子去乘DFT，则它是两个序列的卷积。一是长度为N的输入序列x(n)，另一则是脉冲响应序列 可以看作x(n)通过具有该脉冲响应的滤波器。该滤波器在n＝N点的输出就是频点k处的DFT.49DFT的Goertzel算法单位脉冲响应为hk(n)的滤波器的系统函数就是(9.4.4)式的z变换，容易求出为这个滤波器只有一个位于单位圆上频率为ωk=2πk/N处的极点。Goertzel算法的好处不在于节省时间，而在于节省空间。如果只需要K个DFT样本，可以只用K个并行的单极点滤波器来分别计算这K个样本，，这就可以大大节省硬软件资源。50DFT的Goertzel算法根据这个滤波器的差分方程，因为可以用迭代方法计算yk(n)，如下图。要执行该计算，可以只算一次相位因子WN-k，将其存贮起来。迭代N次后，预期的输出为：X(k)=yk(N)。51DFT的Goertzel算法上式中包含对 的复数运算，想避免它，可将成对的复共轭极点 组合在一起计算，导出双极点滤波器的系统函数:它可以看成用两个差分方程组的联立:

初始条件为

vk(-1)=vk(-2)=0。运算结构图见后：52DFT的Goertzel算法53DFT的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，故计算的最后一步还可以简化，从而完全避开复数运算。54Goertzel算法子程序按(9.4.7)式编成子程序gfft。它根据输入序列x，和规定的DFT样本序号k，计算DFT样本X。functionX=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);%为下一步求WX=v(N+1)-W*v(N);%求出第k个DFT55在双音拨号中的应用DTMF信号参数选择的考虑因素：(1).为了抗干扰，要检测8个频率及其倍频；(2).为了能分辨清这16个频率，要求采集数据的长度N>Fs/D，已知Fs=8000，而根据最小相邻基频间隔（73Hz）大于2D的要求，D<37，得到N应在200以上。而为了提高检测速度，又希望N尽量小些。(3).为了检测尽量精确，希望这16个离散频点靠近k=整数处，这需要仔细选择N。最后选定的规范值是N=205。56在双音拨号中的应用本节用一个MATLAB程序hc941来演示DTMF双频拨号的全过程。包括以下几部分：（1）首先程序要求用户输入一个数字，然后根据这个数字查出它对应的两个频率，并生成相应的双频信号x(n)，并产生声音作为标志。（2）在接收端，程序截取收到的双频信号并进行傅立叶变换，求出它在八个规定频率上的样本幅度，画出频谱图。（3）从频谱中分别取出满足规定电平的高低频分量的两个下标。再由下标找到它的ASCII码和数字。57程序hc941的说明此程序分为三段。第一段是双频信号的生成，从输入一个数字开始，以产生相应的双频信号结尾；第二段是接收端信号处理，从收到双频信号开始，进行截断和量化，并计算八个基频DFT，以画出DFT幅度分布图结尾；第三段是根据基频幅度分布图,截取其两个峰值，查找和显示输入字符。各段之间用pause命令隔开。58d=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码，检测输入码与谁相符forp=1:4;forq=1:4;iftm(p,q)==abs(d);break,end%检测列号qendiftm(p,q)==abs(d);break,end%检测行号pend程序hc941的说明59f1=[697,770,852,941];%行频率向量f2=[1209,1336,1477,1633];%列频率向量n=0:2040;%为了发声，加长序列%构成双频信号x=sin(2*pi*n*f1(p)/8000)+sin(2*pi*n*f2(q)/8000);sound(x);%发出声音disp('双频信号已经生成并发出'),pause％以上是第一段程序，表示发送端的功能。程序hc941的说明60程序hc941的说明（第二段）接收检测端的计算8点DFT的程序N=205; %截取信号的长度 k=[1820222431343842];%待求的DFT序号form=1:8;X(m)=gfft(x(1:N),k(m));%算八点DFT值endval=abs(X);stem(k,val,'.');%画出八点DFT向量disp('图上显示检测到的八个基频的DFT幅度');

％信号处理工具箱提供的goertzel函数能够取代带下划线的三行程序。61程序hc941的说明62在双音拨号中的应用（第三段程序，选出两个峰值基频，查数字）limit=80;%设定门限，这个数与N有关fors=5:8;ifval(s)>limit,break,end%查找列号endforr=1:4;ifval(r)>limit,break,end%查找行号enddisp(['接收端检测到的号码为',setstr(tm(r,s-4))])%显示接收到的字符63在双音拨号中的应用在信号处理工具箱版本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不同，此门限也必须修改。649.5正余弦信号的谱分析

数字信号处理的一个重要用途是在离散时间域中确定一个连续时间信号的频谱，通常称为频谱分析。如果连续时间信号ga(t)是限带的，那末它的离散时间等效物g(n)应当能给出ga(t)频谱的近似的估计。然而，这种近似取决于g(n)在多大程度上与ga(t)相近，在大多数情况下，ga(t)是在-∞<t<∞范围内定义的，要估计一个无限长信号的频谱是不可能的，g(n)必须有有限的数据长度，且定义在有限范围内。这就要涉及很多因素：如采样的周期T，截取信号的长度N等。65正余弦信号的谱分析假定表征正余弦信号的基本参数，如振幅、频率和相位，不随时间故变，则此信号g(n)的付立叶变换G(ejω)可以用计算它的DTFT得到

实际上，把g(n)首先乘以一个长度为N的窗函数w(n)，使它变成一个长为N的有限序列，g1(n)=g(n)w(n)，对g1(n)求出的DTFTG1(ejω)，把它作为原连续模拟信号ga(t)的频谱近似估计。66正余弦信号的谱分析为了快速计算，要用DFT来求DTFT。故在0≤ω≤2π区间等分为R点，求出离散付立叶变换DFT。为了保证足够的分辨率，DFT的长度R选得比窗长度N大，其方法是在截断了的序列后面补上R-N个零。计算采用FFT算法。我们更详细地考察一下上面的方法。特别要分析加窗的效果，以及和由DFT样本来估计DTFT频率采样值的问题。

67正余弦信号的谱分析在讨论由G1(k)来估计频谱G1(ejω)和G(ejω)时，需要探讨一下这些变换和它们所对应的频率之间的关系，R点的DFTG1(k)与它的DTFTG1(ejω)的关系为：

0≤k≤R-1 (9.5.2)归一化的数字频率ωk和DFT样本序号k的关系为

(9.5.3)模拟角频率和k的关系为(其中T是采样周期) (9.5.4)68正余弦信号的谱分析首先考虑单频率正余弦序列的频域分析。设一个具有数字角频率ω0的余弦信号为：

(9.5.5)把这个序列表为指数序列：

(9.5.6)查表得知它的DTFT为：

因此，G(ejω)是一个以2π为周期的ω的周期信号，每个周期中在±ω0处有两个冲击信号。69正余弦信号的谱分析取g(n)的一个有限长序列假如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处。

70正余弦信号的谱分析71正余弦信号的谱分析在上图中，k=5对应于频率10Hz，

k=27对应于频率54Hz（也就是-10Hz），这样DFT确实正确地分辨了余弦信号的频率。但是这样理想的结果是碰巧得到的，因为我们恰好截取了五个完整的余弦周期(f*N/Fs=5)。如果截取的不是整数周期，情况就不同了。例如把频率f改为11Hz，而采样频率仍为64Hz和窗长度仍为32点，用同样的程序计算此余弦信号的频谱，则计算结果见图9.5.1下图。频谱图上k=5和k=27处都有较大的峰值，而其它的点上幅度不再为零。这样的现象称为频率泄漏，来源于截断效应。72正余弦信号的谱分析

从图9.5.1中可以看到，上下两个子图的DTFT形状是很相像的。两个DFT样本点所以有那末大的差别，原因就在于上图的采样点位置正巧都在频谱的零点。实际工程中，输入信号的参数是未知的，而且通常含有丰富的各种频谱，这种理想情况不会出现。为了理解图中的频谱形状，可以把有限序列g1(n)看作无限序列g(n)和长度为N的矩形窗序列的乘积。乘积的频谱应该等于乘子频谱的卷积。结果，序列g1(n)的DTFTG1(ejω)是窗函数wd(n)的DTFTWdR(ejω)移频±ωO并加权后的频谱的的和。

73正余弦信号的谱分析长度为N的矩形窗序列wd(n)的DTFT可从图3.2.2推论得到。其DTFTWdR(ejω)

形状如下图所以移频±ωO叠加就得到前面截断正余弦信号的频谱。74正余弦信号的谱分析对于频率为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的旁瓣所造成。75正余弦信号谱分析算例9.5.2当输入信号有一个以上的正余弦分量时，上述问题变得更加严重。例：9.5.2考察DFT的长度对双频率信号频谱分析的影响。设待分析的信号为X(n)=0.5sin(2πf1n)+sin(2πf2n)0≦n≦N-1令两个长度为N=16的正余弦序列的数字频率为f1=0.22及f2=0.34。解：现在计算x(n)在取不同R值时的DFT值。为此采用下列程序，其中R,N,f1,f2都是可设定的。76正余弦信号谱分析算例9.5.2程序hc952N=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;%求其DFTstem(k,abs(X),'.');gridon%画图77正余弦信号谱分析算例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减小了相邻样本间的频率间距，提高频谱的视在分辨率，因而可以提高样本位置的测定精度。78正余弦信号的谱分析79正余弦信号谱分析算例9.5.2上例中两对冲击函数，本来应该可以分辨的，但由于和有限长的窗函数的频谱进行卷积，形成了较宽的频谱图形。所有的窗函数主瓣宽度都与N成反比。因此增加序列长度N可以有效地提高频谱的实际分辨率。虽然矩形窗的主瓣宽度最小，但它的旁瓣幅度太大，造成严重的频率泄漏，这会使频谱分析的可靠性和精确性下降。因此，人们宁可选一个没有旁瓣的窗函数，而靠加大N来提高频谱的分辨率。下面的例子说明用Hamming窗对改善分辨率的影响。80正余弦信号谱分析算例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即可。81正余弦信号谱分析算例9.5.2程序hc953选择不同N,R和窗函数运行后，可得到图9.5.3。从图上可以看出，信号的两个峰值是能够分辨的。将横坐标放大，可以得到图9.5.4。左边的一列是加矩形窗的结果，右边的一列是加哈明窗的结果，可以看出，使用无旁瓣的窗函数得到的频谱函数比较光滑，便于分辨峰值位置和准确的数值。所以在频谱分析仪中都要采取加窗的措施。此外，为了提高实际的分辨率，应该尽量增加信号长度N及DFT长度R(≥N)，82正余弦信号的谱分析83正余弦信号的谱分析84非平稳信号的谱分析在工程实践中，往往不可能任意增加数据长度N。人们所遇到的大多数信号，例如语音通信、雷达信号…等等，都不是平稳的。在这种情况下，增加取样的长度N未必有实际的意义，很多的时候，要求人们进行非平稳信号的频谱分析，于是出现了多种多样的分析方法。这些都要用比较高深的数学工具，这里简单介绍一种，称为短时间傅立叶变换(STFT－Short-TimeFourierTransfom)方法。。85非平稳信号的谱分析序列x(n)的STFT定义为

(9.5.12)

其中，wd(m)是一个适当选择的窗函数序列。这个式子说明，窗函数wd(m)在x(n)中选择一个有限长度序列。在这段时间内，近似把非平稳过程当作平稳过程来进行傅立叶分析。对于一个非平稳序列，如果窗函数随时间移动，则得出的频谱也将随时间而变化。对于随时间而变化的频谱.时间轴和频率轴就占了两维，表示其幅度的第三维只能用颜色或图形的灰度来表示。86非平稳信号谱分析算例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。87非平稳信号谱分析算例9.5.4上图为调频信号时间曲线。下图中横坐标为时间轴，表示窗的位置；纵坐标为频率，向上为频率增加；灰度表示幅特性，全白意味着幅度最大。向上的白色斜线说明峰值频率随时间的增加向高频端移动。

889.6音乐信号处理由于CD、DVD和家庭视听技术的发展，音乐的录制和加工已经愈来愈多地采用数字技术。音乐的录制大致分以下几个步骤：首先在一个隔音的舞台上把乐队中各个乐器的声音分别录在一个多磁道磁带的各个独立磁道上；然后由音响工程师把各个磁道上的信号进行单独的处理，加入特定的声音效果；最后在一个混音系统中把这些信号进行最后的合成，录制在一个立体声的双磁道磁带上。这些特定的声音效果已愈来愈多地靠数字信号处理技术来实现，本节将作一简单的介绍。89音乐信号处理时域处理方法在一个像音乐厅那样的封闭空间中，人们听到的声音包括直接传播声音、一次反射声音和混响(reverberation)等几种成分。一次反射声音主要是由直接传播声音在近距离内的直接回声构成，而混响则是由密集的回声形成的。所以在一个隔音的舞台上录制的音乐与音乐厅中录制的效果就不一样，人们听起来会感到‘不自然’。这时就需要用数字滤波器来人为地改变录制的信号，通过增加一些回声，使它接近于音乐厅中的效果。90音乐信号处理回声可以用延迟单元来生成。直接声音和它的延迟了R个周期的单个回声可以用如下的差分方程表示：

(9.6.1)其中α<1表示回声的衰减系数。上述差分方程也可以用传递函数为: (9.6.2)的FIR滤波器来实现，这实际上是一个梳状滤波器，其结构如图9.6.1(a)。91音乐信号处理可以用实验来证明这一点。将‘大家好’文件中的变量x作为输入语声x(n)，让它通过(9.6.2)的滤波器，把得到的y(n)放音检验，为了得到0.3秒的延迟和30%的衰减，滤波器参数R选为0.3Fs=0.3×22050=6615，α=0.3。MATLAB程序为hc961如下：loaddajiahaox1=[x;zeros(5000,1)];%把x加长y=filter([1,zeros(1,6615),0.3],1,x1);

％滤波输出sound(y,22050) ％发声为了使输出时间足够长，把x加了5000个样本。92音乐信号处理多重回声可以用一个IIR滤波器来实现。其传递函数为

(9.6.4)在上述程序第二行中，把滤波器的分子分母系数掉换一下，同时把α取大一些，又把x加得更长，写成以下求合成声音y1的语句x1=[x;zeros(50000,1)]y1=filter(1,[1,zeros(1,6615),0.8],x1);sound(y1)这样就可以检验多重回声效果。93音乐信号处理图示为这两种梳状滤波器的结构。不过它们还无法满足音乐处理的要求。首先是因为这种滤波器的幅频特性不是常数，对不同频率的声音谐波响应不均匀；其次是这种回声太单调，每秒中的回声数目太少会导致声音的颤动。

94音乐信号处理频域处理方法把分别录制的各种乐器或歌手的声音进行混合时，通常要由音乐工程师修改它们的频率响应。方法是让信号通过一个均衡器（equalizer），其作用是使这些声音在混合信号中充分表现。另外也需要通过‘扩大’或‘削减’在某些频率范围的信号，以修正低频和高频的信号之间的关系。通常用许多一阶和二阶的参数可调的滤波器级联起来实现这个功能。滤波器的结构选择的主要要求之一是调整方便，最好是调一个参数只影响一个应用指标，而且可调参数要少。95音乐信号处理一阶均衡器的可调滤波器结构图如图9.6.1(a)，它到两个输出端的传递函数分别为

(9.6.6)其中96音乐信号处理若把低频输出乘以K与高频输出相加，如右图，就成为低频均衡滤波器；其传递函数为

(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值，得出的曲线如下左图。97音乐信号处理设定的K值小于1时，得出的曲线如下右图。因此，可以根据需要，设置可调的电位器，来提高低频的增益或高频的增益，并可调节边界频率。98音乐信号处理二阶均衡器也采用图9.6.1的结构图，不同点仅仅在于把其中的A1(z)换成A2(z)。

把A2(z)替换了A1(z)后，公式(9.6.6)到(9.6.9)全部适用。因此可以写出

(9.6.11)它的频率响应见图9.6.5。

99音乐信号处理二阶均衡器的可调参数有三个：K,α和β。可以分别调节幅度、带宽和中心频率。通常不用更高阶的均衡器，宁可用多个均衡器级联。1009.7变采样率数字滤波

前面所讨论的信号处理的各种方法都是把采样率Fs视为固定值，即在一个数字中只有一个采样频率。但在实际系统中，经常会遇到采样率的转换问题，即要求一个数字系统能工作在“多采样率”状态。在数字系统愈来愈普及的情况下，各个数字系统都有自己不同的标准，它们之间的衔接也不可避免地要遇到这个问题。近年来，建立在采样率转换基础上的“多采样率数字信号处理”已成为数字信号处理学科中的主要内容之一。

101变采样率数字滤波直观地想，变换采样率的最简单方法是：先将以采样率Fs1采集的数字信号进行D／A转换，变成模拟信号；再按采样率Fs2进行A／D变换，这就实现了从Fs1到Fs2的采样率转换。但这样做系统复杂，且易使信号受到额外的干扰和崎变，所以在实用的变采样率系统中，改变采样率并不经过模拟信号，而完全是在数字域实现的。

102变采样率数字滤波采样率转换通常分为“抽取(Decimation)”和“插值(Interpolation)”。抽取是降低采样率(所以在英文中也称为Down-sampling)，以去掉多余数据样本的过程，而插值则是提高采样率(Up-sampling)以增加数据样本的过程。变采样率是一个非线性问题，它涉及相当复杂的数学推导，通常不在本科阶段讨论它的设计问题。本书主要介绍抽取和插值的一般概念，一方面可以加深对采样定理的理解，同时可以帮助读者了解更多的应用问题，开扩思路。

103内插过程的数字滤波信号的整数倍内插整数倍内插是在已知的相邻两个采样点之间插入L-1个样本点。由于这些样本点并非已知的值，所以关键问题是如何求出这L-1个样本值。前提是要直接在数字域完成这个插值。用数字方法实现整数内插的步骤如下：在已知采样序列x(n1T1)的相邻两个样点之间等间隔插入L-1个零值样本点，得到输出y，则有104内插过程的数字滤波示例式中L为大于1的整数，称为插入因子。按上式y(n)是取值跳动的序列，让它经H(z)进行低通数字滤波，得到教平滑的内插结果。这种内插方案的框图如图9.7.1所示。各点的波形和频谱见图9.7.2。105内插中的数字滤波示例9.7.1Nx=10;L=3;fsx=2;T=1/fsx;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.0019,0.0203,0.0772,0.1589,0.2,0.1589,0.0772,0.0203,0.0019,0];%用例733的滤波器系数y1=filter(b,1,y);%进行滤波106内插中的数字滤波示例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频率位置序列%以下绘图语句，得出各组曲线如下图。107变采样率数字滤波108内插中的数字滤波示例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倍的冲激频谱。

109抽取过程的数字滤波信号的整数倍抽取整数倍抽取是在已知序列x中每隔M-1个采样点取出一个样本点，组成新的序列y。因此y序列中的所有值都是已知的，可以直接写出下列表达式

也就是说，输入序列x中下标为M的整倍数的点都被取入y序列，而其它的

样本都被舍去不用了。用MATLAB语言来表达，为：Ny=ceil(Nx/M);nx=0:Nx-1; y=x(0:M:Nx);ny=0:Ny-1;

110抽取过程的数字滤波 抽取的数学关系看起来较插入简单。但由于采样率的降低，造成的频率泄漏可能变得相当严重。要解决好这个问题，在进行抽取之前，必须对原来的序列进行预先数字滤波，称为‘抗泄漏滤波’。它的框图可表为图9.7.3

111抽取中的数字滤波算例9.7.2与框图对应的程序为hc972，先用滤波器滤波得到x1，再抽取得到y1。故其核心语句为：nx=0:Nx-1;x=exp(-alpha*nx*T);%给出x序列b=[0,0.0019,0.0203,0.0772,0.1589,0.2,0.1589,0.0772,0.0203,0.0019,0];%例733的滤波器系数x1=filter(b,1,x);%进行滤波，得出x1Ny=ceil(Nx/M);ny=0:Ny-1;y1=x1(1:M:Nx);%实行抽取，得出y1序列以下分别画出x,x1,y1,y2的频谱见下图。112变采样率数字滤波113抽取中的数字滤波算例9.7.2从图中左边一列看出抽取过程的时域波形，指数下降的输入信号x经过低通滤波成为波形x1，然后经过抽取而降低了采样率,成为y1.最下面的图表示不经过滤波直接抽取的结果。再从右边的频谱图上可以看出滤波的必要性。直接抽取波形的频谱abs(Y2)在奈奎斯特频率处（也就是频谱图的左右边界±π处）有很高的幅度，这意味着存在着很大的频率泄漏。而经过滤波再抽取的信号频谱abs(Y1)在该处的幅度就小得多，如果用一个好一点的滤波器，结果会更加理想。

114分数倍变采样率数字滤波分数倍变采样率和相应的MATLAB函数如果变换前后的采样频率的比满足两个整数之比L/M，则不难想象，可以通过先作L倍内插后作M倍抽取方法来实现，也可以先抽取后内插。其滤波器的设计当然有不同的要求和方法。MATLAB信号处理工具箱提供了一个变采样率函数resample，用它可以把变采样率的全过程在一步中完成。这个变采样率函数的典型调用格式为y=resample(x,L,M)。其中L为采样率提高的倍数，M为采样率降低的倍数，L和M必须取整数。

115变采样率数字滤波的应用整数倍抽取和内插的应用举例在数字电话系统中，采样频率取为F＝8kHz，希望传输尽量接近4kHz的音频带宽。但送话器发出的信号x（t）的带宽比4kHz大很多。因此，在A／D变换之前要对其进行模拟预滤波，以防止采样后发生频率混叠失真。这就要求该滤波器在4kHz处幅频特性很低，而为了使能传送的信号频带尽量宽，又要求该滤波器的通带接近于4kHz(例如3.8kHz)，给过渡带宽留的裕度只有几百Hz，这样的要求对用模拟低通滤波器是很难做到的。

116变采样率数字滤波的应用为了降低对模拟预滤波器的技术要求，新方案先用较高的采样率进行采样，比如取采样率Fsx＝16kHz，这样，模拟预滤波器的过渡带可以从3.8kHz到8kHz，这就容易设计了。经过A/D后，再经M=2倍抽取，把采样率降低至8kHz。此时问题的难点变成了抽取时的滤波器设计技术。因为这是是在数字域的滤波，用FIR结构的数字滤波器不难设计成线形相位和过渡带很窄的特性。这种方案最终可有效地利用通信带宽，并增加信号数据传送量。1179.8稀疏天线阵列设计

本书讨论的信号处理都是以时间作为自变量的。而在本例中自变量却是一维的空间位置。从本例可以看到，利用数学模型的相似性，时域的信号处理理论和方法可以推广到空域的信号处理中。

可以发现，等间隔分布的线性相位天线阵元所产生的远场方向图和FIR滤波器的频率响应在数学模型上非常相似。这就使我们可以根据方向图的要求来设计天线阵列。在实际中往往需要去掉其中某些阵元，称为稀疏阵列。118稀疏天线阵列设计图9.8.1所示为N+1个等间隔配置的方向性相同的阵元。阵元间距为d，因此它们在基线上排列的坐标为x(n)=n×d，0≤