数字信号处理第四章介绍

第四章线性时不变离散时间系统的频域分析

一、传输函数和频率响应

例4.1传输函数分析

Q4.1

clear;

M = input('Enter the filter length M: ');

w = 0:2*pi/1023:2*pi;

num = (1/M)*ones(1,M);

den = [1];

h = freqz(num, den, w);

subplot(2,1,1)

plot(w/pi,abs(h));grid

title('Magnitude Spectrum |H(e^{j\omega})|')

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,1,2)

plot(w/pi,angle(h));grid

title('Phase Spectrum arg[H(e^{j\omega})]')

xlabel('\omega /\pi');

ylabel('Phase in radians');

M=2 M=10 M=15

幅度谱为偶对称，相位谱为奇对称，这是一个低通滤波器。M越大，通带越窄且过渡带越陡峭。

Q4.2使用修改后的程序P3.1，计算并画出当w=[0,pi]时传输函数的因果线性时不变离散时间系统的频率响应。它表示哪种类型的滤波器？

w = 0:pi/511:pi;

num = [0.15 0 -0.15];

den = [1 -0.5 0.7];

如下图1这是一个带通滤波器。

图1 图2

Q4.3对下面的传输函数重做习题Q4.2：，式（4.36）和式（4.37）给出的两个滤波器之间的区别是什么？你将选择哪一个滤波器来滤波，为什么？

w = 0:pi/511:pi;

num = [0.15 0 -0.15];

den = [0.7 -0.5 1];

如上图2也是一个带通滤波器，这两个滤波器的幅度谱是一样的，相位谱不太一样，我会选择第一个带通滤波器，因为它的相位谱更加平滑，相位失真小。

Q4.4 使用MATLAB计算并画出当w=[0,pi]时因果线性时不变离散时间系统的群延迟。系统

的传输函数为。

clf;

w = 0:pi/511:pi;

num = [1 -1.2 1];

den = [1 -1.3 1.04 -0.222];

h= grpdelay(num,den,w);

plot(w/pi,h);

xlabel('w/pi');

ylabel('群延迟');

Q4.5 使用Q3.50中编写的程序，分别计算并画出式（4.36）和式（4.37）确定的两个滤波器的冲激响应中的前一百个样本。讨论你的结果。

clf;

num = [0.15 0 -0.15];

den = [0.7 -0.5 1];

L = input('输入样本数 L: ');

[g t] = impz(num,den,L);

stem(t,g);

title(['前 ',num2str(L),' 脉冲响应的样本']);

xlabel('时间序号 n');

ylabel('h[n]');

（4.36）式（4.37）式

由图可知：这些情节由impz给生成的因果的脉冲响应实现的H(z)。我们观察到Q4.3因果滤波器与H(z)在(4.36)稳定,这意味着H[n]是绝对可和,我们看到交替和指数衰减的脉冲响应。在另一方面,因果编档人员与H(z)在(4.37)极点以外的单位圆,是不稳定的。不足为奇的是,相应的h[n]上图显示与n指数增长。

Q4.6 传输函数的极零点图同样能分析线性时不变离散时间系统的性质。使用命令zplane 可以很容易地得到系统的极零点图。使用zplane分别生成式（4.36）和式（4.37）确定的两个滤波器的极零点图。讨论你的结果。

clf;

num = [0.15 0 -0.15];

den = [1 -0.5 0.7];

[z p k] = tf2zpk(num,den);

disp('Zeros:');

disp(z);

disp('Poles:');

disp(p);

input('Hit to continue...');

[sos k] = zp2sos(z,p,k)

input('Hit to continue...');

zplane(z,p);

式（4.36）

式（4.37）

由图可知：过滤器在(4.36)在单位圆和两极因此它的因果实现稳定；较低的图显示过滤器

(4.37)极点在单位圆外,其因果关系的实现是不稳定的。

二、传输函数的类型

例4.2滤波器

Q4.7

clf;

fc = 0.25;

n = [-6.5:1:6.5];

y = 2*fc*sinc(2*fc*n);k = n+6.5;

stem(k,y);title('N = 14');axis([0 13 -0.2 0.6]);

xlabel('Time index n');

ylabel('Amplitude');

grid;

图1 图2

如图1低通有限冲激滤波器的长度为14，决定滤波器长度的语句为n = [-6.5:1:6.5]，而控制截止频率的参数是fc = 0.25。

Q4.8

fc = 0.45;

n = [-9.5:1:9.5];

y = 2*fc*sinc(2*fc*n);k = n+9.5;

stem(k,y);title('N = 20');axis([0 19 -0.2 0.7]);

xlabel('Time index n');

ylabel('Amplitude');

grid;

修改参数fc和n，得到如上图2，可知低通有限冲激滤波器的长度变为20.

Q4.9

clf;

fc = 0.65;

n = [-7.0:1:7.0];

y = 2*fc*sinc(2*fc*n);k = n+7.0;

stem(k,y);title('N = 14');axis([0 14 -0.4 1.4]);

xlabel('Time index n');

ylabel('Amplitude');

grid;

Q4.10

clear;

N = input('Enter the filter time shift N: ');

No2 = N/2;fc = 0.25;n = [-No2:1:No2];

y = 2*fc*sinc(2*fc*n);w = 0:pi/511:pi;h = freqz(y, [1], w);

plot(w/pi,abs(h));

grid;

title(strcat('|H(e^{j\omega})|, N=',num2str(N)));

xlabel('\omega /\pi');

ylabel('Amplitude');

上图依次分别为N=5,10,30,100的四幅图，从这四幅图可以看出随着阶数N的增大，低通滤波器的过渡带越来越窄，阻带衰减越来越快，滤波器越来越接近理想低通滤波器。

Q4.11

clf;

M = 2;

num = ones(1,M)/M;

[g,w] = gain(num,1);

plot(w/pi,g);grid

axis([0 1 -50 0.5])

xlabel('\omega /\pi');ylabel('Gain in dB');

title(['M = ',num2str(M)])

可以验证3dB截止频率在π/2处。

Q4.12

clear;

K = input('Enter the number of sections K: ');

Hz = [1];

for i=1:K;

Hz = conv(Hz,[1 1]);

end;

Hz = (0.5)^K * Hz;

[g,w] = gain(Hz,1);

ThreedB = -3*ones(1,length(g));

t1 = 2*acos((0.5)^(1/(2*K)))*ones(1,512)/pi;

t2 = -50:50.5/511:0.5;

plot(w/pi,g,w/pi,ThreedB,t1,t2);grid;

axis([0 1 -50 0.5])

xlabel('\omega /\pi');ylabel('Gain in dB');

title(['K = ',num2str(K),'; Theoretical \omega_{c} = ',num2str(t1(1))]);

Q4.13

clear;

M = input('Enter the filter length M: ');

n = 0:M-1;

num = (-1).^n .* ones(1,M)/M;

[g,w] = gain(num,1);

plot(w/pi,g);grid;

axis([0 1 -50 0.5]);

xlabel('\omega /\pi');ylabel('Gain in dB');

title(['M = ', num2str(M)]);

其3dB截止频率约为0.82pi

Q4.14 设计一个在0.45pi处具有3dB截止频率wc的一阶无限冲激响应低通滤波器和一阶无限冲激响应高通滤波器。用MATLAB计算并画出它们的增益响应，验证设计的滤波器是否满足指标。用MATLAB证明两个滤波器是全通互补和功率互补的。

Q4.15 级联10个式（4.15）所示一阶无限冲激响应低通滤波器，设计一个在0.3pi处具有3dB截止频率wc的无限冲激响应低通滤波器。把它与一个具有相同截止频率的一阶无限冲激响应低通滤波器的增益响应作比较。

Q4.16 设计一个中心频率wo在0.61pi处、3dB带宽为0.51pi的二阶无限冲激响应带通滤波器。由于式（4.20）是α的二次方程，为了产生相同的3dB带宽，参数α将有两个数值，得到的传输函数HBP（z）也会有两个不同的表达式。使用函数zplane可产生两个传输函数的极零点图，从中可以选择一个稳定的传输函数。用MATLAB计算并画出你所设计的滤波器的增益响应，并验证它确实满足给定的条件。用设计的稳定无限冲激响应带通滤波器的传输函数的参数α和β，生成一个二阶无限冲激响应带阻滤波器的传输函数HBS（z）。用MATLAB 证明HBP(z)和HBS(z)都是全通互补和功率互补的。

Q4.17 用MATLAB计算并画出一个梳状滤波器的幅度响应，该梳状滤波器是在L取不同值的情况下，由式（4.40）给出的原型有限冲激响应低通滤波器得到的。证明新滤波器的幅度响

应在处有L个极小值，在处有L个极大值，

k=0,1,2,3......，L-1.

Q4.18 用MATLAB计算并画出一个梳状滤波器的幅度响应，该梳状滤波器是在L取不同值的情况下，由式（4.42）在M=2时给出的原型有限冲激响应低通滤波器得到的。确定这种梳状滤波器冲激响应的极大值和极小值的位置。

从这些情节我们观察,梳状滤波器极距为1kπ/L,山峰为(2k+1)π/L.

Q4.19

clf;

b = [1 -8.5 30.5 -63];

num1 = [b 81 fliplr(b)];

num2 = [b 81 81 fliplr(b)];

num3 = [b 0 -fliplr(b)];

num4 = [b 81 -81 -fliplr(b)];

n1 = 0:length(num1)-1;

n2 = 0:length(num2)-1;

subplot(2,2,1); stem(n1,num1);

xlabel('Time index n');ylabel('Amplitude'); grid; title('Type 1 FIR Filter');

subplot(2,2,2); stem(n2,num2);

xlabel('Time index n');ylabel('Amplitude'); grid; title('Type 2 FIR Filter');

subplot(2,2,3); stem(n1,num3);

xlabel('Time index n');ylabel('Amplitude'); grid; title('Type 3 FIR Filter');

subplot(2,2,4); stem(n2,num4);

xlabel('Time index n');ylabel('Amplitude'); grid; title('Type 4 FIR Filter');

pause

subplot(2,2,1); zplane(num1,1);

title('Type 1 FIR Filter');

subplot(2,2,2); zplane(num2,1);

title('Type 2 FIR Filter');

subplot(2,2,3); zplane(num3,1);

title('Type 3 FIR Filter');

subplot(2,2,4); zplane(num4,1);

title('Type 4 FIR Filter');

disp('Zeros of Type 1 FIR Filter are');

disp(roots(num1));

disp('Zeros of Type 2 FIR Filter are');

disp(roots(num2));

disp('Zeros of Type 3 FIR Filter are');

disp(roots(num3));

disp('Zeros of Type 4 FIR Filter are');

disp(roots(num4));

1型有限冲激响应滤波器的零点是Zeros of Type 1 FIR Filter are 2.9744

2.0888

0.9790 + 1.4110i

0.9790 - 1.4110i

0.3319 + 0.4784i

0.3319 - 0.4784i

0.4787

0.3362

2型有限冲激响应滤波器的零点是Zeros of Type 2 FIR Filter are 3.7585 + 1.5147i

3.7585 - 1.5147i

0.6733 + 2.6623i

0.6733 - 2.6623i

-1.0000

0.0893 + 0.3530i

0.0893 - 0.3530i

0.2289 + 0.0922i

0.2289 - 0.0922i

3型有限冲激响应滤波器的零点是Zeros of Type 3 FIR Filter are 4.7627

1.6279 + 3.0565i

1.6279 - 3.0565i

-1.0000

1.0000

0.1357 + 0.2549i

0.1357 - 0.2549i

0.2100

4型有限冲激响应滤波器的零点是Zeros of Type 4 FIR Filter are 3.4139

1.6541 + 1.5813i

1.6541 - 1.5813i

-0.0733 + 0.9973i

-0.0733 - 0.9973i

1.0000

0.3159 + 0.3020i

0.3159 - 0.3020i

0.2929

Q4.20

用b=[1.5 -3.25 5.25 -4]代替原来的b，所得的结果为

1型有限冲激响应滤波器的零点是

Zeros of Type 1 FIR Filter are

2.3273 + 2.0140i

2.3273 - 2.0140i

-1.2659 + 2.0135i

-1.2659 - 2.0135i

-0.2238 + 0.3559i

-0.2238 - 0.3559i

0.2457 + 0.2126i

0.2457 - 0.2126i

2型有限冲激响应滤波器的零点是

Zeros of Type 2 FIR Filter are

2.5270 + 2.0392i

2.5270 - 2.0392i

-1.0101 + 2.1930i

-1.0101 - 2.1930i

-1.0000

-0.1733 + 0.3762i

-0.1733 - 0.3762i

0.2397 + 0.1934i

0.2397 - 0.1934i

3型有限冲激响应滤波器的零点是

Zeros of Type 3 FIR Filter are

-1.0000

0.2602 + 1.2263i

0.2602 - 1.2263i

1.0000

0.6576 + 0.7534i

0.6576 - 0.7534i

0.1655 + 0.7803i

0.1655 - 0.7803i

4型有限冲激响应滤波器的零点是

Zeros of Type 4 FIR Filter are

2.0841 + 2.0565i

2.0841 - 2.0565i

-1.5032 + 1.9960i

-1.5032 - 1.9960i

1.0000

-0.2408 + 0.3197i

-0.2408 - 0.3197i

0.2431 + 0.2399i

0.2431 - 0.2399i

Q4.21 用MATLAB 确定如下传输函数是否是有界实函数：它若不是一个有界实函数，求一个与有着相同幅度的有界实函数。

由下图可知：H1(z)不是有界实函数。

故H2(z)为

Q4.22 用MATLAB 确定如下传输函数是否是有界实函数：它若不是一个有界实函数，求一个与有着相同幅度的有界实函数。

使用zplane我们观察到的G1(z)在单位圆,因此传递函数是稳定的。

Q4.23 用MATLAB产生如下两个因果系统传输函数的极零点图：，

，研究生成的极零点图，你可以推断它们的稳定性么？

用Q4.6的程序做H1(z)

Q4.24 用程序P4.4检测Q4,23中两个传输函数的稳定性。这两个传输函数哪一个是稳定的？% Program P4_4

clf;

den = input('分母系数 = ');

ki = poly2rc(den);

disp('稳定性测试参数是');

disp(ki);

由此我们可以总结出H1(z)稳定，H2(z)不稳定

Q4.25 用程序P4.4确定下面这个多项式的所有根是否都在单位圆内：

由此看出，都在单位圆内。Q4.26 用程序P4.4确定下面这个多项式的所有根是否都在单位圆内：

由此看出，都在单位圆内。