当前位置:首页 >> 计算机软件及应用 >>

数字信号处理上机实验及答案(第三版,第十章)


上机实验

第十章

上机实验

数字信号处理是一门理论和实际密切结合的课程,为深入掌握课程内容,最好在学习理 论的同时,做习题和上机实验。上机实验不仅可以帮助读者深入的理解和消化基本理论,而 且能锻炼初学者的独立解决问题的能力。 本章在第二版的基础上编写了六个实验, 前五个实 验属基础理论实验,第六个属应用综合实验。 实验一 系统响应及系统稳定性。 实验二 时域采样与频域采样。 实验三 用 FFT 对信号作频谱分析。 实验四 IIR 数字滤波器设计及软件实现。 实验五 FIR 数字滤波器设计与软件实现 实验六 应用实验——数字信号处理在双音多频拨号系统中的应用 任课教师根据教学进度,安排学生上机进行实验。建议自学的读者在学习完第一章后作 实验一;在学习完第三、四章后作实验二和实验三;实验四 IIR 数字滤波器设计及软件实现 在。学习完第六章进行;实验五在学习完第七章后进行。实验六综合实验在学习完第七章或 者再后些进行;实验六为综合实验,在学习完本课程后再进行。

10.1

实验一: 实验一 系统响应及系统稳定性

1.实验目的 (1)掌握 求系统响应的方法。 (2)掌握时域离散系统的时域特性。 (3)分析、观察及检验系统的稳定性。 2.实验原理与方法 在时域中,描写系统特性的方法是差分方程和单位脉冲响应,在频域可以用系统函数 描述系统特性。 已知输入信号可以由差分方程、 单位脉冲响应或系统函数求出系统对于该输 入信号的响应,本实验仅在时域求解。在计算机上适合用递推法求差分方程的解,最简单的 方法是采用 MATLAB 语言的工具箱函数 filter 函数。也可以用 MATLAB 语言的工具箱函数 conv 函数计算输入信号和系统的单位脉冲响应的线性卷积,求出系统的响应。 系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。重点分析实验系统 的稳定性,包括观察系统的暂态响应和稳定响应。 系统的稳定性是指对任意有界的输入信号,系统都能得到有界的系统响应。或者系统 的单位脉冲响应满足绝对可和的条件。系统的稳定性由其差分方程的系数决定。 实际中检查系统是否稳定,不可能检查系统对所有有界的输入信号,输出是否都是有 界输出, 或者检查系统的单位脉冲响应满足绝对可和的条件。 可行的方法是在系统的输入端 加入单位阶跃序列,如果系统的输出趋近一个常数(包括零) ,就可以断定系统是稳定的[19]。 系统的稳态输出是指当 n → ∞ 时,系统的输出。如果系统稳定,信号加入系统后,系统输 出的开始一段称为暂态效应,随 n 的加大,幅度趋于稳定,达到稳态输出。 注意在以下实验中均假设系统的初始状态为零。 3.实验内容及步骤 (1)编制程序,包括产生输入信号、单位脉冲响应序列的子程序,用 filter 函数或 conv 函数求解系统输出响应的主程序。程序中要有绘制信号波形的功能。 (2)给定一个低通滤波器的差分方程为

y (n) = 0.05 x(n) + 0.05 x(n ? 1) + 0.9 y (n ? 1)
1

上机实验

输入信号 x1 (n) = R8 ( n)

x 2 ( n) = u ( n)
a) 分别求出系统对 x1 (n) = R8 ( n) 和 x 2 ( n) = u ( n) 的响应序列,并画出其波形。 b) 求出系统的单位冲响应,画出其波形。 (3)给定系统的单位脉冲响应为

h1 (n) = R10 (n)

h2 (n) = δ (n) + 2.5δ (n ? 1) + 2.5δ (n ? 2) + δ (n ? 3)
并画出波形。 用线性卷积法分别求系统 h1(n)和 h2(n)对 x1 (n) = R8 ( n) 的输出响应, (4)给定一谐振器的差分方程为

y (n) = 1.8237 y (n ? 1) ? 0.9801 y (n ? 2) + b0 x(n) ? b0 x(n ? 2)
令 b0 = 1 / 100.49 ,谐振器的谐振频率为 0.4rad。 a) 用实验方法检查系统是否稳定。输入信号为 u (n) 时,画出系统输出波形。 b) 给定输入信号为

x(n) = sin(0.014n) + sin(0.4n)
求出系统的输出响应,并画出其波形。 4.思考题 (1) 如果输入信号为无限长序列,系统的单位脉冲响应是有限长序列,可否用线性卷积 法求系统的响应? 如何求? (2)如果信号经过低通滤波器,把信号的高频分量滤掉,时域信号会有何变化,用前面 第一个实验结果进行分析说明。 5.实验报告要求 (1)简述在时域求系统响应的方法。 (2)简述通过实验判断系统稳定性的方法。分析上面第三个实验的稳定输出的波形。 (3)对各实验所得结果进行简单分析和解释。 (4)简要回答思考题。 (5)打印程序清单和要求的各信号波形。

10.1.2 实验程序清单
%实验 1:系统响应及系统稳定性 实验 : close all;clear all %======内容 1:调用 filter 解差分方程,由系统对 u(n)的响应判断稳定性 解差分方程, 的响应判断稳定性====== 内容 : 的响应判断稳定性 A=[1,-0.9];B=[0.05,0.05]; %系统差分方程系数向量 B 和 A 系统差分方程系数向量 x1n=[1 1 1 1 1 1 1 1 zeros(1,50)]; %产生信号 x1(n)=R8(n) 产生信号 x2n=ones(1,128); %产生信号 x2(n)=u(n) 产生信号

2

上机实验

hn=impz(B,A,58); %求系统单位脉冲响应 h(n) 求系统单位脉冲响应 subplot(2,2,1);y='h(n)';tstem(hn,y); %调用函数 tstem 绘图 调用函数 title('(a) 系统单位脉冲响应 h(n)');box on y1n=filter(B,A,x1n); %求系统对 x1(n)的响应 y1(n) 求系统对 的响应 subplot(2,2,2);y='y1(n)';tstem(y1n,y); title('(b) 系统对 R8(n)的响应 y1(n)');box on 的响应 y2n=filter(B,A,x2n); %求系统对 x2(n)的响应 y2(n) 求系统对 的响应 subplot(2,2,4);y='y2(n)';tstem(y2n,y); title('(c) 系统对 u(n)的响应 y2(n)');box on 的响应 %===内容 2:调用 conv 函数计算卷积 函数计算卷积============================ 内容 : x1n=[1 1 1 1 1 1 1 1 ]; %产生信号 x1(n)=R8(n) 产生信号 h1n=[ones(1,10) zeros(1,10)]; h2n=[1 2.5 2.5 1 zeros(1,10)]; y21n=conv(h1n,x1n); y22n=conv(h2n,x1n); figure(2) subplot(2,2,1);y='h1(n)';tstem(h1n,y); %调用函数 tstem 绘图 调用函数 title('(d) 系统单位脉冲响应 h1(n)');box on subplot(2,2,2);y='y21(n)';tstem(y21n,y); title('(e) h1(n)与 R8(n)的卷积 y21(n)');box on 与 的卷积 subplot(2,2,3);y='h2(n)';tstem(h2n,y); %调用函数 tstem 绘图 调用函数 title('(f) 系统单位脉冲响应 h2(n)');box on subplot(2,2,4);y='y22(n)';tstem(y22n,y); title('(g) h2(n)与 R8(n)的卷积 y22(n)');box on 与 的卷积 %=========内容 3:谐振器分析 内容 :谐振器分析======================== un=ones(1,256); %产生信号 u(n) 产生信号 n=0:255; xsin=sin(0.014*n)+sin(0.4*n); %产生正弦信号 产生正弦信号 A=[1,-1.8237,0.9801];B=[1/100.49,0,-1/100.49]; %系统差分方程系数向量 B 和 A 系统差分方程系数向量 y31n=filter(B,A,un); %谐振器对 u(n)的响应 y31(n) 谐振器对 的响应 y32n=filter(B,A,xsin); %谐振器对 u(n)的响应 y31(n) 谐振器对 的响应 figure(3) subplot(2,1,1);y='y31(n)';tstem(y31n,y); title('(h) 谐振器对 u(n)的响应 y31(n)');box on 的响应 subplot(2,1,2);y='y32(n)';tstem(y32n,y); title('(i) 谐振器对正弦信号的响应 y32(n)');box on

10.1.3 实验程序运行结果及分析讨论 实验程序运行结果及分析讨论
程序运行结果如图 10.1.1 所示。
实验内容(2)系统的单位冲响应、系统对 x1 (n) = R8 ( n) 和 x 2 ( n) = u ( n) 的响应序列 分别如图(a)、(b)和(c)所示; 实验内容(3)系统 h1(n)和 h2(n)对 x1 (n) = R8 ( n) 的输出响应分别如图(e)和(g)所示;

3

上机实验

实验内容(4)系统对 u (n) 和 x ( n) = sin(0.014n) + sin(0.4n) 的响应序列分别如图(h) 和(i)所示。由图(h)可见,系统对 u (n) 的响应逐渐衰减到零,所以系统稳定。由图(i)可见, 系统对 x ( n) = sin(0.014n) + sin(0.4n) 的稳态响应近似为正弦序列 sin(0.4n) , 这一结论 验证了该系统的谐振频率是 0.4 rad。
(a) 系 系 系 系 系 系 系 系 h(n) 0.1 0.08 0.06 0.04 0.02 0 20 n 40 0 40 n (c) 系 系 谐 u(n)的 系 系 y2(n) 20 y1(n) h(n) 0.4 0.2 0.6 (b) 系 系 谐 R8(n)的 系 系 y1(n)

1 0.8 y2(n) 0.6 0.4 0.2 0 50 n 100

4

上机实验

(d) 系 系 系 系 系 系 系 系 h1(n) 1 y21(n) 0 10 15 n (f) 系 系 系 系 系 系 系 系 h2(n) 5 h1(n) 8 6 4 2 0 0

(e) h1(n)与 R8(n)的 的 的 y21(n)

0.5

0

10

20

n (g) h2(n)与 R8(n)的 的 的 y22(n) 8 6 y22(n) 4 2

3

2 h2(n) 1 0

0

5 n

10

0

0

5

10 n

15

20

(h) 谐 谐 谐 谐 u(n)的 系 系 y31(n) 0.04 y31(n) 0.02 0 -0.02 -0.04 0 50 100 150 200 250 n (i) 谐 谐 谐 谐 谐 谐 谐 谐 的 系 系 y32(n) 1 0.5 y32(n) 0 -0.5 0 50 100 n 150 200 250

图 10.1.1

10.1.4

简答思考题

(1) 如果输入信号为无限长序列,系统的单位脉冲响应是有限长序列,可否用线性卷积

5

上机实验

法求系统的响应。①对输入信号序列分段;②求单位脉冲响应 h(n)与各段的卷积;③将各段 卷积结果相加。具体实现方法有第三章介绍的重叠相加法和重叠保留法。 (2)如果信号经过低通滤波器,把信号的高频分量滤掉,时域信号的剧烈变化将被平滑, 由实验内容(1)结果图 10.1.1(a)、(b)和(c)可见,经过系统低通滤波使输入信号 δ ( n) 、

x1 (n) = R8 (n) 和 x 2 (n) = u (n) 的阶跃变化变得缓慢上升与下降。

10.2

实验二

时域采样与频域采样

10.2.1 实验指导
1. 实验目的 时域采样理论与频域采样理论是数字信号处理中的重要理论。要求掌握模拟信号采样 前后频谱的变化, 以及如何选择采样频率才能使采样后的信号不丢失信息; 要求掌握频率域 采样会引起时域周期化的概念,以及频率域采样定理及其对频域采样点数选择的指导作用。 2. 实验原理与方法 时域采样定理的要点是: a) 对模拟信号 x a (t ) 以间隔 T 进行时域等间隔理想采样,形成的采样信号的频谱

? X ( j?) 是原模拟信号频谱 X a ( j?) 以采样角频率 ? s ( ? s = 2π / T )为周
期进行周期延拓。公式为:


1 ? ? X a ( j?) = FT [ x a (t )] = T

n = ?∞

∑X

a

( j? ? jn? s )

b) 采样频率 ? s 必须大于等于模拟信号最高频率的两倍以上,才能使采样信号的 频谱不产生频谱混叠。 利用计算机计算上式并不方便,下面我们导出另外一个公式,以便用计算机上进行实验。

? 理想采样信号 x a (t ) 和模拟信号 x a (t ) 之间的关系为: ? x a (t ) = x a (t ) ∑ δ (t ? nT )
n = ?∞ ∞

对上式进行傅立叶变换,得到:
∞ ? X a ( j?) = ∫ [ x a (t ) ∑ δ (t ? nT )]e ? j?t dt ?∞ n = ?∞ ∞

=∑



n = ?∞






?∞

x a (t )δ (t ? nT )e ? j?t dt

在上式的积分号内只有当 t = nT 时,才有非零值,因此:

? X a ( j?) =

n = ?∞

∑x

a

(nT )e ? j?nT
6

上机实验

上式中,在数值上 x a (nT ) = x (n) ,再将 ω = ?T 代入,得到:

? X a ( j? ) =

n = ?∞

∑ x ( n) e



? jωn

上式的右边就是序列的傅立叶变换 X (e



) ,即

? X a ( j?) = X (e jω ) ω =?T
上式说明理想采样信号的傅立叶变换可用相应的采样序列的傅立叶变换得到, 只要将自变量 ω用 ?T 代替即可。 频域采样定理的要点是: a) 对信号 x(n)的频谱函数 X(e


)在[0,2π]上等间隔采样 N 点,得到

X N (k ) = X (e jω )

ω = 2π k N

, k = 0,1, 2,L , N ? 1

则 N 点 IDFT[ X N ( k ) ]得到的序列就是原序列 x(n)以 N 为周期进行周期延拓后的主值区 序列,公式为:

xN (n) = IDFT[ X N (k )]N = [ ∑ x(n + iN )]RN (n)
i =?∞



b) 由上式可知,频域采样点数 N 必须大于等于时域离散信号的长度 M(即 N≥M),才 能使时域不产生混叠,则 N 点 IDFT[ X N ( k ) ]得到的序列 xN ( n) 就是原序列 x(n),即

xN (n) =x(n)。如果 N>M, xN (n) 比原序列尾部多 N-M 个零点;如果 N<M,z 则 xN (n) =IDFT[ X N (k ) ]发生了时域混叠失真,而且 xN (n) 的长度 N 也比 x(n)的长度
M 短,因此。 xN ( n) 与 x(n)不相同。 在数字信号处理的应用中,只要涉及时域或者频域采样,都必须服从这两个采样理论的 要点。 对比上面叙述的时域采样原理和频域采样原理,得到一个有用的结论,这两个采样理论 具有对偶性: “时域采样频谱周期延拓,频域采样时域信号周期延拓” 。因此放在一起进行实 验。 3. 实验内容及步骤 (1)时域采样理论的验证。 给定模拟信号, x a (t ) = Ae
?αt

sin(? 0 t )u (t )

式中 A=444.128, =50 2 π, 0 =50 2 πrad/s, α ? 它的幅频特性曲线如图 10.2.1

7

上机实验

图 10.2.1

x a (t ) 的幅频特性曲线

现用 DFT(FFT)求该模拟信号的幅频特性,以验证时域采样理论。 安照 x a (t ) 的幅频特性曲线,选取三种采样频率,即 Fs =1kHz,300Hz,200Hz。观 测时间选 T p = 50ms 。 为使用 DFT,首先用下面公式产生时域离散信号,对三种采样频率,采样序列按顺 序用 x1 ( n) , x 2 ( n) , x3 ( n) 表示。

x(n) = x a (nT ) = Ae ?αnT sin(? 0 nT )u (nT )
因为采样频率不同,得到的 x1 ( n) , x 2 ( n) , x3 ( n) 的长度不同, 长度(点数)用 公式 N = T p × Fs 计算。选 FFT 的变换点数为 M=64,序列长度不够 64 的尾部加零。 X(k)=FFT[x(n)] , k=0,1,2,3,-----,M-1

式中 k 代表的频率为 ω k =

2π k。 M

要求: 编写实验程序,计算 x1 ( n) 、 x 2 ( n) 和 x3 ( n) 的幅度特性,并绘图显示。观察分 析频谱混叠失真。 (2)频域采样理论的验证。 )频域采样理论的验证。 给定信号如下:

? n + 1 0 ≤ n ≤ 13 ? x(n) = ?27 ? n 14 ≤ n ≤ 26 ? 0 其它 ?
编写程序分别对频谱函数 X (e jω ) = FT[ x ( n)] 在区间 [0,2π ] 上等间隔采样 32 和 16 点,得到 X 32 ( k )和X 16 ( k ) :

X 32 (k ) = X (e jω )

ω=

2π k 32

, k = 0,1, 2,L 31

8

上机实验

X 16 (k ) = X (e jω )

ω=

2π k 16

,

k = 0,1, 2,L15

再分别对 X 32 ( k )和X 16 ( k ) 进行 32 点和 16 点 IFFT,得到 x32 ( n)和x16 ( n) :

x32 (n) = IFFT[ X 32 (k )]32 , n = 0,1, 2,L ,31 x16 (n) = IFFT[ X 16 (k )]16 , n = 0,1, 2,L ,15
分别画出 X (e jω ) 、 X 32 ( k )和X 16 ( k ) 的幅度谱,并绘图显示 x(n)、 x32 ( n)和x16 ( n) 的波形, 进行对比和分析,验证总结频域采样理论。 提示:频域采样用以下方法容易变程序实现。 ① 直接调用 MATLAB 函数 fft 计算 X 32 ( k ) = FFT[ x( n)]32 就得到 X (e jω ) 在 [0,2π ] 的 32 点频率域采样 ② 抽取 X 32 ( k ) 的偶数点即可得到 X (e jω ) 在 [0,2π ] 的 16 点频率域采样 X 16 ( k ) ,即

X 16 (k ) = X 32 (2k ) , k = 0,1, 2,L ,15 。
3 ○ 当然也可以按照频域采样理论, 先将信号 x(n)以 16 为周期进行周期延拓, 取其主值

区(16 点) ,再对其进行 16 点 DFT(FFT),得到的就是 X (e jω ) 在 [0,2π ] 的 16 点频率域采样

X 16 (k ) 。
4.思考题: 如果序列 x(n)的长度为 M,希望得到其频谱 X (e jω ) 在 [0,2π ] 上的 N 点等间隔采样, 当 N<M 时, 如何用一次最少点数的 DFT 得到该频谱采样? 5. 实验报告及要求 a) 运行程序打印要求显示的图形, 。 b) 分析比较实验结果,简述由实验得到的主要结论 c) 简要回答思考题 d) 附上程序清单和有关曲线。

10.2.2 实验程序清单 实验程序清单
1 时域采样理论的验证程序清单 % 时域采样理论验证程序 exp2a.m Tp=64/1000; %观察时间 Tp=64 微秒 %产生 M 长采样序列 x(n) % Fs=1000;T=1/Fs; Fs=1000;T=1/Fs; M=Tp*Fs;n=0:M-1;

9

上机实验

A=444.128;alph=pi*50*2^0.5;omega=pi*50*2^0.5; xnt=A*exp(-alph*n*T).*sin(omega*n*T); Xk=T*fft(xnt,M); %M 点 FFT[xnt)] yn='xa(nT)';subplot(3,2,1); tstem(xnt,yn); %调用自编绘图函数 tstem 绘制序列图 box on;title('(a) Fs=1000Hz'); k=0:M-1;fk=k/Tp; subplot(3,2,2);plot(fk,abs(Xk));title('(a) T*FT[xa(nT)],Fs=1000Hz'); xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))]) %================================================= % Fs=300Hz 和 Fs=200Hz 的程序与上面 Fs=1000Hz 完全相同。 2 频域采样理论的验证程序清单 %频域采样理论验证程序 exp2b.m M=27;N=32;n=0:M; %产生 M 长三角波序列 x(n) xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=[xa,xb]; Xk=fft(xn,1024); %1024 点 FFT[x(n)], 用于近似序列 x(n)的 TF X32k=fft(xn,32) ;%32 点 FFT[x(n)] x32n=ifft(X32k); %32 点 IFFT[X32(k)]得到 x32(n) X16k=X32k(1:2:N); %隔点抽取 X32k 得到 X16(K) x16n=ifft(X16k,N/2); %16 点 IFFT[X16(k)]得到 x16(n) subplot(3,2,2);stem(n,xn,'.');box on title('(b) 三角波序列 x(n)');xlabel('n');ylabel('x(n)');axis([0,32,0,20]) k=0:1023;wk=2*k/1024; % subplot(3,2,1);plot(wk,abs(Xk));title('(a)FT[x(n)]'); xlabel('\omega/\pi');ylabel('|X(e^j^\omega)|');axis([0,1,0,200]) k=0:N/2-1; subplot(3,2,3);stem(k,abs(X16k),'.');box on title('(c) 16 点频域采样');xlabel('k');ylabel('|X_1_6(k)|');axis([0,8,0,200]) n1=0:N/2-1; subplot(3,2,4);stem(n1,x16n,'.');box on title('(d) 16 点 IDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)');axis([0,32,0,20]) k=0:N-1; subplot(3,2,5);stem(k,abs(X32k),'.');box on title('(e) 32 点频域采样');xlabel('k');ylabel('|X_3_2(k)|');axis([0,16,0,200]) n1=0:N-1; subplot(3,2,6);stem(n1,x32n,'.');box on title('(f) 32 点 IDFT[X_3_2(k)]');xlabel('n');ylabel('x_3_2(n)');axis([0,32,0,20])

10.2.3 实验程序运行结果
时域采样理论的验证程序运行结果 1 时域采样理论的验证程序运行结果 exp2a.m 如图 10.3.2 所示。由图可见,采样序列 的频谱的确是以采样频率为周期对模拟信号频谱的周期延拓。当采样频率为 1000Hz 时频谱 混叠很小;当采样频率为 300Hz 时,在折叠频率 150Hz 附近频谱混叠很严重;当采样频率

10

上机实验

为 200Hz 时,在折叠频率 110Hz 附近频谱混叠更很严重。

图 10.2.2 2 时域采样理论的验证程序 exp2b.m 运行结果如图 10.3.3 所示。

11

上机实验

图 10.3.3 该图验证了频域采样理论和频域采样定理。 对信号 x(n)的频谱函数 X(e


)在[0,2π]上等间

隔采样 N=16 时, N 点 IDFT[ X N ( k ) ]得到的序列正是原序列 x(n)以 16 为周期进行周期延 拓后的主值区序列:

xN (n) = IDFT[ X N (k )]N = [ ∑ x(n + iN )]RN (n)
i =?∞



由于 N<M,所以发生了时域混叠失真,因此。 xN ( n) 与 x(n)不相同,如图图 10.3.3(c)和(d) 所示。当 N=32 时,如图图 10.3.3(c)和(d)所示,由于 N>M,频域采样定理,所以不存在时 域混叠失真,因此。 xN ( n) 与 x(n)相同。

10.2.4

简答思考题

先对原序列 x(n)以 N 为周期进行周期延拓后取主值区序列,

xN (n) = [ ∑ x(n + iN )]RN (n)
i =?∞



再计算 N 点 DFT 则得到 N 点频域采样:

X N (k ) = DFT[ xN (n)]N =X (e jω )

ω=

2π k N

,

k = 0,1, 2,L , N ? 1

10.3

实验三: 实验三:用 FFT 对信号作频谱分析

10.3.1 实验指导 实验指导
1.实验目的 学习用 FFT 对连续信号和时域离散信号进行谱分析的方法, 了解可能出现的分析 误差及其原因,以便正确应用 FFT。 2. 实验原理 用 FFT 对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行谱分析的信 号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率 D 和分析误差。 频谱分辨率直接和 FFT 的变换区间 N 有关, 因为 FFT 能够实现的频率分辨率是 2π / N , 因 此要求 2π / N ≤ D 。可以根据此式选择 FFT 的变换区间 N。误差主要来自于用 FFT 作频谱 分析时,得到的是离散谱,而信号(周期信号除外)是连续谱,只有当 N 较大时离散谱的 包络才能逼近于连续谱,因此 N 要适当选择大一些。 周期信号的频谱是离散谱, 只有用整数倍周期的长度作 FFT, 得到的离散谱才能代表周 期信号的频谱。如果不知道信号周期,可以尽量选择信号的观察时间长一些。 对模拟信号进行谱分析时,首先要按照采样定理将其变成时域离散信号。如果是模拟 周期信号,也应该选取整数倍周期的长度,经过采样后形成周期序列,按照周期序列的谱分 析进行。 3.实验步骤及内容

12

上机实验

(1)对以下序列进行谱分析。

x1 (n) = R4 (n) ?n + 1, ? x 2 (n) = ?8 ? n, ?0 , ? ? 4 ? n, ? x3 (n) = ?n ? 3, ?0, ? 0≤n≤3 4≤n≤7

其它n
0≤n≤3 4≤n≤7

其它n

选择 FFT 的变换区间 N 为 8 和 16 两种情况进行频谱分析。分别打印其幅频特性曲线。 并进行对比、分析和讨论。 (2)对以下周期序列进行谱分析。

x4 (n) = cos

π
4

n

x5 (n) = cos(π n / 4) + cos(π n / 8)
选择 FFT 的变换区间 N 为 8 和 16 两种情况分别对以上序列进行频谱分析。分别打印 其幅频特性曲线。并进行对比、分析和讨论。 (3)对模拟周期信号进行谱分析

x6 (t ) = cos8π t + cos16π t + cos 20π t
选择 采样频率 Fs = 64 Hz ,变换区间 N=16,32,64 三种情况进行谱分析。分别打印其幅频 特性,并进行分析和讨论。 4.思考题 (1)对于周期序列,如果周期不知道,如何用 FFT 进行谱分析? (2)如何选择 FFT 的变换区间?(包括非周期信号和周期信号) (3)当 N=8 时, x 2 ( n) 和 x3 ( n) 的幅频特性会相同吗?为什么?N=16 呢? 5.实验报告要求 (1)完成各个实验任务和要求。附上程序清单和有关曲线。 (2)简要回答思考题。

10.3.2 实验程序清单
%第 10 章实验 3 程序 exp3.m % 用 FFT 对信号作频谱分析 clear all;close all %实验内容(1)=================================================== x1n=[ones(1,4)]; M=8;xa=1:(M/2); x3n=[xb,xa]; %产生序列向量 x1(n)=R4(n) xb=(M/2):-1:1; x2n=[xa,xb]; %产生长度为 8 的三角波序列 x2(n)

13

上机实验 %计算 x1n 的 8 点 DFT %计算 x1n 的 16 点 DFT %计算 x1n 的 8 点 DFT %计算 x1n 的 16 点 DFT %计算 x1n 的 8 点 DFT %计算 x1n 的 16 点 DFT

X1k8=fft(x1n,8); X1k16=fft(x1n,16); X2k8=fft(x2n,8); X2k16=fft(x2n,16); X3k8=fft(x3n,8); X3k16=fft(x3n,16); %以下绘制幅频特性曲线

subplot(2,2,1);mstem(X1k8); %绘制 8 点 DFT 的幅频特性图 title('(1a) 8 点 DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X1k8))]) subplot(2,2,3);mstem(X1k16); %绘制 16 点 DFT 的幅频特性图 title('(1b)16 点 DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X1k16))]) figure(2) subplot(2,2,1);mstem(X2k8); %绘制 8 点 DFT 的幅频特性图 title('(2a) 8 点 DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X2k8))]) subplot(2,2,2);mstem(X2k16); %绘制 16 点 DFT 的幅频特性图 title('(2b)16 点 DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X2k16))]) subplot(2,2,3);mstem(X3k8); %绘制 8 点 DFT 的幅频特性图 title('(3a) 8 点 DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X3k8))]) subplot(2,2,4);mstem(X3k16); %绘制 16 点 DFT 的幅频特性图 title('(3b)16 点 DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X3k16))]) %实验内容(2) 周期序列谱分析================================== N=8;n=0:N-1; x4n=cos(pi*n/4); x5n=cos(pi*n/4)+cos(pi*n/8); X4k8=fft(x4n); X5k8=fft(x5n); N=16;n=0:N-1; x4n=cos(pi*n/4); x5n=cos(pi*n/4)+cos(pi*n/8); X4k16=fft(x4n); X5k16=fft(x5n); figure(3) subplot(2,2,1);mstem(X4k8); %绘制 8 点 DFT 的幅频特性图 title('(4a) 8 点 DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X4k8))]) subplot(2,2,3);mstem(X4k16); %绘制 16 点 DFT 的幅频特性图 title('(4b)16 点 DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度'); %计算 x4n 的 16 点 DFT %计算 x5n 的 16 点 DFT %计算 x4n 的 8 点 DFT %计算 x5n 的 8 点 DFT %FFT 的变换区间 N=16 %FFT 的变换区间 N=8

14

上机实验

axis([0,2,0,1.2*max(abs(X4k16))]) subplot(2,2,2);mstem(X5k8); %绘制 8 点 DFT 的幅频特性图 title('(5a) 8 点 DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X5k8))]) subplot(2,2,4);mstem(X5k16); %绘制 16 点 DFT 的幅频特性图 title('(5b)16 点 DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X5k16))]) %实验内容(3) 模拟周期信号谱分析=============================== figure(4) Fs=64;T=1/Fs; N=16;n=0:N-1; %FFT 的变换区间 N=16 %对 x6(t)16 点采样 %计算 x6nT 的 16 点 DFT %将零频率移到频谱中心 %产生 16 点 DFT 对应的采样点频率(以零频率为中心)

x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); X6k16=fft(x6nT); X6k16=fftshift(X6k16); Tp=N*T;F=1/Tp;

%频率分辨率 F

k=-N/2:N/2-1;fk=k*F;

subplot(3,1,1);stem(fk,abs(X6k16),'.');box on %绘制 8 点 DFT 的幅频特性图 title('(6a) 16 点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度'); axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16))]) N=32;n=0:N-1; %FFT 的变换区间 N=16 %对 x6(t)32 点采样 %计算 x6nT 的 32 点 DFT %将零频率移到频谱中心 %产生 16 点 DFT 对应的采样点频率(以零频率为中心)

x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); X6k32=fft(x6nT); X6k32=fftshift(X6k32); Tp=N*T;F=1/Tp;

%频率分辨率 F

k=-N/2:N/2-1;fk=k*F;

subplot(3,1,2);stem(fk,abs(X6k32),'.');box on %绘制 8 点 DFT 的幅频特性图 title('(6b) 32 点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度'); axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32))]) N=64;n=0:N-1; %FFT 的变换区间 N=16 %对 x6(t)64 点采样 %计算 x6nT 的 64 点 DFT %将零频率移到频谱中心 %产生 16 点 DFT 对应的采样点频率(以零频率为中心)

x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); X6k64=fft(x6nT); X6k64=fftshift(X6k64); Tp=N*T;F=1/Tp;

%频率分辨率 F

k=-N/2:N/2-1;fk=k*F;

subplot(3,1,3);stem(fk,abs(X6k64),'.'); box on%绘制 8 点 DFT 的幅频特性图 title('(6a) 64 点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度'); axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64))])

10.3.3 实验程序运行结果
实验 3 程序 exp3.m 运行结果如图 10.3.1 所示。

15

上机实验

16

上机实验

图 10.3.1 程序运行结果分析讨论: 程序运行结果分析讨论: 分析讨论 请读者注意,用 DFT(或 FFT)分析频谱,绘制频谱图时,最好将 X(k)的自变量 k 换 算成对应的频率,作为横坐标便于观察频谱。

ωk =

2π k , k = 0,1, 2,L, N ? 1 N

为了便于读取频率值,最好关于π归一化,即以 1、实验内容(1) 图(1a)和(1b)说明 x1 ( n) ) ) 函数的 8 点和 16 点采样; 因为 x3 ( n)

ω / π 作为横坐标。

= R4 (n) 的 8 点 DFT 和 16 点 DFT 分别是 x1 (n) 的频谱

= x2 ((n + 3))8 R8 (n) , 所以,x3 ( n) 与 x2 ( n) 的 8 点 DFT 的模相等, x (n) 与 x2 (n) 不满足循环移位关系,所

如图(2a)和(3a) 但是,当 N=16 时, 3 。但是 ) ) 但是, 。

以图(2b)和(3b)的模不同。 ) )的模不同。 2、实验内容(2) ,对周期序列谱分析

x4 (n) = cos

π
4

n 的周期为 8,所以 N=8 和 N=16 均是其周期的整数倍,得到正确

的单一频率正弦波的频谱,仅在 0.25π处有 1 根单一谱线。如图(4b)和(4b)所示。

x5 (n) = cos(π n / 4) + cos(π n / 8) 的周期为 16,所以 N=8 不是其周期的整
数倍,得到的频谱不正确,如图(5a)所示。N=16 是其一个周期,得到正确的频谱, 仅在 0.25π和 0.125π处有 2 根单一谱线, 如图(5b)所示。

17

上机实验

3、实验内容(3) ,对模拟周期信号谱分析

x6 (t ) = cos8π t + cos16π t + cos 20π t x6 (t ) 有 3 个频率成分, f1 = 4 Hz , f 2 = 8 Hz , f 3 = 10 Hz 。所以 x6 (t ) 的周
期为 0.5s。 采样频率 Fs

= 64 Hz = 16 f1 = 8 f 2 = 6.4 f 3 。变换区间 N=16 时,观察

时间 Tp=16T=0.25s,不是 x6 (t ) 的整数倍周期,所以所得频谱不正确,如图(6a)所示。 变换区间 N=32,64 时,观察时间 Tp=0.5s,1s,是 x6 (t ) 的整数周期,所以所得频谱正确, 如图(6b)和(6c)所示。图中 3 根谱线正好位于 4 Hz ,8 Hz ,10 Hz 处。变换区间 N=64 时频谱幅度是变换区间 N=32 时 2 倍,这种结果正好验证了用 DFT 对中期序列谱分析的理 论。 注意: 注意: (1)用 DFT(或 FFT)对模拟信号分析频谱时,最好将 X(k)的自变量 k 换算成对应的 ) 模拟频率 fk,作为横坐标绘图,便于观察频谱。这样,不管变换区间 N 取信号周期的几倍, 画出的频谱图中有效离散谐波谱线所在的频率值不变,如图(6b)和(6c)所示。

fk =

Fs 1 1 k= k = k , k = 0,1, 2,L , N ? 1 N NT Tp

(2)本程序直接画出采样序列 N 点 DFT 的模值,实际上分析频谱时最好画出归一化幅 度谱,这样就避免了幅度值随变换区间 N 变化的缺点。本实验程序这样绘图只要是为了验 证了用 DFT 对中期序列谱分析的理论。

10.3.4 简答思考题
思考题(1)和(2)的答案请读者在教材 3.?节找,思考题(3)的答案在程序运行结 果分析讨论已经详细回答。

10.4 实验四 IIR 数字滤波器设计及软件实现
10.4.1 实验指导
1.实验目的 (1)熟悉用双线性变换法设计 IIR 数字滤波器的原理与方法; (2)学会调用 MATLAB 信号处理工具箱中滤波器设计函数(或滤波器设计分析工具 fdatool)设计各种 IIR 数字滤波器,学会根据滤波需求确定滤波器指标参数。 (3)掌握 IIR 数字滤波器的 MATLAB 实现方法。 (3)通过观察滤波器输入输出信号的时域波形及其频谱,建立数字滤波的概念。 2.实验原理 设计 IIR 数字滤波器一般采用间接法(脉冲响应不变法和双线性变换法) ,应用最广泛 的是双线性变换法。 基本设计过程是: ①先将给定的数字滤波器的指标转换成过渡模拟滤波 器的指标; ②设计过渡模拟滤波器;③将过渡模拟滤波器系统函数转换成数字滤波器的系 统函数。 MATLAB 信号处理工具箱中的各种 IIR 数字滤波器设计函数都是采用双线性变换法。

18

上机实验

第六章介绍的滤波器设计函数 butter、cheby1 、cheby2 和 ellip 可以分别被调用来直接 设计巴特沃斯、切比雪夫 1、切比雪夫 2 和椭圆模拟和数字滤波器。本实验要求读者调用如 上函数直接设计 IIR 数字滤波器。 本实验的数字滤波器的 MATLAB 实现是指调用 MATLAB 信号处理工具箱函数 filter 对给 定的输入信号 x(n)进行滤波,得到滤波后的输出信号 y(n) 。 3. 实验内容及步骤 (1)调用信号产生函数 mstg 产生由三路抑制载波调幅信号相加构成的复合信号 st, 该函数还会自动绘图显示 st 的时域波形和幅频特性曲线,如图 10.4.1 所示。由图可见,三 路信号时域混叠无法在时域分离。但频域是分离的,所以可以通过滤波的方法在频域分离, 这就是本实验的目的。

图 10.4.1 三路调幅信号 st 的时域波形和幅频特性曲线 (2)要求将 st 中三路调幅信号分离,通过观察 st 的幅频特性曲线,分别确定可以分 离 st 中三路抑制载波单频调幅信号的三个滤波器(低通滤波器、带通滤波器、高通滤波器) 的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为 0.1dB,阻带最小衰减为 60dB。 提示:抑制载波单频调幅信号的数学表示式为

1 s (t ) = cos(2π f 0t ) cos(2π f c t ) = [cos(2π ( f c ? f 0 )t ) + cos(2π ( f c + f 0 )t )] 2
其中, cos(2π f c t ) 称为载波,fc 为载波频率, cos(2π f 0t ) 称为单频调制信号,f0 为调制正 弦波信号频率,且满足 f c > f 0 。由上式可见,所谓抑制载波单频调幅信号,就是 2 个正弦 信号相乘,它有 2 个频率成分:和频 f c + f 0 和差频 f c ? f 0 ,这 2 个频率成分关于载波频率 fc 对称。所以,1 路抑制载波单频调幅信号的频谱图是关于载波频率 fc 对称的 2 根谱线,其 中没有载频成分,故取名为抑制载波单频调幅信号。容易看出,图 10.4.1 中三路调幅信号 的载波频率分别为 250Hz、500Hz、1000Hz。如果调制信号 m(t)具有带限连续频谱,无直流 成分,则 s (t ) = m(t ) cos(2π f c t ) 就是一般的抑制载波调幅信号。其频谱图是关于载波频率

19

上机实验

fc 对称的 2 个边带(上下边带) ,在专业课通信原理中称为双边带抑制载波 (DSB-SC) 调幅 信号,简称双边带 (DSB) 信号。 如果调制信号 m(t)有直流成分, s (t ) = m(t ) cos(2π f c t ) 就 则 是一般的双边带调幅信号。其频谱图是关于载波频率 fc 对称的 2 个边带(上下边带) ,并包 含载频成分。 (3)编程序调用 MATLAB 滤波器设计函数 ellipord 和 ellip 分别设计这三个椭圆滤波 器,并绘图显示其幅频响应特性曲线。 (4)调用滤波器实现函数 filter,用三个滤波器分别对信号产生函数 mstg 产生的信 号 st 进行滤波,分离出 st 中的三路不同载波频率的调幅信号 y1(n)、y2(n)和 y3(n), 并绘 图显示 y1(n)、y2(n)和 y3(n)的时域波形,观察分离效果。 4.信号产生函数 mstg 清单 function st=mstg %产生信号序列向量 st,并显示 st 的时域波形和频谱 %st=mstg 返回三路调幅信号相加形成的混合信号,长度 N=1600 N=1600 %N 为信号 st 的长度。 Fs=10000;T=1/Fs;Tp=N*T; %采样频率 Fs=10kHz,Tp 为采样时间 t=0:T:(N-1)*T;k=0:N-1;f=k/Tp; fc1=Fs/10; %第 1 路调幅信号的载波频率 fc1=1000Hz, fm1=fc1/10; %第 1 路调幅信号的调制信号频率 fm1=100Hz fc2=Fs/20; %第 2 路调幅信号的载波频率 fc2=500Hz fm2=fc2/10; %第 2 路调幅信号的调制信号频率 fm2=50Hz fc3=Fs/40; %第 3 路调幅信号的载波频率 fc3=250Hz, fm3=fc3/10; %第 3 路调幅信号的调制信号频率 fm3=25Hz xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第 1 路调幅信号 xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第 2 路调幅信号 xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第 3 路调幅信号 st=xt1+xt2+xt3; %三路调幅信号相加 fxt=fft(st,N); %计算信号 st 的频谱 %====以下为绘图部分,绘制 st 的时域波形和幅频特性曲线==================== subplot(3,1,1) plot(t,st);grid;xlabel('t/s');ylabel('s(t)'); axis([0,Tp/8,min(st),max(st)]);title('(a) s(t)的波形') subplot(3,1,2) stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b) s(t)的频谱') axis([0,Fs/5,0,1.2]); xlabel('f/Hz');ylabel('幅度') 5.实验程序框图如图 10.4.2 所示,供读者参考。

20

上机实验

调用函数 mstg 产生 st,自动绘图 显示 st 的时域波形和幅频特性曲线

调用 ellipord 和 ellip 分别设计三个椭圆滤 波器,并绘图显示其幅频响应特性曲线。

调用 filter,用三个滤波器分别对信号 st 进行滤波,分离 出三路不同载波频率的调幅信号 y1(n)、y2(n)和 y3(n)

绘图显示 y1(n)、y2(n)和 y3(n)的时域波形和幅频特性曲线

End 图 10.4.2 实验 4 程序框图 6.思考题 (1)请阅读信号产生函数 mstg,确定三路调幅信号的载波频率和调制信号频率。 (2)信号产生函数 mstg 中采样点数 N=800,对 st 进行 N 点 FFT 可以得到 6 根理想谱 线。 如果取 N=1000, 可否得到 6 根理想谱线?为什么?N=2000 呢?请改变函数 mstg 中采样 点数 N 的值,观察频谱图验证您的判断是否正确。 (3)修改信号产生函数 mstg,给每路调幅信号加入载波成分,产生调幅(AM)信号, 重复本实验,观察 AM 信号与抑制载波调幅信号的时域波形及其频谱的差别。 提示:AM 信号表示式: s (t ) = [1 + cos(2π f 0t )]cos(2π f ct ) 。 7.实验报告要求 (1)简述实验目的及原理。 (2)画出实验主程序框图,打印程序清单。 (3)绘制三个分离滤波器的损耗函数曲线。 (4)绘制经过滤波分理出的三路调幅信号的时域波形。 (5)简要回答思考题。

10.4.2 滤波器参数及实验程序清单 滤波器参数及实验程序清单 实验
1、滤波器参数选取 、 观察图 10.4.1 可知, 三路调幅信号的载波频率分别为 250Hz、 500Hz、 1000Hz。 带宽 (也 可以由信号产生函数 mstg 清单看出)分别为 50Hz、100Hz、200Hz。所以,分离混合信号 st 中三路抑制载波单频调幅信号的三个滤波器(低通滤波器、带通滤波器、高通滤波器)的指 标参数选取如下: 对载波频率为 250Hz 的条幅信号,可以用低通滤波器分离,其指标为

21

上机实验

带截止频率

f p = 280 Hz,通带最大衰减α p = 0.1dB dB;
= 450 Hz,阻带最小衰减α s = 60dB dB,

阻带截止频率 f s

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

f pl = 440 Hz, f pu = 560 Hz,通带最大衰减α p = 0.1dB dB; = 275 Hz, f su = 900 Hz,Hz,阻带最小衰减α s = 60dB dB,

阻带截止频率 f sl

对载波频率为 1000Hz 的条幅信号,可以用高通滤波器分离,其指标为 带截止频率

f p = 890 Hz,通带最大衰减α p = 0.1dB dB; = 550 Hz,阻带最小衰减α s = 60dB dB,

阻带截止频率 f s

说明: (1)为了使滤波器阶数尽可能低,每个滤波器的边界频率选择原则是尽量使滤波 器过渡带宽尽可能宽。 (2)与信号产生函数 mstg 相同,采样频率 Fs=10kHz。 (3)为了滤波器阶数最低,选用椭圆滤波器。 按照图 10.4.2 所示的程序框图编写的实验程序为 exp4.m。 2、实验程序清单
%实验 4 程序 exp4.m % IIR 数字滤波器设计及软件实现 clear all;close all Fs=10000;T=1/Fs; st=mstg; %低通滤波器设计与实现========================================= fp=280;fs=450; wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; [B,A]=ellip(N,rp,rs,wp); y1t=filter(B,A,st); % 低通滤波器设计与实现绘图部分 figure(2);subplot(3,1,1); myplot(B,A); yt='y_1(t)'; subplot(3,1,2);tplot(y1t,T,yt); %调用绘图函数 tplot 绘制滤波器输出波形 %带通滤波器设计与实现==================================================== fpl=440;fpu=560;fsl=275;fsu=900; wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60; [N,wp]=ellipord(wp,ws,rp,rs); y2t=filter(B,A,st); %调用 ellipord 计算椭圆 DF 阶数 N 和通带截止频率 wp [B,A]=ellip(N,rp,rs,wp); %调用 ellip 计算椭圆带通 DF 系统函数系数向量 B 和 A %滤波器软件实现 % 带通滤波器设计与实现绘图部分(省略) %调用绘图函数 myplot 绘制损耗函数曲线 %DF 指标(低通滤波器的通、阻带边界频) [N,wp]=ellipord(wp,ws,rp,rs); %调用 ellipord 计算椭圆 DF 阶数 N 和通带截止频率 wp %调用 ellip 计算椭圆带通 DF 系统函数系数向量 B 和 A %滤波器软件实现 %采样频率 %调用信号产生函数 mstg 产生由三路抑制载波调幅信号相加构成的复合信号 st

22

上机实验 %高通滤波器设计与实现================================================ fp=890;fs=600; wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; [N,wp]=ellipord(wp,ws,rp,rs); y3t=filter(B,A,st); %DF 指标(低通滤波器的通、阻带边界频) %调用 ellipord 计算椭圆 DF 阶数 N 和通带截止频率 wp

[B,A]=ellip(N,rp,rs,wp,'high'); %调用 ellip 计算椭圆带通 DF 系统函数系数向量 B 和 A %滤波器软件实现 % 高低通滤波器设计与实现绘图部分(省略)

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

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

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

23

上机实验

(c)高通滤波器损耗函数及其分离出的调幅信号 y3(t) 图 104. 实验 4 程序 exp4.m 运行结果

10.4.4 简要回答思考题 简要回答思考题
思考题(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 数字滤波器设计与软件实现 实验五: 数字滤波器设计与软件 软件实现
10.5.1 实验指导
1.实验目的 (1)掌握用窗函数法设计 FIR 数字滤波器的原理和方法。 (2)掌握用等波纹最佳逼近法设计 FIR 数字滤波器的原理和方法。 (3)掌握 FIR 滤波器的快速卷积实现原理。 (4)学会调用 MATLAB 函数设计与实现 FIR 滤波器。 2. 实验内容及步骤 (1)认真复习第七章中用窗函数法和等波纹最佳逼近法设计 FIR 数字滤波器的原理; (2)调用信号产生函数 xtg 产生具有加性噪声的信号 xt,并自动显示 xt 及其频谱,如 图 10.5.1 所示;

24

上机实验

图 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 数字滤波器。并比较两种设计方法设计的滤波器阶数。 提示: 提示: 1 MATLAB 函数 fir1 和 fftfilt 的功能及其调用格式请查阅本书第 7 章和第?章; ○ 2 ○采样频率 Fs=1000Hz,采样周期 T=1/Fs; 3 ○根据图 10.6.1(b)和实验要求,可选择滤波器指标参数:通带截止频率 fp=120Hz,阻 带截至频率 fs=150Hz,换算成数字频率,通带截止频率 ωp = 2π f p Τ = 0.24π ,通带最大衰为
0.1dB,阻带截至频率 ωs = 2π f s Τ = 0.3π ,阻带最小衰为 60dB。]

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

25

上机实验

Fs=1000, T=1/Fs

xt=xtg 产生信号 xt, 并显示 xt 及其频谱

用窗函数法或等波纹最佳逼近法 设计 FIR 滤波器 hn

对信号 xt 滤波:yt=fftfilt(hn,xt)

1、计算并绘图显示滤波器损耗函数 2、绘图显示滤波器输出信号 yt

End

图 10.5.2 实验程序框图 4.思考题 (1)如果给定通带截止频率和阻带截止频率以及阻带最小衰减,如何用窗函数法设计线 性相位低通滤波器?请写出设计步骤. (2)如果要求用窗函数法设计带通滤波器,且给定通带上、下截止频率为

ωpl 和 ωpu ,

阻带上、下截止频率为 ωsl 和 ωsu ,试求理想带通滤波器的截止频率 ωcl 和ωcu 。 (3)解释为什么对同样的技术指标,用等波纹最佳逼近法设计的滤波器阶数低? 5.实验报告要求 (1)对两种设计 FIR 滤波器的方法(窗函数法和等波纹最佳逼近法)进行分析比较, 简述其优缺点。 (2)附程序清单、打印实验内容要求绘图显示的曲线图。 (3)分析总结实验结果。 (4)简要回答思考题。 6.信号产生函数 xtg 程序清单 function xt=xtg(N) %实验五信号 x(t)产生,并显示信号的幅频特性曲线 %xt=xtg(N) 产生一个长度为 N,有加性高频噪声的单频调幅信号 xt,采样频率 Fs=1000Hz %载波频率 fc=Fs/10=100Hz,调制正弦波频率 f0=fc/10=10Hz. Fs=1000;T=1/Fs;Tp=N*T; t=0:T:(N-1)*T;

26

上机实验

fc=Fs/10;f0=fc/10; %载波频率 fc=Fs/10,单频调制信号频率为 f0=Fc/10; mt=cos(2*pi*f0*t); %产生单频正弦波调制信号 mt,频率为 f0 ct=cos(2*pi*fc*t); %产生载波正弦波信号 ct,频率为 fc xt=mt.*ct; %相乘产生单频调制信号 xt nt=2*rand(1,N)-1; %产生随机噪声 nt %=======设计高通滤波器 hn,用于滤除噪声 nt 中的低频成分,生成高通噪声======= fp=150; fs=200;Rp=0.1;As=70; % 滤波器指标 fb=[fp,fs];m=[0,1]; % 计算 remezord 函数所需参数 f,m,dev dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)]; [n,fo,mo,W]=remezord(fb,m,dev,Fs); % 确定 remez 函数所需参数 hn=remez(n,fo,mo,W); % 调用 remez 函数进行设计,用于滤除噪声 nt 中的低频成分 yt=filter(hn,1,10*nt); %滤除随机噪声中低频成分,生成高通噪声 yt %================================================================ xt=xt+yt; %噪声加信号 fst=fft(xt,N);k=0:N-1;f=k/Tp; subplot(3,1,1);plot(t,xt);grid;xlabel('t/s');ylabel('x(t)'); axis([0,Tp/5,min(xt),max(xt)]);title('(a) 信号加噪声波形') subplot(3,1,2);plot(f,abs(fst)/max(abs(fst)));grid;title('(b) 信号加噪声的频谱') axis([0,Fs/2,0,1.2]);xlabel('f/Hz');ylabel('幅度')

10.5.2 滤波器参数及实验程序清单 滤波器参数及实验程序清单 实验
1、滤波器参数选取 、 根据 10.5.1 节实验指导的提示③选择滤波器指标参数:通带截止频率 fp=120Hz,阻带 截 至 频 率 fs=150Hz 。 代 入 采 样 频 率 Fs=1000Hz , 换 算 成 数 字 频 率 , 通 带 截 止 频 率

ωp = 2π f p Τ = 0.24π ,通带最大衰为 0.1dB,阻带截至频率 ωs = 2π f s Τ = 0.3π ,阻带最小衰为
60dB。所以选取 blackman 窗函数。与信号产生函数 xtg 相同,采样频率 Fs=1000Hz。 按照图 10.5.2 所示的程序框图编写的实验程序为 exp5.m。 2、实验程序清单
%《数字信号处理(第三版)学习指导》第 10 章实验 5 程序 exp5.m % FIR 数字滤波器设计及软件实现 clear all;close all; %==调用 xtg 产生信号 xt, xt 长度 N=1000,并显示 xt 及其频谱,========= N=1000;xt=xtg(N); fp=120; fs=150;Rp=0.2;As=60;Fs=1000; % (1) 用窗函数法设计滤波器 wc=(fp+fs)/Fs; B=2*pi*(fs-fp)/Fs; Nb=ceil(11*pi/B); %理想低通滤波器截止频率(关于 pi 归一化) %过渡带宽度指标 %blackman 窗的长度 N % 输入给定指标

hn=fir1(Nb-1,wc,blackman(Nb)); Hw=abs(fft(hn,1024)); % 求设计的滤波器频率特性 ywt=fftfilt(hn,xt,N); %调用函数 fftfilt 对 xt 滤波

%以下为用窗函数法设计法的绘图部分(滤波器损耗函数,滤波器输出信号波形) %省略 % (2) 用等波纹最佳逼近法设计滤波器

27

上机实验 % 确定 remezord 函数所需参数 f,m,dev % 确定 remez 函数所需参数

fb=[fp,fs];m=[1,0];

dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; [Ne,fo,mo,W]=remezord(fb,m,dev,Fs); hn=remez(Ne,fo,mo,W); Hw=abs(fft(hn,1024)); yet=fftfilt(hn,xt,N); % 调用 remez 函数进行设计 % 求设计的滤波器频率特性 % 调用函数 fftfilt 对 xt 滤波

%以下为用等波纹设计法的绘图部分(滤波器损耗函数,滤波器输出信号 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)可以直观 地看出时延差别。

图 10.5.3

28

上机实验

10.5.4 简答思考题 简答思考题
(1) 用窗函数法设计线性相位低通滤波器的设计步骤教材中有详细的介绍. (2) 希望逼近的理想带通滤波器的截止频率 ωcl 和ωcu 分别为:

ωcl = (ωsl + ωpl ) / 2, ωcu = (ωsu + ωpu ) / 2
(3)解释为什么对同样的技术指标,用等波纹最佳逼近法设计的滤波器阶数低? ①用窗函数法设计的滤波器,如果在阻带截止频率附近刚好满足,则离开阻带截止频 率越远,阻带衰减富裕量越大,即存在资源浪费; ② 几种常用的典型窗函数的通带最大衰减和阻带最小衰减固定,且差别较大,又不能 分别控制。 所以设计的滤波器的通带最大衰减和阻带最小衰减通常都存在较大富裕。 如本实 验所选的 blackman 窗函数,其阻带最小衰减为 74dB,而指标仅为 60dB。 ③ 用等波纹最佳逼近法设计的滤波器,其通带和阻带均为等波纹特性,且通带最大衰 减和阻带最小衰减可以分别控制, 所以其指标均匀分布, 没有资源浪费, 所以期阶数低得多。

10.6
10.6.1

实验六 数字信号处理在双音多频拨号系统中的应用
实验指导

1、引言 、 双音多频(Dual Tone Multi Frequency, 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 两个频率, 信号用 sin( 2πf1t ) + sin( 2πf 2 t ) 表示, 其中 f 1 = 679 Hz ,

f 2 = 1209 Hz 。这样 8 个频率形成 16 种不同的双频信号。具体号码以及符号对应的频率如
表 10.6.1 所示。表中最后一列在电话中暂时未用。

表 10.6.1 列 行

双频拨号的频率分配 1336Hz 1477Hz 633Hz

1209Hz

29

上机实验

697Hz 770Hz 852Hz 942Hz

1 4 7 *

2 5 8 0

3 6 9 #

A B C D

DTMF 信号在电话中有两种作用,一个是用拨号信号去控制交换机接通被叫的用户电话 机,另一个作用是控制电话机的各种动作,如播放留言、语音信箱等。 2 电话中的双音多频(DTMF)信号的产生与检测 电话中的双音多频( ) (1)双音多频信号的产生 假设时间连续的 DTMF 信号用 x(t ) = sin( 2πf 1t ) + sin( 2πf 2 t ) 表示, 式中 f 1和f 2 是按 照表 10.10.1 选择的两个频率, f1 代表低频带中的一个频率, f 2 代表高频带中的一个频率。 显然采用数字方法产生 DTMF 信号,方便而且体积小。下面介绍采用数字方法产生 DTMF 信号。规定用 8KHz 对 DTMF 信号进行采样,采样后得到时域离散信号为

x(n) = sin(2πf 1 n / 8000) + sin(2πf 2 n / 8000)
形成上面序列的方法有两种,即计算法和查表法。用计算法求正弦波的序列值容易,但实 际中要占用一些计算时间,影响运行速度。查表法是预先将正弦波的各序列值计算出来,寄 存在存储器中, 运行时只要按顺序和一定的速度取出便可。 这种方法要占用一定的存储空间, 但是速度快。 因为采样频率是 8000Hz,因此要求每 125ms 输出一个样本,得到的序列再送到 D/A 变 换器和平滑滤波器,输出便是连续时间的 DTMF 信号。DTMF 信号通过电话线路送到交换 机。 (2)双音多频信号的检测 在接收端,要对收到的双音多频信号进行检测,检测两个正弦波的频率是多少,以判 断所对应的十进制数字或者符号。 显然这里仍然要用数字方法进行检测, 因此要将收到的时 间连续 DTMF 信号经过 A/D 变换,变成数字信号进行检测。检测的方法有两种,一种是用 一组滤波器提取所关心的频率, 根据有输出信号的 2 个滤波器判断相应的数字或符号。 另一 种是用 DFT (FFT) 对双音多频信号进行频谱分析, 由信号的幅度谱, 判断信号的两个频率, 最后确定相应的数字或符号。当检测的音频数目较少时,用滤波器组实现更合适。FFT 是 DFT 的快速算法,但当 DFT 的变换区间较小时,FFT 快速算法的效果并不明显,而且还要 占用很多内存,因此不如直接用 DFT 合适。下面介绍 Goertzel 算法,这种算法的实质是直 接计算 DFT 的一种线性滤波方法。这里略去 Goertzel 算法的介绍(请参考文献[19]) ,可以 直接调用 MATLAB 信号处理工具箱中戈泽尔算法的函数 Goertzel, 计算 N 点 DFT 的几个感 兴趣的频点的值。 3 检测 DTMF 信号的 DFT 参数选择 用 DFT 检测模拟 DTMF 信号所含有的两个音频频率,是一个用 DFT 对模拟信号进行频 谱分析的问题。根据第三章用 DFT 对模拟信号进行谱分析的理论,确定三个参数: (1)采 样频率 Fs , (2)DFT 的变换点数 N, (3)需要对信号的观察时间的长度 T p 。这三个参数不 能随意选取,要根据对信号频谱分析的要求进行确定。这里对信号频谱分析也有三个要求: (1)频率分辨率, (2)谱分析的频谱范围, (3)检测频率的准确性。

30

上机实验

1. 频谱分析的分辨率。 观察要检测的 8 个频率,相邻间隔最小的是第一和第二个频率,间隔是 73Hz,要求 DFT 最少能够分辨相隔 73Hz 的两个频率,即要求 Fmin = 73Hz 。DFT 的分辨率和对信号的 观察时间 T p 有关, T p min = 1 / F = 1 / 73 = 13.7 ms 。考虑到可靠性,留有富裕量,要求按 键的时间大于 40ms。 2 频谱分析的频率范围 要检测的信号频率范围是 697~1633Hz,但考虑到存在语音干扰,除了检测这 8 个频率 外, 还要检测它们的二次倍频的幅度大小, 波形正常且干扰小的正弦波的二次倍频是很小的, 如果发现二次谐波很大,则不能确定这是 DTMF 信号。这样频谱分析的频率范围为 697~ 3266Hz。按照采样定理,最高频率不能超过折叠频率,即 0.5 Fs ≥ 3622 Hz ,由此要求最小 的采样频率应为 7.24KHz。因为数字电话总系统已经规定 Fs =8KHz,因此对频谱分析范围 的要求是一定满足的。按照 T p min = 13.7 ms , Fs =8KHz,算出对信号最少的采样点数为

N min = T p min ? Fs ≈ 110 。
检测频率的准确性 这是一个用 DFT 检测正弦波频率是否准确的问题。序列的 N 点 DFT 是对序列频谱函 数在 0~ 2π 区间的 N 点等间隔采样,如果是一个周期序列,截取周期序列的整数倍周期, 进行 DFT,其采样点刚好在周期信号的频率上,DFT 的幅度最大处就是信号的准确频率。 分析这些 DTMF 信号,不可能经过采样得到周期序列,因此存在检测频率的准确性问题。 3 DFT 的频率采样点频率为 ω k = 2πk / N (k=0,1,2,---,N-1) ,相应的模拟域采样点频率为

f k = Fs k / N (k=0,1,2,---,N-1) ,希望选择一个合适的 N,使用该公式算出的 f k 能接近要检
测的频率,或者用 8 个频率中的任一个频率 f k 代入公式 f k
' '

= Fs k / N 中时,得到的

k

值最接近整数值,这样虽然用幅度最大点检测的频率有误差,但可以准确判断所对应的 DTMF 频率,即可以准确判断所对应的数字或符号。经过分析研究认为 N=205 是最好的。 按照 Fs =8KHz,N=205,算出 8 个频率及其二次谐波对应 k 值,和 k 取整数时的频率误差 见表 10.6.2。 表 10.6.2 8 个基频 Hz 697 770 852 941 1209 1336 最 近 的 整数 k 值 17.861 19.531 21.833 24.113 30.981 34.235 DFT 的 k值 18 20 22 24 31 34 绝对误差 0.139 0.269 0.167 0.113 0.019 0.235 二次谐波 Hz 1394 1540 1704 1882 2418 2672 对应的 k值 35.024 38.692 42.813 47.285 60.752 67.134 最近的 整数 k 值 35 39 43 47 61 67 绝对误差 0.024 0.308 0.187 0.285 0.248 0.134
31

上机实验

1477 1633

37.848 41.846

38 42

0.152 0.154

2954 3266

74.219 82.058

74 82

0.219 0.058

通过以上分析,确定 Fs =8KHz,N=205, T p ≥ 40ms 。 4 DTMF 信号的产生与识别仿真实验 信号的产生与识别仿真 产生与识别仿真实验 下面先介绍 MATLAB 工具箱函数 goertzel,然后介绍 DTMF 信号的产生与识别仿真实 验程序。Goerztel 函数的调用格式额为 Xgk=goertzel(xn,K) xn 是被变换的时域序列,用于 DTMF 信号检测时,xn 就是 DTMF 信号的 205 个采样值。 K 是要求计算的 DFT[xn]的频点序号向量, N 表示 xn 的长度, 用 则要求 1≤K≤N。 由表 10.2.2 可知,如果只计算 DTMF 信号 8 个基频时, K=[18,20,22,24,31,34,38,42], 如果同时计算 8 个基频及其二次谐波时, K=[18,20,22,24,31,34,35,38,39,42,43,47,61,67,74,82]。 Xgk 是变换结果向量,其中存放的是由 K 指定的频率点的 DFT[x(n)]的值。设 X(k)= DFT[x(n)],则 Xgk (i )

= X ( K (i)), i = 1,2,L,length( K ) 。

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 双频拨号信号的生成和检测程序 %clear all;clc; tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68]; N=205;K=[18,20,22,24,31,34,38,42]; f1=[697,770,852,941]; f2=[1209,1336,1477,1633]; TN=input('键入 6 位电话号码= '); TNr=0; for l=1:6; d=fix(TN/10^(6-l)); TN=TN-d*10^(6-l); for p=1:4; for q=1:4; if tm(p,q)==abs(d); break,end end
32

% DTMF 信号代表的 16 个数 % 行频率向量 % 列频率向量 % 输入 6 位数字 %接收端电话号码初值为零

% 检测码相符的列号 q

上机实验

if tm(p,q)==abs(d); break,end end n=0:1023; sound(x,8000); pause(0.1) % 接收检测端的程序 X=goertzel(x(1:205),K+1); val = abs(X); subplot(3,2,l);

% 检测码相符的行号 p % 为了发声,加长序列 % 发出声音

x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000);% 构成双频信号

% 用 Goertzel 算法计算八点 DFT 样本 % 列出八点 DFT 向量

stem(K,val,'.');grid;xlabel('k');ylabel('|X(k)|') % 画出 DFT(k)幅度 axis([10 50 0 120]) limit = 80; for s=5:8; if val(s) > limit, break, end end for r=1:4; if val(r) > limit, break, end end TNr=TNr+tm(r,s-4)*10^(6-l); end disp('接收端检测到的号码为:') disp(TNr) 运行程序,根据提示键入 6 位电话号码 123456,回车后可以听见 6 位电话号码对应的 DTMF 信号的声音,并输出相应的 6 幅频谱图如图 10.10.1 所示,左上角的第一个图在 k=18 和 k=31 两点出现峰值,所以对应第一位号码数字 1。最后显示检测到的电话号码 123456。 % 显示接收到的字符 % 查找行号 % 查找列号 %

33

上机实验

图 10.6.1 6 位电话号码 123456 的 DTMF 信号在 8 个近似基频点的 DFT 幅度 (1) 实验内容 ① 运行仿真程序 exp6.m,任意送入 6 位电话号码,打印出相应的幅度谱。观察程 序运行结果,对照表 10.10.1 和表 10.10.2,判断程序谱分析的正确性。 ② 分析该仿真程序,将产生、检测和识别 6 位电话号码的程序改为能产生、检测和 识别 8 位电话号码的程序,并运行一次,打印出相应的幅度谱和 8 位电话号码。 5. 实验报告 (1) 分析程序 exp8.m,画出仿真程序流程图。 (2) 打印 6 位和 8 位电话号码 DTMF 信号的幅度谱。 (3) 简述 DTMF 信号的参数:采样频率、DFT 的变换点数以及观测时间的确定原则。

10.6.2 实验程序清单及运行结果 实验程序清单及运行结果 程序清单
1、实验内容① 6 位电话号码的 DTMF 双频拨号信号的生成和检测程序清单 exp6.m 已 经在实验指导中给出。运行程序,并输入 6 位电话号码 123456,则输出相应的 6 幅频谱图 如图 10.10.1 所示,左上角的第一个图在 k=18 和 k=31 两点出现峰值,所以对应第一位号码 数字 1。其他 5 个图请读者对照表 10.10.1 和表 10.10.2,确定确定其对应的数字,验证程序 输出的电话号码“123456”是正确的。 2、实验内容② 只要对 6 位电话号码检测程序 exp6.m 作如下修改,即可产生、检测和 识别 8 位电话号码。 (1)将第 8 行改为 TN=input('键入 8 位电话号码= '); (2)将第 10~12 行改为 for l=1:8; d=fix(TN/10^(8-l)); TN=TN-d*10^(8-l); (3)将第 26 行改为 subplot(4,2,l);

34

上机实验

(4)将第 36 行改为 TNr=TNr+tm(r,s-4)*10^(8-l); 修改后的程序为 exp6_8.m,程序清单见程序集。运行程序 exp6_8.m,输入输入 8 位 电话号码 87654321,则输出相应的 8 幅频谱图如图 10.10.2 所示。最后显示检测到的 电话号码 87654321。

100 50 0 10 100 50 0 10 100 50 0 10

20

30 k

40

50

100 50 0 10 100 50 0 10 100 50 0 10

|X(k)|

|X(k)|

20

30 k

40

50

|X(k)|

20

30 k

40

50

|X(k)|

20

30 k

40

50

|X(k)|

20

30 k

40

50

|X(k)|

20

30 k

40

50

100 |X(k)| 50 0 10 20 30 k 40 50 |X(k)|

100 50 0 10 20 30 k 40 50

图 10.6.1 8 位电话号码 87654321 的 DTMF 信号在 8 个近似基频点的 DFT 幅度

35


赞助商链接
相关文章:
数字信号处理第三版上机实验答案 (1)
数字信号处理第三版上机实验答案 (1)_理学_高等教育_教育专区。实验指导 及分析 第十章上机实验 10.1 实验一: 系统响应及系统稳定性 1.实验目的 (1)掌握 求...
数字信号处理上机实验答案(第三版,第十章)
数字信号处理上机实验答案(第三版,第十章) - 10.1.2 实验程序清单 %实验 1:系统响应及系统稳定性 实验 : close all;clear all %===内容 1:调用 fi...
数字信号处理上机实验答案(第三版,第十章)
第十章上机实验 第十章 上机实验数字信号处理是一门理论和实际密切结合的课程,为深入掌握课程内容,最好在学习理 论的同时,做习题和上机实验上机实验不仅可以帮助...
数字信号处理上机实验答案(全)1
数字信号处理上机实验答案(全)1_数学_自然科学_专业资料。第十章上机实验 第十...2、实验程序清单 %《数字信号处理(第三版)学习指导》第 10 章实验 5 程序 ...
数字信号处理上机答案(含程序及图片)第三版高西全著
数字信号处理上机答案(含程序及图片)第 三版高西全 丁玉美著 数字信号处理实验一内容一 a=0.8;ys=0; A=[1,-0.9];B=[0.05,0.05]; xn=[1,zeros(1,...
数字信号处理第三版高西泉上机实验完整版(含子程序) - ...
数字信号处理第三版高西泉上机实验完整版(含子程序) - 副本_理学_高等教育_教育专区 暂无评价|0人阅读|0次下载|举报文档 数字信号处理第三版高西泉上机实验完整...
数字信号处理第三版高西泉上机实验完整版(含子程序)
数字信号处理第三版高西泉上机实验完整版(含子程序)_理学_高等教育_教育专区。实验一:系统响应及系统稳定性一、 实验原理与方法 1、在时域求系统响应的方法有两种...
数字信号处理第三版高西泉上机实验完整版(含子程序)_-_...
数字信号处理第三版高西泉上机实验完整版(含子程序)_-_副本 - 实验一:系统响应及系统稳定性 一、 实验原理与方法 1、在时域求系统响应的方法有两种,第一种是...
数字信号处理上机实验 作业结果与说明 实验三、四、五
数字信号处理上机实验 作业结果与说明 实验三、四、五_工学_高等教育_教育专区。数字信号处理上机实验_作业结果与说明(实验三、四、五)_赵晓磊_1120130376 ...
数字信号处理第三版用MATLAB上机实验
数字信号处理第三版用MATLAB上机实验 - 实验二:时域采样与频域采样。实验三:用FFT对信号做频谱分析。
更多相关标签: