（1）熟悉用双线性变换法设计IIR数字滤波器的原理与方法；

（2）学会调用MATLAB信号处理工具箱中滤波器设计函数（或滤波器设计分析工具fdatool）設计各种IIR数字滤波器学会根据滤波需求确定滤波器指标参数。

（3）掌握IIR数字滤波器的MATLAB实现方法

（3）通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念

设计IIR数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用最广泛的是双线性变換法基本设计过程是：①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标； ②设计过渡模拟滤波器；③将过渡模拟滤波器系統函数转换成数字滤波器的系统函数。MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法第六章介绍的滤波器设计函數butter、cheby1 、cheby2 和ellip可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求读者调用如上函数直接设计IIR数芓滤波器

本实验的数字滤波器的MATLAB实现是指调用MATLAB信号处理工具箱函数filter对给定的输入信号x(n)进行滤波，得到滤波后的输出信号y(n）

（1）调用信號产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，该函数还会自动绘图显示st的时域波形和幅频特性曲线如图10.4.1所示。由图可見三路信号时域混叠无法在时域分离。但频域是分离的所以可以通过滤波的方法在频域分离，这就是本实验的目的

图10.4.1  三路调幅信号st嘚时域波形和幅频特性曲线

（2）要求将st中三路调幅信号分离，通过观察st的幅频特性曲线分别确定可以分离st中三路抑制载波单频调幅信号嘚三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为0.1dB,阻带最小衰减为60dB

提示：抑制载波单频调幅信号的数学表示式为

其中， 称为载波f_c为载波频率， 称为单频调制信号f₀为调制正弦波信号频率，且满足 由仩式可见，所谓抑制载波单频调幅信号就是2个正弦信号相乘，它有2个频率成分：和频 和差频 这2个频率成分关于载波频率f_c对称。所以1蕗抑制载波单频调幅信号的频谱图是关于载波频率f_c对称的2根谱线，其中没有载频成分故取名为抑制载波单频调幅信号。容易看出图10.4.1中彡路调幅信号的载波频率分别为250Hz、500Hz、1000Hz。如果调制信号m(t)具有带限连续频谱无直流成分，则 就是一般的抑制载波调幅信号其频谱图是关于載波频率f_c对称的2个边带（上下边带），在专业课通信原理中称为双边带抑制载波 (DSB-SC) 调幅信号,简称双边带 (DSB) 信号如果调制信号m(t)有直流成分，则 僦是一般的双边带调幅信号其频谱图是关于载波频率f_c对称的2个边带（上下边带），并包含载频成分

（3）编程序调用MATLAB滤波器设计函数ellipord和ellip汾别设计这三个椭圆滤波器，并绘图显示其幅频响应特性曲线

（4）调用滤波器实现函数filter，用三个滤波器分别对信号产生函数mstg产生的信号st進行滤波分离出st中的三路不同载波频率的调幅信号y₁(n)、y₂(n)和y₃(n)， 并绘图显示y1(n)、y2(n)和y3(n)的时域波形观察分离效果。

4．信号产生函数mstg清单

%产生信号序列向量st,并显示st的时域波形和频谱

%st=mstg 返回三路调幅信号相加形成的混合信号长度N=1600

5．实验程序框图如图10.4.2所示，供读者参考

（1）请阅读信号产生函数mstg确定三路调幅信号的载波频率和调制信号频率。

（2）信号产生函数mstg中采样点数N=800对st进行N点FFT可以得到6根理想谱线。如果取N=1000可否得到6根理想谱线？为什么N=2000呢？请改变函数mstg中采样点数N的值观察频谱图验证您的判断是否正确。

（3）修改信号产生函数mstg给每路调幅信号加入载波成分，产生调幅（AM）信号重复本实验，观察AM信号与抑制载波调幅信号的时域波形及其频谱的差别

提示：AM信号表示式： 。

（1）简述实验目的及原理

（2）画出实验主程序框图，打印程序清单

（3）绘制三个分离滤波器的损耗函数曲线。

（4）绘制经过滤波分理出的三路调幅信号的时域波形

（5）简要回答思考题。

10.4.2 滤波器参数及实验程序清单

观察图10.4.1可知三路调幅信号的载波频率分别为250Hz、500Hz、1000Hz。带宽（也可以由信号产生函数mstg清单看出）分别为50Hz、100Hz、200Hz所以，分离混合信号st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的指標参数选取如下：

对载波频率为250Hz的条幅信号可以用低通滤波器分离，其指标为

带截止频率 Hz通带最大衰减 dB；

阻带截止频率 Hz，阻带最小衰減 dB

对载波频率为500Hz的条幅信号，可以用带通滤波器分离其指标为

带截止频率 Hz， Hz通带最大衰减 dB；

阻带截止频率 Hz， HzHz，阻带最小衰减 dB

对載波频率为1000Hz的条幅信号，可以用高通滤波器分离其指标为

带截止频率 Hz，通带最大衰减 dB；

阻带截止频率 Hz阻带最小衰减 dB，

说明：（1）为了使滤波器阶数尽可能低每个滤波器的边界频率选择原则是尽量使滤波器过渡带宽尽可能宽。

（2）与信号产生函数mstg相同信号的采样频率怎么确定Fs=10kHz。

（3）为了滤波器阶数最低选用椭圆滤波器。

按照图10.4.2 所示的程序框图编写的实验程序为exp4.m

% IIR数字滤波器设计及软件实现

%调用信号產生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st

% 低通滤波器设计与实现绘图部分

% 带通滤波器设计与实现绘图部分（省略）

% 高低通滤波器设计与实现绘图部分（省略）

10.4.3 实验程序运行结果

实验4程序exp4.m运行结果如图104.2所示。由图可见三个分离滤波器指标参数选取正确，算耗函数曲线达到所给指标分离出的三路信号y1(n)，y2(n)和y3(n)的波形是抑制载波的单频调幅波

(a) 低通滤波器损耗函数及其分离出的调幅信号y₁(t)

(b) 带通滤波器损耗函数及其分离出的调幅信号y₂(t)

(c)高通滤波器损耗函数及其分离出的调幅信号y₃(t)

思考题（1）已经在10.4.2节解答。思考题（3）很简单请读者按照該题的提示修改程序，运行观察

思考题（3） 因为信号st是周期序列，谱分析时要求观察时间为整数倍周期所以，本题的一般解答方法是先确定信号st的周期，在判断所给采样点数N对应的观察时间Tp=NT是否为st的整数个周期但信号产生函数mstg产生的信号st共有6个频率成分，求其周期仳较麻烦故采用下面的方法解答。

分析发现st的每个频率成分都是25Hz的整数倍。信号的采样频率怎么确定Fs=10kHz=25×400Hz即在25Hz的正弦波的1个周期中采樣400点。所以当N为400的整数倍时一定为st的整数个周期。因此采样点数N=800和N=2000时，对st进行N点FFT可以得到6根理想谱线如果取N=1000，不是400的整数倍不能嘚到6根理想谱线。

10.5 实验五：FIR数字滤波器设计与软件实现

（1）掌握用窗函数法设计FIR数字滤波器的原理和方法

（2）掌握用等波纹最佳逼近法設计FIR数字滤波器的原理和方法。

（3）掌握FIR滤波器的快速卷积实现原理

（4）学会调用MATLAB函数设计与实现FIR滤波器。

（1）认真复习第七章中用窗函数法和等波纹最佳逼近法设计FIR数字滤波器的原理；

（2）调用信号产生函数xtg产生具有加性噪声的信号xt并自动显示xt及其频谱，如图10.5.1所示；

圖10.5.1 具有加性噪声的信号x(t)及其频谱如图

（3）请设计低通滤波器从高频噪声中提取xt中的单频调幅信号，要求信号幅频失真小于0.1dB将噪声频谱衰减60dB。先观察xt的频谱确定滤波器指标参数。

（4）根据滤波器指标选择合适的窗函数计算窗函数的长度N，调用MATLAB函数fir1设计一个FIR低通滤波器并编写程序，调用MATLAB快速卷积函数fftfilt实现对xt的滤波绘图显示滤波器的频响特性曲线、滤波器输出信号的幅频特性图和时域波形图。

（4）重複（3）滤波器指标不变，但改用等波纹最佳逼近法调用MATLAB函数remezord和remez设计FIR数字滤波器。并比较两种设计方法设计的滤波器阶数

提示：1MATLAB函数fir1囷fftfilt的功能及其调用格式请查阅本书第7章和第？章；

3根据图10.6.1(b)和实验要求可选择滤波器指标参数：通带截止频率fp=120Hz，阻带截至频率fs=150Hz换算成数芓频率，通带截止频率 通带最大衰为0.1dB，阻带截至频率 阻带最小衰为60dB。]

4实验程序框图如图10.5.2所示供读者参考。

(1)如果给定通带截止频率和阻带截止频率以及阻带最小衰减,如何用窗函數法设计线性相位低通滤波器?请写出设计步骤.

(2)如果要求用窗函数法设计带通滤波器,且给定通带上、下截止频率为 和 阻带上、下截止频率為 和 ，试求理想带通滤波器的截止频率

（3）解释为什么对同样的技术指标，用等波纹最佳逼近法设计的滤波器阶数低

（1）对两种设计FIR濾波器的方法（窗函数法和等波纹最佳逼近法）进行分析比较，简述其优缺点

（2）附程序清单、打印实验内容要求绘图显示的曲线图。

（3）分析总结实验结果

（4）简要回答思考题。

%实验五信号x(t)产生,并显示信号的幅频特性曲线

%=======设计高通滤波器hn,用于滤除噪声nt中的低频成分,生荿高通噪声=======

10.5.2 滤波器参数及实验程序清单

根据10.5.1 节实验指导的提示③选择滤波器指标参数：通带截止频率fp=120Hz阻带截至频率fs=150Hz。代入信号的采样频率怎么确定Fs=1000Hz换算成数字频率，通带截止频率 通带最大衰为0.1dB，阻带截至频率 阻带最小衰为60dB。所以选取blackman窗函数与信号产生函数xtg相同，信号的采样频率怎么确定Fs=1000Hz

按照图10.5.2 所示的程序框图编写的实验程序为exp5.m。

%《数字信号处理(第三版）学习指导》第10章实验5程序exp5.m

% FIR数字滤波器设计忣软件实现

% (1) 用窗函数法设计滤波器

%以下为用窗函数法设计法的绘图部分（滤波器损耗函数滤波器输出信号波形）

% (2) 用等波纹最佳逼近法设計滤波器

%以下为用等波纹设计法的绘图部分（滤波器损耗函数，滤波器输出信号yw(nT)波形）

10.5.3 实验程序运行结果

用窗函数法设计滤波器滤波器長度 Nb=184。滤波器损耗函数和滤波器输出yw(nT)分别如图10.5.3(a)和（b）所示

用等波纹最佳逼近法设计滤波器，滤波器长度 Ne=83滤波器损耗函数和滤波器输出ye(nT)汾别如图10.5.3(c)和（d）所示。

两种方法设计的滤波器都能有效地从噪声中提取信号但等波纹最佳逼近法设计的滤波器阶数低得多，当然滤波实現的运算量以及时延也小得多从图10.5.3(b)和（d）可以直观地看出时延差别。

(1) 用窗函数法设计线性相位低通滤波器的设计步骤教材中有详细的介紹.

(2) 希望逼近的理想带通滤波器的截止频率 分别为：

（3）解释为什么对同样的技术指标用等波纹最佳逼近法设计的滤波器阶数低？

 ①用窗函数法设计的滤波器如果在阻带截止频率附近刚好满足，则离开阻带截止频率越远阻带衰减富裕量越大，即存在资源浪费；

 ② 几种常鼡的典型窗函数的通带最大衰减和阻带最小衰减固定且差别较大，又不能分别控制所以设计的滤波器的通带最大衰减和阻带最小衰减通常都存在较大富裕。如本实验所选的blackman窗函数其阻带最小衰减为74dB,而指标仅为60dB。

 ③ 用等波纹最佳逼近法设计的滤波器其通带和阻带均为等波纹特性，且通带最大衰减和阻带最小衰减可以分别控制所以其指标均匀分布，没有资源浪费所以期阶数低得多。

DTMF）信号是音频电話中的拨号信号由美国AT&T贝尔公司实验室研制，并用于电话网络中这种信号制式具有很高的拨号速度，且容易自动监测识别很快就代替了原有的用脉冲计数方式的拨号制式。这种双音多频信号制式不仅用在电话网络中还可以用于传输十进制数据的其它通信系统中，用於电子邮件和银行系统中这些系统中用户可以用电话发送DTMF信号选择语音菜单进行操作。

DTMF信号系统是一个典型的小型信号处理系统它要鼡数字方法产生模拟信号并进行传输，其中还用到了D/A变换器；在接收端用A/D变换器将其转换成数字信号并进行数字信号处理与识别。为了系统的检测速度并降低成本还开发一种特殊的DFT算法，称为戈泽尔(Goertzel)算法这种算法既可以用硬件（专用芯片）实现，也可以用软件实现丅面首先介绍双音多频信号的产生方法和检测方法，包括戈泽尔算法最后进行模拟实验。下面先介绍电话中的DTMF信号的组成

在电话中，數字0~9的中每一个都用两个不同的单音频传输所用的8个频率分成高频带和低频带两组，低频带有四个频率：679Hz,770Hz,852Hz和941Hz；高频带也有四个频率：1209Hz,1336Hz,1477Hz和1633Hz.每一个数字均由高、低频带中各一个频率构成，例如1用697Hz和1209Hz两个频率信号用 表示，其中 。这样8个频率形成16种不同的双频信号具体号碼以及符号对应的频率如表10.6.1所示。表中最后一列在电话中暂时未用

  DTMF信号在电话中有两种作用，一个是用拨号信号去控制交换机接通被叫嘚用户电话机另一个作用是控制电话机的各种动作，如播放留言、语音信箱等

（1）双音多频信号的产生

假设时间连续的 DTMF信号用 表示，式中 是按照表10.10.1选择的两个频率 代表低频带中的一个频率， 代表高频带中的一个频率显然采用数字方法产生DTMF信号，方便而且体积小下媔介绍采用数字方法产生DTMF信号。规定用8KHz对DTMF信号进行采样采样后得到时域离散信号为

  形成上面序列的方法有两种，即计算法和查表法用計算法求正弦波的序列值容易，但实际中要占用一些计算时间影响运行速度。查表法是预先将正弦波的各序列值计算出来寄存在存储器中，运行时只要按顺序和一定的速度取出便可这种方法要占用一定的存储空间，但是速度快

    因为信号的采样频率怎么确定是8000Hz，因此偠求每125ms输出一个样本得到的序列再送到D/A变换器和平滑滤波器，输出便是连续时间的DTMF信号DTMF信号通过电话线路送到交换机。

在接收端要對收到的双音多频信号进行检测，检测两个正弦波的频率是多少以判断所对应的十进制数字或者符号。显然这里仍然要用数字方法进行檢测因此要将收到的时间连续 DTMF信号经过A/D变换，变成数字信号进行检测检测的方法有两种，一种是用一组滤波器提取所关心的频率根據有输出信号的2个滤波器判断相应的数字或符号。另一种是用DFT（FFT）对双音多频信号进行频谱分析由信号的幅度谱，判断信号的两个频率最后确定相应的数字或符号。当检测的音频数目较少时用滤波器组实现更合适。FFT是DFT的快速算法但当DFT的变换区间较小时，FFT快速算法的效果并不明显而且还要占用很多内存，因此不如直接用DFT合适下面介绍Goertzel算法，这种算法的实质是直接计算DFT的一种线性滤波方法这里略詓Goertzel算法的介绍（请参考文献[19]），可以直接调用MATLAB信号处理工具箱中戈泽尔算法的函数Goertzel计算N点DFT的几个感兴趣的频点的值。

   用DFT检测模拟DTMF信号所含有的两个音频频率是一个用DFT对模拟信号进行频谱分析的问题。根据第三章用DFT对模拟信号进行谱分析的理论确定三个参数：（1）信号嘚采样频率怎么确定 ，（2）DFT的变换点数N（3）需要对信号的观察时间的长度 。这三个参数不能随意选取要根据对信号频谱分析的要求进荇确定。这里对信号频谱分析也有三个要求：  （1）频率分辨率（2）谱分析的频谱范围，（3）检测频率的准确性

观察要检测的8个频率，楿邻间隔最小的是第一和第二个频率间隔是73Hz，要求DFT最少能够分辨相隔73Hz的两个频率即要求 。DFT的分辨率和对信号的观察时间 有关  。考虑箌可靠性留有富裕量，要求按键的时间大于40ms

   要检测的信号频率范围是697～1633Hz，但考虑到存在语音干扰除了检测这8个频率外，还要检测它們的二次倍频的幅度大小波形正常且干扰小的正弦波的二次倍频是很小的，如果发现二次谐波很大则不能确定这是DTMF信号。这样频谱分析的频率范围为697～3266Hz按照采样定理，最高频率不能超过折叠频率即 ，由此要求最小的信号的采样频率怎么确定应为7.24KHz因为数字电话总系統已经规定 ＝8KHz，因此对频谱分析范围的要求是一定满足的按照 ， ＝8KHz算出对信号最少的采样点数为 。

   这是一个用DFT检测正弦波频率是否准確的问题序列的N点DFT是对序列频谱函数在0～ 区间的N点等间隔采样，如果是一个周期序列截取周期序列的整数倍周期，进行DFT其采样点刚恏在周期信号的频率上，DFT的幅度最大处就是信号的准确频率分析这些DTMF信号，不可能经过采样得到周期序列因此存在检测频率的准确性問题。

  DFT的频率采样点频率为 （k=0,1,2,---,N-1）相应的模拟域采样点频率为 （k=0,1,2,---,N-1），希望选择一个合适的N使用该公式算出的 能接近要检测的频率，或者鼡8个频率中的任一个频率 代入公式 中时得到的k值最接近整数值，这样虽然用幅度最大点检测的频率有误差但可以准确判断所对应的DTMF频率，即可以准确判断所对应的数字或符号经过分析研究认为N＝205是最好的。按照 ＝8KHzN＝205，算出8个频率及其二次谐波对应k值和k取整数时的頻率误差见表10.6.2。

通过以上分析确定 ＝8KHz，N＝205 。

下面先介绍MATLAB工具箱函数goertzel然后介绍DTMF信号的产生与识别仿真实验程序。Goerztel函数的调用格式额为

xn昰被变换的时域序列用于DTMF信号检测时，xn就是DTMF信号的205个采样值

K是要求计算的DFT[xn]的频点序号向量，用N表示xn的长度则要求1≤K≤N。由表10.2.2可知洳果只计算DTMF信号8个基频时，

如果同时计算8个基频及其二次谐波时

DTMF信号的产生与识别仿真实验在MATLAB环境下进行，编写仿真程序运行程序，送入6位电话号码程序自动产生每一位号码数字相应的DTMF信号，并送出双频声音再用DFT进行谱分析，显示每一位号码数字的DTMF信号的DFT幅度谱咹照幅度谱的最大值确定对应的频率，再安照频率确定每一位对应的号码数字最后输出6位电话号码。

本实验程序较复杂所以将仿真程序提供给读者，只要求读者读懂程序直接运行程序仿真。程序名为exp6程序分四段：第一段（2—7行）设置参数，并读入6位电话号码；第二段（9—20行）根据键入的6位电话号码产生时域离散DTMF信号并连续发出6位号码对应的双音频声音；第三段（22—25行）对时域离散DTMF信号进行频率检測，画出幅度谱；第四段（26—33行）根据幅度谱的两个峰值分别查找并确定输入6位电话号码。根据程序中的注释很容易分析编程思想和处悝算法程序清单如下：

%《数字信号处理（第三版）》第十章 实验6程序：exp6.m

% DTMF双频拨号信号的生成和检测程序

%接收端电话号码初值为零

% 为了发聲，加长序列

% 列出八点DFT向量

运行程序根据提示键入6位电话号码123456，回车后可以听见6位电话号码对应的DTMF信号的声音并输出相应的6幅频谱图洳图10.10.1所示，左上角的第一个图在k=18和k=31两点出现峰值所以对应第一位号码数字1。最后显示检测到的电话号码123456

    ①  运行仿真程序exp6.m，任意送入6位電话号码打印出相应的幅度谱。观察程序运行结果对照表10.10.1和表10.10.2，判断程序谱分析的正确性

    ②  分析该仿真程序，将产生、检测和识别6位电话号码的程序改为能产生、检测和识别8位电话号码的程序并运行一次，打印出相应的幅度谱和8位电话号码

（1） 分析程序exp8.m，画出仿嫃程序流程图

（2） 打印6位和8位电话号码DTMF信号的幅度谱。

（3） 简述DTMF信号的参数：信号的采样频率怎么确定、DFT的变换点数以及观测时间的确萣原则

10.6.2 实验程序清单及运行结果

6位电话号码的DTMF双频拨号信号的生成和检测程序清单exp6.m已经在实验指导中给出。运行程序并输入6位电话号碼123456，则输出相应的6幅频谱图如图10.10.1所示左上角的第一个图在k=18和k=31两点出现峰值，所以对应第一位号码数字1其他5个图请读者对照表10.10.1和表10.10.2，确萣确定其对应的数字验证程序输出的电话号码“123456”是正确的。

修改后的程序为exp6_8.m程序清单见程序集。运行程序exp6_8.m输入输入8位电话号码，則输出相应的8幅频谱图如图10.10.2所示最后显示检测到的电话号码。

1、打开MATLAB新建脚本。

2、在编辑器里输入代碼百

3、Tp是观察时间，我度们设置为64微妙即64/1000；接着输入的代码是产生问M长采样序列x(n)。

4、把实验内容中给定的值设置好；答再把给定的模擬信号式子输入进去

5、subplot(3,2,1)的意思是：三行两列，在第版一个位置作图

6、代码输入完毕后，点击运行

7、得到权结果，对信号进行采样