数字信号处理MATLAB实现及综合应用

第7章数字信号处理MATLAB实现及综合应用7.1数字信号处理MATLAB实现7.2在双音拨号系统中的应用7.3在wav信号分析方面的应用7.4在自适应滤波器方面的应用随着信息学科和计算机学科的飞速发展，数字信号处理的重要性日益显著，数字信号的处理已成为电子信息学科中的工程技术人员必不可少的知识。由于数字信号处理的概念比较抽象，再加上其数值计算又比较繁琐，因此，国外很早就开始把MATLAB用于数字信号处理的教学过程，并取得了很好的经验，现在MATLAB已经成为解决数字信号处理问题的公认的标准软件。为了使读者更好地掌握数字信号处理的基本概念和基本理论，本章首先在7.1节给出前面各章的MATLAB典型例题。通过这些例题，帮助读者弄清一些重要的概念，加深对课程基本概念和基本原理的理解，学会利用计算机软件解决在理论学习中不易解决的问题。然后在7.2节~7.4节分别给出了数字信号处理技术在通信、信号分析、自适应滤波等方面的实际应用，并给出了MATLAB仿真程序和结果分析，这些内容可以扩宽读者数字信号处理技术应用的知识面，促进理论与实际相结合。本章和下一章上机实验都是课程学习的重要实践环节，它不仅能帮助同学们理解和掌握课程的重要概念和基本内容，而且也有益于读者深入学习和应用信号处理知识来解决实际问题。7.1数字信号处理MATLAB实现7.1.1MATLAB简介7.1.2离散时间信号与系统分析应用举例7.1.3离散傅里叶变换应用举例7.1.4IIR滤波器设计常用函数及举例7.1.5FIR滤波器设计常用函数及举例7.1.1MATLAB简介MATLAB是由美国MathWorks公司推出的软件产品，“MATLAB”是MatrixLaboratory（矩阵实验室）的缩写。MATLAB是一个完整的并可扩展的计算机环境，是一种进行科学和工程计算的交互式程序语言。它的基本数据单元是不需要指定维数的矩阵，解决同样的数值计算问题，使用MATLAB效率要提高许多倍。MATLAB采用开放式的环境，你可以读到它的原码、了解它的算法，并能改变当前的函数或增添自己编写的函数。目前，在国内外大学和研究机构中，MATLAB是一种非常流行的计算机语言，许多重要的学术刊物上发表的论文均是用MATLAB来分析计算和绘制各种图形。它还是一种有力的教学工具，在大学的线性代数、自动控制理论、数字信号处理、动态系统仿真等课程的教学中，已成为标准的教学工具。1．变量2．基本操作符3．输入、输出语句4．库函数5．绘图操作表7-1 MATLAB的常用绘图命令及功能实际的图形显示有多种情况，下面分别加以说明。（1）需要用多个图形窗显示多个图形时，应在图形显示语句之前加上图形窗指令figure；（2）在一个图形窗中显示多个图形时，应在图形显示语句之前调用图形窗分块函数subplot()；（3）在一个图形坐标系中显示多条曲线时，要在plot()语句前先执行指令holdon用以保持当前曲线，显示完成之后，再执行指令holdoff；（4）需要选择坐标轴尺寸时，在plot()语句之后调用axis()函数。信号处理工具箱覆盖了经典信号处理的大多数内容，是一个非常优秀的算法研究、辅助设计工具。它提供了大量的工具函数、滤波器分析和设计工具，常用的函数有以下几种。（1）设计IIR滤波器相关函数buttord、butter、cheblord、cheb2ord、ellipord、cheby1、cheby2等。6．用于数字信号处理的工具箱——SIGNAL（2）设计FIR滤波器相关函数fir1、fir2、firls、remezord、remez等；窗函数：boxcar、hanning、hamming、blackman、kaiser。（3）傅里叶变换相关函数FFT（一维快速傅里叶变换）、IFFT（一维快速傅里叶逆变换）、FFT2（二维快速傅里叶变换）、abs（绝对值或模）、angle（相角）等。（4）频率变换相关函数freqz、freqs、lp2lp、lp2hp、lp2bp、lp2bs。（5）频率响应相关函数freqz、freqs。7.1.2离散时间信号与系统分析应用举例1．zplane(b,a)函数或pzplotz(b,a)函数2．freqz函数3．filter函数图7-1零极点图图7-2频率响应图图7-3仿真波形7.1.3离散傅里叶变换应用举例在MATLAB中，可用函数dft来实现离散傅里叶变换运算。而且，在实用中，MATLAB还提供了内部函数来快速进行离散傅里叶变换运算。它采用了优化的算法，并且程序是用目的码编写的，使得它有极高的运算速度，命名为fft，反变换是ifft。调用方式为y=fft(x)或y=fft(x,N)。图7-44点DTFT和DFT图图7-5DFT图（N=8，N=16）图7-6x(n)的前10点数据对应的x(n)、X(ejw)、X(k)

图7-7x(n)补零至100点对应的x(n)、X(ejw)、X(k)

图7-9128点x(n)对应的x(n)、X(ejw)、X(k)7.1.4IIR滤波器设计常用函数及举例1．最小阶数选择函数[n,wn]=buttord/cheb1ord/cheb2ord/ellipord(wp,ws,Rp,Rs,'s')1．最小阶数选择函数[z,p,k]=buttap(n)/cheb1ap(n,Rp)/cheb2ap(n,Rs)/ellipap(n,Rp,Rs)2．模拟低通原型函数（1）低通到低通的频率变换[b1,a1]=lp2lp(b,a,w0)，其中，w0为低通滤波器的截止频率（rad/s）。（2）低通到高通的频率变换[b1,a1]=lp2hp(b,a,w0)，其中，w0为高通滤波器的截止频率（rad/s）。3．频率变换函数（3）低通到带通的频率变换[b1,a1]=lp2bp(b,a,w0,Bw)，其中，w0为带通滤波器的中心频率，Bw为带通滤波器的带宽。当滤波器通带的下截止频率为w1，上截止频率为w2时，w0=sqrt(w1*w2)，Bw=w2-w1。（4）低通到带阻的频率变换[b1,a1]=lp2bs(b,a,w0,Bw)，其中，w0为带阻滤波器的中心频率，Bw为带阻滤波器的带宽。当滤波器通带的下截止频率为w1，上截止频率为w2时，w0=sqrt(w1*w2)，Bw=w2-w1。MATLAB信号处理工具箱还提供了模拟滤波器设计的完全工具函数：butter、cheby1、cheby2、ellip、besself。用户只需一次调用就可自动完成全部设计过程，编程十分简单。这些工具函数既适用于模拟滤波器设计，也适用于数字滤波器。4．完全设计函数（1）巴特沃思滤波器：[b,a]=butter(n,wn,'ftype','s')其中，n为滤波器阶数；wn为滤波器截止频率；'s'为模拟滤波器，缺省时为数字滤波器。（2）切比雪夫滤波器：[b,a]=cheby1(n,Rp,wn,'ftype','s')或[b,a]=cheby2(n,Rs,wn,'ftype','s')（3）椭圆滤波器：[b,a]=ellip(n,Rp,Rs,wn,'ftype','s')[bz,az]=impinvar/bilinear(b,a,fs)5．脉冲响应和双线性变换函数其中，b、a分别为模拟滤波器的分子和分母多项式系数向量；fs为采样频率，单位为Hz，缺省时fs默认为1Hz；bz，az分别为数字滤波器分子和分母多项式系数向量。图7-10不同采样频率对所设计数字滤波器频率响应的影响图7-11ChebyshevⅠ型数字高通滤波器幅频特性图7-127阶ChebyshevⅡ型数字低通滤波器幅度和相位图7.1.5FIR滤波器设计常用函数及举例1．窗函数法设计FIR滤波器（2）计算截止频率为wc的理想低通滤波器的单位脉冲响应hd(n)，语句如下hd=ideal_lp(wc,N)（3）h(n)为h=hd.*win。（4）已知h(n)，为了得到频率响应的幅度、相位及群时延，需用语句如下[db,mag,pha,grd,w]=freqz_m(h,1)图7-13设计的实际低通滤波器结果（1）h(n)为h=real(ifft(H,N))。（2）MATLAB程序中，可以用freqz计算频率响应，但不能从中确定振幅响应。这是因为在MATLAB函数中，不存在与abs函数功能相当而能找到振幅的函数。但在MATLAB中，提供了4种扩展函数，可以计算FIR4种类型滤波器的振幅响应。●Hr_type1●Hr_type2●Hr_type3●Hr_type42．频率采样法设计FIR滤波器前面介绍的窗函数法设计FIR滤波器，完全是根据窗函数法的设计思路和设计流图一步一步进行的，这有助于掌握方法的原理。在实际应用中，还可以采用另外一种方法——直接利用MATLAB中的内部函数来进行FIR滤波器的设计。MATLAB信号处理工具箱提供了两个用窗函数法设计FIR滤波器的函数，见表7-2。值得一提的是，FIR滤波器的窗函数设计法不仅可以设计标准频率响应（低通、高通、带通及带阻）滤波器，还可设计任意频率响应的多带滤波器。3．直接采用MATLAB内部函数设计FIR滤波器表7-2 用窗函数法设计FIR滤波器的内部函数函数名函

数

功

能fir1用窗函数法设计具有标准频率响应的FIR滤波器fir2用窗函数法设计具有任意频率响应的多带FIR滤波器图7-14例7-10运行结果图（1）函数fir1例7-11设计一个低通滤波器，通带0～200Hz，阻带300～1000Hz，通带波动7%，阻带衰减60dB，采样频率为2000Hz。图7-15例7-11滤波器幅度和相位图（2）函数fir2图7-16被污染信号的频谱特性图7-17设计的多带滤波器幅频特性图7-18原始信号、污染信号及滤波后信号的比较7.2在双音拨号系统中的应用7.2.1电话的双音拨号系统7.2.2Goertzel算法7.2.3检测DTMF信号的DFT参数选择双音多频（DualToneMultiFrequency，DTMF）是用按键进行电话拨号的制式，它不单单用在电话中，还可以用于传输十进制数据的其他通信系统中。DTMF也广泛应用于电子邮件和银行系统。在这些系统中，用户可以从电话发送DTMF信号来选择菜单语音进行操作。图7-19双频拨号频率分配7.2.1电话的双音拨号系统7.2.2Goertzel算法图7-20计算X(k)的系统结构7.2.3检测DTMF信号的DFT参数选择DTMF信号的参数选择要考虑多方面的因素。首先电话数字化的采样频率已知为8kHz，实际上为了抑止语音干扰，检测系统除了规定的8个频率之外，还要检测它们的2次倍频处的DFT幅度。这些2倍频处幅度一般很小，如果基频和倍频分量同时都大，那就可能是外来声音的干扰。表7-4 各个拨号频率及N＝205时其对应的DFT序号8个基频Hz准确k值最近整数k值绝对误差2次谐波频率Hz准确k值最近整数k值绝对误差69717.861180.139139435.024350.02477019.531200.269154038.692390.30885221.833220.167170442.813430.18794124.113240.113188247.285470.285120930.981310.019241860.752610.248133634.235340.235267267.134670.134147737.848380.152295474.219740.219163341.846420.154326682.058820.0587.2.4DTMF信号检测的MATLAB仿真图7-21发送端产生的双频信号时域波形图7-22接收端双频信号在8个近似基频上的DFT幅度7.3在wav信号分析方面的应用7.3.1wav文件的一次性傅里叶 变换7.3.2wav文件的分段傅里叶分析本节将通过分析计算机中的wav文件来进一步讨论数字信号处理中的信号分析方法。首先介绍几个常用的函数：[x，fs，bits]＝waveread（‘filename’）——这是一个MATLAB中读取wav文件的数据函数。其中的x表示一长串的数据，一般是两列（立体声）；fs是该wav文件在采集时用的采样频率；bits是指在进行A/D转换时用的量化位长（一般是8bits或16bits）。[d]＝fft(w,N)——这是MATLAB中快速傅里叶变换（FFT）函数的一种输入输出形式。w是一列波形数据；N是指用多少点的FFT，一般选择2的整次方数（如16，128，1024等），因为这样就可以使用优化的蝶形算法；d是频域的输出。由于FFT（DFT）的对称性，且输入的是实数，FFT结果的复数序列是共轭反对称的，所以它们模是大小对称的，一般来说只用其一半的数据就可以了。sound（w，fs，bits）——和前面的waveread一样的参数表示，它将序列的数据通过声卡转化为声音。下面先介绍wav文件的一次性傅里叶变换。7.3.1wav文件的一次性傅里叶变换1．声波主要频率的分析图7-23信号的时域波形图7-24信号频域的幅值2．声音信号的重构前面介绍了对wav文件直接进行一次性傅里叶变换的分析，但存在两个问题：第一，对分析的wav文件直接一次性做FFT，点数太大，而在实际运用中DSP的FFT的点数是有限的，一般只能达到千点；第二，语音是分音节的，应该对它分段分析，所以下面引入分段处理数据，并连接每个段的频域结果。7.3.2