数字信号处理 实验3 离散系统的频域分析

function［db,mag,pha,grd,w］＝freqz_m(b,a,N)； %离散系统响应求解函数,包括幅度响应、相位响应、群时延; %a,b,N 为输入参数; %b 为系统分子多项式系数向量; %b 为系统分子多项式系数向量; N 为频率离散化的点数; %db,mag,pha,grd,w 为输出参数; %db 为系统[0，π]的度响应,单位是 dB; %ag 也是系统[0，π]的幅度响应,单位是伏; %pha 为系统的[0，π]上的相位响应,单位为 rad; %grd 为[0，π]上的群延时响应; %w 为[0，π]上的频率采样点. [H，w]=freqz(b,a,N,’whole’); H=(H(1:N/2));w=(w(1:N/2)); mag＝abs(H); db＝20*log10((mag＋eps)/max(mag)); pha＝angle(H); grd＝grpdelay(b,a,w);
MATLAB 为求解离散系统的频率响应和连续系统的频率响应，分别提供了 freqz 和 freq（s 求
连续系统的频率响应函数）两个函数，使用方法类似。本实验主要讨论离散系统的频率响应。
例 3-1 已知离散时间系统的系统函数为
H(z)
=
0.1321− 0.3963 z−2 + 0.3963 z−4 − 0.1321z−6 1+ 0.34319 z−2 + 0.60439 z−4 + 0.20407 z−6
求该系统在 0～π频率范围内的绝对幅频响应、相对幅度响应、相位响率响应及群迟延。
解 MATLAB 程序如下：
b＝[0.1321，0，0.3963，0，0.3963，0，0.1321]；
a＝[1，0，－0.34319，0，0.60439，0，－0.20407]；
N=1024;
[db，mag，pha，grd，w]＝freqz_m(b,a,N)；
n =1
M
∏ Cme jαm
=
K
m =1 N
= H(e jω ) e jϕ (ω)
∏ Dne jβn
n =1
其中，系统的幅度频响特性为
M
∏Cm
H(e jω )
=K
m=1 N
∏ Dn
n =1
系统的相位频响特性为
N
N
ϕ
(
ω
)亖∑
α
m亅∑
β
亄
n
ω(N亅M)
m=1
n =1
由公式可见，系统函数与频率响应有着密切的联系。适当地控制系统函数极点、零点的分
布，可以改变离散系统的频率响应特性：
(1)在原点(z＝0)处的零点或极点至单位圆的距离始终保持不变，其值|ejω|＝1，所以对幅度
响应不起作用。
(2)单位圆附近的零点对系统幅度响应的凹谷的位置及深度有明显的影响。
(3)单位圆内且靠近单位圆附近的极点对系统幅度响应的凸峰的位置及峰度有明显的影响。
2.系统的频率响应特性
实验 3 离散系统的频域分析
一、实验目的
(1)加深对离散系统的频率响应特性基本概念以及快速傅里叶变换(FFT)基本理论的理解。 (2)了解离散系统的零极点与频响特性之间关系。 (3)了解使用快速傅里叶变换(FFT)计算有限长序列和无限长序列信号频谱的方法。 (4)熟悉 MATLAB 中进行离散系统分析频响特性的常用子函数，掌握离散系统幅频响应和 相频响应的求解方法；掌握用 MATLAB 语言进行快速傅里叶变换时常用的子函数。
subplot(2，2，1)，plot(w/pi，mag)；grid %作绝对幅度频响图
axis([0，1，1.1*min(mag)，1.1*max(mag)])；
title(‘幅频特性(V)’)；
subplot(2，2，2)，plot(w/pi，pha)；grid%作相位频响图
subplot(2，2，1)，plot(w/pi，abs(h))；grid %作系统的绝对幅度频响图 axis([0，1，1.1*min(abs(h))，1.1*max(abs(h))])； title(‘幅频特性(V)’)； subplot(2，2，2)，plot(w/pi，angle(h))；grid %作系统的相位频响图 axis([0，1，1.1*min(angle(h))，1.1*max(angle(h))])； title(‘相频特性’)； subplot(2，2，3)，plot(w/pi，db)；grid %作系统的相对幅度频响图 axis([0，1，－100，5])； title(‘幅频特性( dB)’)； subplot(2，2，4)，zplane(b，a)；%作零极点分布图 title(‘零极点分布’)； 执行结果如图 3.3 所示。
求该系统在 0～π频率范围内的相对幅度频率响应与相位频率响应。
MATLAB 程序如下：
b＝［0.1321，0，－0.3963，0，0.3963，0，－0.1321］； a＝［1，0，0.34319，0，0.60439，0，0.20407］； freqz(b，a)； 以上程序采用了 freqz 不带输出向量的形式，直接出图。执行结果如图 3.1 所示。 由图 3.1 可见，该系统是一个 IIR 数字带通滤波器。其中幅频特性采用归一化的相对幅度值， 以分贝( dB)为单位。
图 3.2 例 3-2 所描述系统的幅度响应与相位响应
由图 10-2 可见，该系统是一个低通滤波器。其中，幅频特性采用归一化的绝对幅度值。
例 3-3
已知离散时间系统的系统函数为
H(z)
=
0.1− 0.4 z−1+ 0.4 z−2 − 0.1z−3 1+ 0.3 z−1+ 0.55 z−2 + 0.2 z−3
求该系统在 0～π频率范围内，归一化的绝对幅度频率响应、相对幅度频率响应、相位频率
响应及零极点分布图。
解 MATLAB 程序如下：
b＝［0.1，－0.4，0.4，－0.1］；
a＝［1，0.3，0.55，0.2］；
w1＝(0：499)*pi/500；
[h，w]＝freqz(b，a，w1)；
db＝20*log10(abs(h))； %求系统的相对幅频响应值
图 3.1 例 3-1 系统的幅度频率响应与相位频率响应 例 3-2 已知离散时间系统的系统函数，求该系统在 0～π频率范围内归一化的绝对幅度频率
响应与相位频率响应。 解 MATLAB 程序如下： b＝［0.2，0.1，0.3，0.1，0.2］； a＝［1，－1.1，1.5，－0.7，0.3］； w＝(0：499)*pi/500； %在 pi 的范围内取 500 个采样点 h＝freqz(b，a，w)； %求系统的频率响应 subplot(2，1，1)，plot(w/pi，abs(h))；grid %作系统的幅度频响图 axis(［0，1，1.1*min(abs(h))，1.1*max(abs(h))］)； ylabel(‘幅度’)； subplot(2，1，2)，plot(w/pi，angle(h))；grid %作系统的相位频响图 axis(［0，1，1.1*min(angle(h))，1.1*max(angle(h))］)； ylabel(‘相位’)；xlabel(‘以 pi 为单位的频率’)； 执行结果如图 3.2 所示。
二、实验涉及的 MATLAB 子函数
1. freqz
功能：用于求解离散时间系统的频率响应函数 H(ejω)。 调用格式： [h，w]＝freqz(b，a，n)；可得到数字滤波器的 n 点频率响应值，这 n 个点均匀地分布在 ［0，π］上，并将这 n 个频点的频率记录在 w 中，相应的频率响应值记录在 h 中。n 的缺省值 为 512。 [h，f]＝freqz(b，a，n，Fs)；用于对 H(ejω)在[0，Fs/2]上等间隔采样 n 点，每一个采样 点频率及相应频率响应值分别记录在 f 和 h 中。Fs(以 Hz 为单位) 可由用户指定任意指定。 h＝freqz(b，a，w)；用于对 H(e jω)在［0，2π］上进行采样，采样频率点由向量 w 指定。 h＝freqz(b，a，f，Fs)；用于对 H(e jω)在［0，Fs］上采样，采样频率点由向量 f 指定。 freqz(b，a，n)；用于在当前图形窗口中绘制幅度响应和相位响应特性曲线。 2. angle 功能：求相角。 调用格式： p＝angle(h)；用于求复矢量或复矩阵 H 的相角(以弧度为单位)，相角介于－π和＋π之间。 3. fft 功能：一维快速傅里叶变换(FFT)。 调用格式： y＝fft(x)；利用 FFT 算法计算矢量 x 的离散傅里叶变换，当 x 为矩阵时，y 为矩阵 x 每一 列的 FFT。当 x 的长度为 2 的幂次方时，则 fft 函数采用基 2 的 FFT 算法，否则采用稍慢的混 合基算法。 y＝fft(x，n)；采用 n 点 FFT。当 x 的长度小于 n 时，fft 函数在 x 的尾部补零，以构成 n 点数据；当 x 的长度大于 n 时，fft 函数会截断序列 x。当 x 为矩阵时，fft 函数按类似的方式 处理列长度。 4. ifft
w 中记录了对应［0，π］频率区域的 500 个频点的频率值。
下面举例说明 freqz_m 函数的使用方法。
例 3-4 已知离散时间系统的系统函数为
H(z)
=
0.1321+ 0.3963 z−2 + 0.3963 z−4 + 0.1321z−6 1− 0.34319 z−2 + 0.60439 z−4 − 0.20407 z−6
三、实验原理
1.离散系统频率响应的基本概念
已知稳定系统传递函数的零-极点增益(zpk)模型为
则系统的频响函数为
M
∏ (z− cm )
H(z)
=

K
m=1 N
∏ (z− dn )
n =1