2024年3月8日发(作者：联想威6 2021款笔记本怎么样)

精品文档

南京邮电大学

实 验 报 告

实验名称：离散时间信号与系统的时、频域表示

离散傅立叶变换和z变换

数字滤波器的频域分析和实现

数字滤波器的设计

课程名称

数字信号处理A(双语)

班级学号

B13011025

姓 名 陈志豪

开课时间 2015/2016学年，第1学期

精品文档

精品文档

实验名称：离散时间信号与系统的时、频域表示

实验目的和任务：

熟悉Matlab基本命令，理解和掌握离散时间信号与系统的时、频域表示及简单应用。在Matlab环境中，按照要求产生序列，对序列进行基本运算；对简单离散时间系统进行仿真，计算线性时不变（LTI）系统的冲激响应和卷积输出；计算和观察序列的离散时间傅立叶变换（DTFT）幅度谱和相位谱。

实验内容：

基本序列产生和运算： Q1.1～1.3，Q1.23，Q1.30～1.33

离散时间系统仿真： Q2.1～2.3

LTI系统：Q2.19，Q2.21，Q2.28

DTFT：Q3.1，Q3.2，Q3.4

实验过程与结果分析：

Q1.1运行程序P1.1，以产生单位样本序列u[n]并显示它。

clf;

n = -10:20;

u = [zeros(1,10) 1 zeros(1,20)];

stem(n,u);

xlabel('Time index n');

ylabel('Amplitude');

title('Unit Sample Sequence');

axis([-10 20 0 1.2]);

精品文档

精品文档

Q1.2 命令clf，axis，title，xlabel和ylabel 命令的作用是什么？

答：clf命令的作用：清除图形窗口上的图形；

axis命令的作用：设置坐标轴的范围和显示方式；

title命令的作用：给当前图片命名；

xlabel命令的作用：添加x坐标标注；

ylabel c命令的作用：添加y坐标标注；

Q1.3修改程序P1.1，以产生带有延时11个样本的延迟单位样本序列ud[n]。运行修改的程序并显示产生的序列。

clf;

n = -10:20;

u = [zeros(1,21) 1 zeros(1,9)];

stem(n,u);

xlabel('Time index n');

ylabel('Amplitude');

title('Unit Sample Sequence');

axis([-10 20 0 1.2]);

Q1.23修改上述程序，以产生长度为50、频率为0.08、振幅为2.5、相移为90度的一个正弦序列并显示它。该序列的周期是多少？

n = 0:50;

精品文档

精品文档

f = 0.08;

phase = 90;

A = 2.5;

arg = 2*pi*f*n - phase;

x = A*cos(arg);

clf;

stem(n,x);

axis([0 50 -3 3]);

grid;

title('Sinusoidal Sequence');

xlabel('Time index n');

ylabel('Amplitude');

axis;

答：周期为：T=2=11==22.5。

f0.08Q1.30未污染的信号s[n]是什么样的形式？加性噪声d[n]是什么样的形式？

答：未污染的信号：s[n]=2×0.9。

加性噪声d[n]是均匀分布在-04到+0.4之间的随机序列。

Q1.31使用语句x=s+d能产生被噪声污染的信号吗？若不能，为什么？

答：不能，因为d是列向量，s是行向量。

Q1.32信号x1，x2和x3与信号x之间的关系是什么？

答：X1是x的延时一个单位，x2和x相等，x3超前于x一个单位。

精品文档

nn

精品文档

Q1.33legend命令的作用是什么？

答：产生图例说明。

Q2.1对M = 2 ，运行上述程序，生成输入 x[n] = s1[n]+s2[n]的输出信号。输出x[n]的哪个分量被该离散时间系统抑制？

答：输入 x[n] 被该离散时间系统抑制的分量为Signal2的高频分量。

Q2.2若线性时不变系统由y[n] = 0.5(x[n]+x[n–1])变成y[n] = 0.5(x[n]–x[n–1])，输入 x[n] = s1[n]+s2[n]的影响是什么？

n = 0:100;

s1 = cos(2*pi*0.05*n);

s2 = cos(2*pi*0.47*n); x = s1+s2;

M = input('Desired length of the filter = ');

num = (-1).^[0:M-1];

y = filter(num,1,x)/M;

clf;

subplot(2,2,1);

plot(n, s1);

axis([0, 100, -2, 2]);

xlabel('Time index n');

ylabel('Amplitude');

title('Signal #1');

subplot(2,2,2);

plot(n, s2);

精品文档

对

精品文档

axis([0, 100, -2, 2]);

xlabel('Time index n');

ylabel('Amplitude');

title('Signal #2');

subplot(2,2,3);

plot(n, x);

axis([0, 100, -2, 2]);

xlabel('Time index n');

ylabel('Amplitude');

title('Input Signal');

subplot(2,2,4);

plot(n, y);

axis([0, 100, -2, 2]);

xlabel('Time index n');

ylabel('Amplitude');

title('Output Signal');

axis;

答：对于输入的影响是- 该系统是一个高通滤波器，它通过高频率的输入分量S2，而不是低频的输入分量S1。

Q2.3对滤波器长度M和正弦信号s1[n]和s2[n]的频率取其他值，运行程序P2.1，算出结果。

M=3，f1=0.1，f2=0.2

精品文档

精品文档

M=8，f1=0.25，f2=0.5

Q2.19运行 P2_5，生成式（2.15）所给离散时间系统的冲激响应。

精品文档

精品文档

Q2.21利用filter命令编写一个MATLAB程序，生成式（2.17）给出的因果线性时不变系统的冲激响应，计算并画出前40个样本。把你的结果和习题Q2.20中得到的结果相比较。

clf;

N = 40;

num = [0.9 -0.45 0.35 0.002];

den = [1.0 0.71 -0.46 -0.62];

x = [1 zeros(1,N-1)];

y = filter(num,den,x);

stem(y);

xlabel('Time index n');

ylabel('Amplitude');

title('Impulse Response');

grid;

程序产生的40个样本如下所示:

精品文档

精品文档

Q2.28运行程序P2.7，对序列h[n]和x[n]求卷积，生成y[n]，并用FIR滤波器h[n]对输入x[n]滤波，求得y1[n]。y[n]和y1[n]有差别吗？为什么要使用对x[n]补零后得到的x1[n]作为输入来产生y1[n]？

精品文档

精品文档

答：y[n] 和 y1[n] 的差别为 - 没有差别

将x[n]补零后得到 x1[n]作为输入,产生y1[n]的原因是 – 对于长度N1和N2的两个序列，转化率返回得到的序列长度N1+N2-1。与此相反，滤波器接收一个输入信号和系统规范，返回的结果是相同的长度作为输入信号。因此，为了从转化率和滤波器得到直接比较的结果，有必要供应滤波器的输入已经填充为长度L(x)+L(h)-1。

Q3.1在程序P3.1中，计算离散时间傅里叶变换的原始序列是什么？MATLAB命令pause的作用是什么？

2z1答：离散时间傅里叶变换的原始序列：；

110.6zpause 命令的作用：暂停程序，直至用户按任意键，程序才可以开始。

Q3.2运行程序 P3.1，求离散时间傅里叶变换的实部、虚部以及幅度和相位谱。离散时间傅里叶变换是的周期函数吗？若是，周期是多少？描述这四个图表示的对称性。

精品文档

精品文档

答：DTFT 是关于 的周期函数么？是周期函数

周期是 - 2π

四个图形的对称性为：实部是偶对称，虚部是奇对称，幅度谱是偶对称相位谱是奇对称。

Q3.4修改程序 P3_1 重做Q3.2的程序如下:

clf;

w = -4*pi:8*pi/511:4*pi;

num = [1 3 5 7 9 11 13 15 17];den = [1];

h = freqz(num, den, w);

subplot(2,1,1)

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

title('Real part of H(e^{jomega})')

xlabel('omega /pi');

ylabel('Amplitude');

subplot(2,1,2)

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

title('Imaginary part of H(e^{jomega})')

xlabel('omega /pi');

ylabel('Amplitude');

pause

subplot(2,1,1)

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

精品文档

精品文档

title('Magnitude Spectrum |H(e^{jomega})|')

xlabel('omega /pi');

ylabel('Amplitude');

subplot(2,1,2)

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

title('Phase Spectrum arg[H(e^{jomega})]')

xlabel('omega /pi');

ylabel('Phase in radians');

修改程序后的运行结果为:

精品文档

精品文档

答：DTFT 是关于 的周期函数么？是周期函数

周期是 - 2π

相位谱中跳变的原因是 - 对相位进行归一化。

实验名称：离散傅立叶变换和z变换

实验目的和任务：

掌握离散傅立叶变换（DFT）及逆变换（IDFT）、z变换及逆变换的计算和分析。利用Matlab语言，完成DFT和IDFT的计算及常用性质的验证，用DFT实现线性卷积，实现z变换的零极点分析，求有理逆z变换。

实验内容：

DFT和IDFT计算： Q3.23～3.24

DFT的性质： Q3.26~3.29，Q3.36，Q3.38，Q3.40

z变换分析：Q3.46～3.48

逆z变换：Q3.50

实验过程与结果分析：…参见实验一格式…

Q3.23编写一个MATLAB程序，计算并画出长度为N的L点离散傅里叶变换X[k]的值，

其中L≥N，然后计算并画出L点离散傅里叶逆变换X[k]。对不同长度N和不同的离散傅里叶变换长度L，运行程序。讨论你的结果。

clf;

N=200;

L=256;

nn=[0:N-1];

精品文档

精品文档

kk=[0:L-1];

xR=[0.1*(1:100) zeros(1,N-100)];

xI=[zeros(1,N)];

x=xR+i*xI;

XF=fft(x,L);

subplot(3,2,1);

grid;

plot(nn,xR);

grid;

title('Re{x[n]}');

xlabel('Time index n');

ylabel('Amplitude');

subplot(3,2,2);

plot(nn,xI);

grid;

title('Im{x[n]}');

xlabel('Time index n');

ylabel('Amplitude');

subplot(3,2,3);

plot(kk,real(XF));

grid;

title('Re{x[n]}');

xlabel('Frequency index k');

ylabel('Amplitude');

subplot(3,2,4);

plot(kk,imag(XF));

grid;

title('Im{x[n]}');

xlabel('Frequency index k');

ylabel('Amplitude');

xx=ifft(XF,L);

subplot(3,2,5);

plot(kk,real(xx));

grid;

title('Real part of IDFT{x[n]}');

xlabel('Time index n');

ylabel('Amplitude');

subplot(3,2,6);

plot(kk,imag(xx));

grid;

title('Imag part of IDFT{x[n]}');

xlabel('Time index n');

ylabel('Amplitude');

精品文档

精品文档

Q3.26在函数circshift中，命令rem的作用是什么？

答：Rem(x,y)是用y对x求余数函数。

Q3.27 解释函数circshift怎样实现圆周移位运算。

答：在输入序列x由M的位置开始被循环移位，如果M>0，则circlshift删除从矢量x最左边开始的M个元素和它们附加在右侧的剩余元素，以获得循环移位序列。如果M<0，则circlshift首先通过x的长度来弥补M，即序列x最右边的长度的M样本从x中删除和附加在其余的M个样本的右侧，以获得循环移位序列。

Q3.28在函数circonv中，运算符~=的作用是什么？

答：“~=”表示不等于。

Q3.29解释函数circonv怎样实现圆周卷积运算。

答：输入时两个长度都为L的向量x1和x2，它是非常有用的定期延长X2的函数。让x2p成为x2延长无限长的周期的序列。从概念上讲，在定点时间上通过时序交换后的x2p的长度L交换x2p序列和x2tr等于1的元素。然后元素1至L的输出向量y是通过取x1和获得的长度为L的sh矢量之间的内积得到通过循环右移的时间反转向量x2tr。对于输出样本Y[n]的1≤N≤L时，右循环移位的量为n-1个位置上。

Q3.30通过加入合适的注释语句和程序语句，修改程序P3.7，对程序生成的图形中的两个轴加标记。哪一个参数决定时移量？若时移量大于序列长度，将会发生什么？

function y = circshift(x,M)

if abs(M) > length(x)

M = rem(M,length(x));

end

if M < 0

M = M + length(x);

end

y = [x(M+1:length(x)) x(1:M)];

精品文档

精品文档

clf;

M = 6;

a = [0 1 2 3 4 5 6 7 8 9];

b = circshift(a,M);

L = length(a)-1;

n = 0:L;

subplot(2,1,1);

stem(n,a);axis([0,L,min(a),max(a)]);

title('Original Sequence');

xlabel('time index k');

ylabel('a[n]');

subplot(2,1,2);

stem(n,b);axis([0,L,min(a),max(a)]);

title(['Sequence Obtained by Circularly Shifting

by',num2str(M),'Samples']);

xlabel('time index n');

ylabel('b[n]');

D决定时移量，左移M位。

Q3.31运行修改后的程序并验证圆周时移运算。

修改后的程序如下：

Q3.32通过加入合适的注释语句和程序语句，修改程序P3.8，对程序生成的图形中的两个轴加标记。时移量是多少？

clf;

精品文档

精品文档

x = [0 2 4 6 8 10 12 14 16];

N = length(x)-1; n = 0:N;

y = circshift(x,5);

XF = fft(x);

YF = fft(y);

subplot(2,2,1)

stem(n,abs(XF));grid

title('Magnitude of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('|X[k]|');

subplot(2,2,2)

stem(n,abs(YF));grid

title('Magnitude of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('|Y[k]|');

subplot(2,2,3)

stem(n,angle(XF));grid

title('Phase of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('arg(X[k])');

subplot(2,2,4)

stem(n,angle(YF));grid

title('Phase of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('arg(Y[k])');

时移量是5。

Q3.33运行修改后的程序并验证离散傅里叶变换的圆周时移性质。

精品文档

精品文档

Q3.34选取两个不同的时移量，重做习题Q3.33。

（1）修改后的程序如下：

clf;

x = [0 2 4 6 8 10 12 14 16];

N = length(x)-1; n = 0:N;

y = circshift(x,3);

XF = fft(x);

YF = fft(y);

subplot(2,2,1)

stem(n,abs(XF));grid

title('Magnitude of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('|X[k]|');

subplot(2,2,2)

stem(n,abs(YF));grid

title('Magnitude of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('|Y[k]|');

subplot(2,2,3)

stem(n,angle(XF));grid

title('Phase of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('arg(X[k])');

subplot(2,2,4)

stem(n,angle(YF));grid

title('Phase of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('arg(Y[k])');

精品文档

精品文档

（2）修改后的程序如下：

clf;

x = [0 2 4 6 8 10 12 14 16];

N = length(x)-1; n = 0:N;

y = circshift(x,6);

XF = fft(x);

YF = fft(y);

subplot(2,2,1)

stem(n,abs(XF));grid

title('Magnitude of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('|X[k]|');

subplot(2,2,2)

stem(n,abs(YF));grid

title('Magnitude of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('|Y[k]|');

subplot(2,2,3)

stem(n,angle(XF));grid

title('Phase of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('arg(X[k])');

subplot(2,2,4)

stem(n,angle(YF));grid

title('Phase of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('arg(Y[k])');

Q3.35选取两个不同长度的序列，重做习题Q3.33。

精品文档

精品文档

修改后的程序如下：

clf;

x = [2 4 6 8 10 12 14 16];

N = length(x)-1; n = 0:N;

y = circshift(x,5);

XF = fft(x);

YF = fft(y);

subplot(2,2,1)

stem(n,abs(XF));grid

title('Magnitude of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('|X[k]|');

subplot(2,2,2)

stem(n,abs(YF));grid

title('Magnitude of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('|Y[k]|');

subplot(2,2,3)

stem(n,angle(XF));grid

title('Phase of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('arg(X[k])');

subplot(2,2,4)

stem(n,angle(YF));grid

title('Phase of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('arg(Y[k])');

修改后的程序如下：

精品文档

精品文档

clf;

x = [0 2 4 6 8 10 12];

N = length(x)-1; n = 0:N;

y = circshift(x,5);

XF = fft(x);

YF = fft(y);

subplot(2,2,1)

stem(n,abs(XF));grid

title('Magnitude of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('|X[k]|');

subplot(2,2,2)

stem(n,abs(YF));grid

title('Magnitude of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('|Y[k]|');

subplot(2,2,3)

stem(n,angle(XF));grid

title('Phase of DFT of Original Sequence');

xlabel('Frequency index k'); ylabel('arg(X[k])');

subplot(2,2,4)

stem(n,angle(YF));grid

title('Phase of DFT of Circularly Shifted Sequence');

xlabel('Frequency index k'); ylabel('arg(Y[k])');

Q3.36运行程序P3.9并验证离散傅里叶变换的圆周卷积性质

精品文档

精品文档

g1 = [1 2 3 4 5 6]; g2 = [1 -2 3 3 -2 1];

ycir = circonv(g1,g2);

disp('Result of circular convolution = ');disp(ycir)

G1 = fft(g1); G2 = fft(g2);

yc = real(ifft(G1.*G2));

disp('Result of IDFT of the DFT products = ');disp(yc)

结果：

x1 =

1 2 3 4 5 6

x2 =

1 -2 3 3 -2 1

Result of circular convolution =

12 28 14 0 16 14

Result of IDFT of the DFT products =

12 28 14 0 16 14

Q3.37选取另外两组等长序列重做习题Q3.36。

g1 = [2 3 7 1 6 -5]; g2 = [-4 6 7 0 2 9];

ycir = circonv(g1,g2);

disp('Result of circular convolution = ');disp(ycir)

G1 = fft(g1); G2 = fft(g2);

yc = real(ifft(G1.*G2));

disp('Result of IDFT of the DFT products = ');disp(yc)

x1 =

2 3 7 1 6 -5

x2 =

-4 6 7 0 2 9

Result of circular convolution =

45 30 25 103 -10 87

Result of IDFT of the DFT products =

45.0000 30.0000 25.0000 103.0000 -10.0000 87.0000

Q3.38运行程序P3.10并验证线性卷积可通过圆周卷积得到。

g1 = [1 2 3 4 5];g2 = [2 2 0 1 1];

g1e = [g1 zeros(1,length(g2)-1)];

g2e = [g2 zeros(1,length(g1)-1)];

ylin = circonv(g1e,g2e);

disp('Linear convolution via circular convolution = ');disp(ylin);

y = conv(g1, g2);

disp('Direct linear convolution = ');disp(y)

x1 =

1 2 3 4 5 0 0 0 0

x2 =

精品文档

精品文档

2 2 0 1 1 0 0 0 0

Linear convolution via circular convolution =

2 6 10 15 21 15 7 9 5

Direct linear convolution =

2 6 10 15 21 15 7 9 5

Q3.39选取两组长的不等的序列重做习题Q3.38。

g1 = [1 2 3 4 5];g2 = [2 2 0 1 1 5];

g1e = [g1 zeros(1,length(g2)-1)];

g2e = [g2 zeros(1,length(g1)-1)];

ylin = circonv(g1e,g2e);

disp('Linear convolution via circular convolution = ');disp(ylin);

y = conv(g1, g2);

disp('Direct linear convolution = ');disp(y)

x1 =

1 2 3 4 5 0 0 0 0 0

x2 =

2 2 0 1 1 5 0 0 0 0

Linear convolution via circular convolution =

2 6 10 15 21 20 17 24 25 25

Direct linear convolution =

2 6 10 15 21 20 17 24 25 25

Q3.40编写一个MATLAB程序，对两个序列做离散傅里叶变换，以生成他们的线性卷积。用此程序验证习题Q3.38和习题Q3.39的结果。

验证Q3.38:

g1=[1 2 3 4 5];

g2=[2 2 0 1 1];

g1e = [g1 zeros(1,length(g2)-1)];

g2e = [g2 zeros(1,length(g1)-1)];

G1EF=fft(g1e);

G2EF=fft(g2e);

ylin=real(ifft(G1EF.*G2EF));

disp('Linear convolution via DFT =');

disp(ylin);

精品文档

精品文档

Linear convolution via DFT =

2.0000 6.0000 10.0000 15.0000 21.0000 15.0000 7.0000 9.0000

5.0000

验证Q3.39:

g1=[1 2 3 4 5];

g2=[2 2 0 1 1 5];

g1e = [g1 zeros(1,length(g2)-1)];

g2e = [g2 zeros(1,length(g1)-1)];

G1EF=fft(g1e);

G2EF=fft(g2e);

ylin=real(ifft(G1EF.*G2EF));

disp('Linear convolution via DFT =');

disp(ylin);

Linear convolution via DFT =

2.0000 6.0000 10.0000 15.0000 21.0000 20.0000

17.0000 24.0000 25.0000 25.0000

25z19z25z33z4Q3.46使用程序P3.1在单位圆生求下面的z变换：G(z)

1234545z2zzzclf;

w = -4*pi:8*pi/511:4*pi;

num = [2 5 9 5 3];den = [5 45 2 1 1];

h = freqz(num, den, w);

subplot(2,1,1)

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

title('Real part of H(e^{jomega})')

xlabel('omega /pi');

ylabel('Amplitude');

subplot(2,1,2)

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

title('Imaginary part of H(e^{jomega})')

xlabel('omega /pi');

ylabel('Amplitude');

pause

subplot(2,1,1)

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

title('Magnitude Spectrum |H(e^{jomega})|')

xlabel('omega /pi');

ylabel('Amplitude');

subplot(2,1,2)

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

title('Phase Spectrum arg[H(e^{jomega})]')

精品文档

精品文档

xlabel('omega /pi');

ylabel('Phase in radians');

精品文档

精品文档

Q3.47编写一个MATLAB程序，计算并显示零点和极点，计算并显示其因式形式，并产生以z的两个多项式之比的形式表示的z变换的零点图。使用该程序，分析式（3.32）的z变换G(z)。

1clf;

num=[2 5 9 5 3];

den=[5 45 2 1 1];

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

disp('Zeros:');

disp(z);

disp('Poles:');

disp(p);

input('Hit ');

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

input('Hit ');

zplane(z,p);

Zeros:

-1.0000 + 1.4142i

-1.0000 - 1.4142i

-0.2500 + 0.6614i

-0.2500 - 0.6614i

Poles:

-8.9576

-0.2718

0.1147 + 0.2627i

0.1147 - 0.2627i

Hit

sos =

1.0000 2.0000 3.0000 1.0000

1.0000 0.5000 0.5000 1.0000

k =

精品文档

9.2293

-0.2293

2.4344

0.0822

精品文档

0.4000

Hit

Q3.48通过习题Q3.47产生的极零点图，求出G(z)的收敛域的数目。清楚地显示所有的收敛域。由极零点图，说明离散时间傅里叶变换是否存在。

答：R1：|z|<0.2718（左边序列，不稳定）

R2：0.2718<|z|<0.2866（双边序列，不稳定）

R3：0.2866<|z|<8.9576（左边序列，稳定）

R4：|z|>8.9576（右边序列，不稳定）

不能从极零点图肯定地说DTFT是否存在，因为其收敛域一定要指定，当收敛域在上述R3内所获得的序列确实证明了DTFT的存在，它是一个具有双面冲击响应的稳定系统。

Q3.50编写一个MATLAB程序，计算一个有理逆z变换的前L个样本，其中L的值由用户通过命令input提供。用该程序计算并画出式（3.32）中G(z)的逆变换的前50个样本。使用命令stem画出由逆变换产生的序列。

clf;

num=[2 5 9 5 3];

den=[5 45 2 1 1];

精品文档

精品文档

L=input('Enter the number of samples L:');

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

stem(t,g);

title(['First',num2str(L),'samples of impulse response']);

xlabel('Time Index n');

ylabel('h[n]');

实验名称：数字滤波器的频域分析和实现

实验目的：

掌握滤波器的传输函数和频率响应的关系，能够从频率响应和零极点模式分析滤波器特性。掌握滤波器的常用结构。

实验任务：

求滤波器的幅度响应和相位响应，观察对称性，判断滤波器类型，判断稳定性。验证FIR线性相位滤波器的特点。实现数字滤波器的直接型、级联型和并联型结构。

实验内容：

传输函数和频率响应、滤波器稳定性：Q4.1～4.3，Q4.5，Q4.6，Q4.19

线性相位滤波器：Q4.19

数字滤波器结构：Q6.1，Q6.3，Q6.5

数字滤波器仿真：Q8.1，Q 8.3，Q 8.5，Q 8.9，Q 8.10，Q 8.14

实验过程与结果分析：…参见实验一格式…

Q4.1修改程序P3.1，取三个不同的M值，当02时计算并画出式（2.13）所示滑动平均滤波器的幅度和相位谱。证明由幅度和相位谱表现出的对称类型。它表示了哪种类型的滤波器？你现在能解释Q2.1的结果吗？

答：

精品文档

精品文档

clf;

w = 0:8*pi/511:2*pi;

M=input('滤波器所需长度=');

num = ones(1,M);

y = freqz(num, 1, w)/M;

subplot(2,1,1)

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

title('Magnitude Spectrum |Y(e^{jomega})|')

xlabel('omega /pi');

ylabel('Amplitude');

subplot(2,1,2)

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

title('Phase Spectrum arg[Y(e^{jomega})]')

xlabel('omega /pi');

ylabel('Phase in radians');

Ｍ＝２

Ｍ＝５

精品文档

精品文档

Ｍ＝８

它表示了Ⅱ型滤波器。

Q4.2使用修改后的程序P3.1，计算并画出当0时传输函数精品文档

精品文档

0.15(1z2)

H(z)10.5z10.7z2的因果线性时不变离散时间系统的频率响应。它表示哪种类型的滤波器？

答：clf;

w = 0:8*pi/511:pi;

M=input('滤波器所需长度=');

num =[0.15 0 -0.15];

den=[1 -0.5 0.7];

y = freqz(num, den, w)/M;

subplot(2,1,1)

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

title('Magnitude Spectrum |Y(e^{jomega})|')

xlabel('omega /pi');

ylabel('Amplitude');

subplot(2,1,2)

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

title('Phase Spectrum arg[Y(e^{jomega})]')

xlabel('omega /pi');

ylabel('Phase in radians');

滤波器所需长度=1

它表示带通滤波器。

Q4.3对下面的传输函数重做习题Q4.2：

精品文档

精品文档

0.15(1z2)

G(z)0.70.5z1z2式（4.36）和式（4.37）给出的两个滤波器之间的区别是什么？你将选择哪一个滤波器来滤波，为什么？

clf;

w = 0:8*pi/511:pi;

M=input('滤波器所需长度=');

num =[0.15 0 -0.15];

den=[0.7 -0.5 1];

y = freqz(num, den, w)/M;

subplot(2,1,1)

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

title('Magnitude Spectrum |Y(e^{jomega})|')

xlabel('omega /pi');

ylabel('Amplitude');

subplot(2,1,2)

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

title('Phase Spectrum arg[Y(e^{jomega})]')

xlabel('omega /pi');

ylabel('Phase in radians');

滤波器所需长度=1

区别在于相位谱，我会选择前者来滤波，因为前者的相位谱相对平缓，后者在通带内有较大的相位突变。

精品文档

精品文档

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

式（4.36）：

clf;

num=[0.15 0 -0.15];

den=[1 -0.5 0.7];

L=input('Enter the number of samples L:');

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

stem(t,g);

title(['First',num2str(L),'samples of impulse response']);

xlabel('Time Index n');

ylabel('h[n]');

Enter the number of samples L:100

clf;

num=[0.15 0 -0.15];

den=[0.7 -0.5 1];

L=input('Enter the number of samples L:');

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

stem(t,g);

title(['First',num2str(L),'samples of impulse response']);

xlabel('Time Index n');

ylabel('h[n]');

Enter the number of samples L:100

精品文档

精品文档

两者反转。

Q4.6使用zplane分别生成式（4.36）和式（4.37）确定的两个滤波器的极零点图。讨论你的结果。

式（4.36）：num = [0.15 0 -1];

den = [1 -0.5 0.7];

zplane(num,den);

精品文档

精品文档

式（4.36）：num = [0.15 0 -1];

den = [0.7 -0.5 1];

zplane(num,den);

两者的零点一样，极点前者在圆内，后者在圆外。

Q4.19运行程序P4.3，生成每一类线性相位有限冲激响应滤波器的冲激响应。每一个有限冲精品文档

精品文档

激响应滤波器的长度是多少？验证冲激响应序列的对称性。接着验证这些滤波器的零点位置。使用MATLAB计算并绘制出这些滤波器的相位响应，验证它们的线性相位特性。这些滤波器的群延迟的多少？

（１）

长度分别是9,10,9,10；

（２）

精品文档

精品文档

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

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

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

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

（３）

pause

w = -4*pi:8*pi/511:4*pi;

h1=freqz(num1,1,w);

h2=freqz(num2,1,w);

h3=freqz(num3,1,w);

h4=freqz(num4,1,w);

subplot(2,2,1);

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

title('Phase Spectrum')

xlabel('omega /pi');

ylabel('Phase in radians'); grid;

title('Type 1 FIR Filter');

subplot(2,2,2);

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

title('Phase Spectrum')

xlabel('omega /pi');

ylabel('Phase in radians'); grid;

title('Type 2 FIR Filter');

subplot(2,2,3);

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

title('Phase Spectrum')

xlabel('omega /pi');

ylabel('Phase in radians'); grid;

title('Type 3 FIR Filter');

subplot(2,2,4);

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

精品文档

精品文档

title('Phase Spectrum')

xlabel('omega /pi');

ylabel('Phase in radians'); grid;

title('Type 4 FIR Filter');

（４）

pause

subplot(2,2,1);

grpdelay(h1);

title('Type 1 FIR Filter');

subplot(2,2,2);

grpdelay(h1);

title('Type 2 FIR Filter');

subplot(2,2,3);

grpdelay(h1);

title('Type 3 FIR Filter');

subplot(2,2,4);

grpdelay(h1);

title('Type 4 FIR Filter');

精品文档

精品文档

Q6.1使用程序P6.1，生成如下有限冲激响应传输函数的一个级联实现：

H1(z)210z123z234z331z416z54z6

画出级联实现的框图。H1(z)是一个线性相位传输函数吗？

num = input('Numerator coefficient vector = ');

den = input('Denominator coefficient vector = ');

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

sos = zp2sos(z,p,k)

Numerator coefficient vector = [2 10 23 34 31 16 4]

Denominator coefficient vector = [1 1 1 1 1 1 1]

sos =

2.0000 6.0000 4.0000 1.0000 -1.2470 1.0000

1.0000 1.0000 0.5000 1.0000 1.8019 1.0000

1.0000 1.0000 2.0000 1.0000 0.4450 1.0000

Q6.3使用程序P6.1，生成如下因果无限冲激响应传输函数的级联实现：

精品文档

精品文档

38z112z27z32z42z5

H1(z)1624z124z214z35z4z5画出级联实现的框图。

Numerator coefficient vector = [3 8 12 7 2 -2]

Denominator coefficient vector = [16 24 24 14 5 1]

sos =

0.1875 -0.0625 0 1.0000 0.5000 0

1.0000 2.0000 2.0000 1.0000 0.5000 0.2500

1.0000 1.0000 1.0000 1.0000 0.5000 0.5000

Q6.5使用程序P6.2生成式（６．２７）所示因果无限冲激响应传输函数的两种不同并联形式实现。画出两种实现的框图。

num = input('Numerator coefficient vector = ');

den = input('Denominator coefficient vector = ');

[r1,p1,k1] = residuez(num,den);

[r2,p2,k2] = residue(num,den);

disp('Parallel Form I')

disp('Residues are');disp(r1);

disp('Poles are at');disp(p1);

disp('Constant value');disp(k1);

disp('Parallel Form II')

disp('Residues are');disp(r2);

disp('Poles are at');disp(p2);

disp('Constant value');disp(k2);

Numerator coefficient vector = [3 8 12 7 2 -2]

Denominator coefficient vector = [16 24 24 14 5 1]

Parallel Form I

Residues are

-0.4219 + 0.6201i

-0.4219 - 0.6201i

2.3437

0.3437 - 2.5079i

0.3437 + 2.5079i

Poles are at

-0.2500 + 0.6614i

-0.2500 - 0.6614i

精品文档

精品文档

-0.5000

-0.2500 + 0.4330i

-0.2500 - 0.4330i

Constant value

-2

Parallel Form II

Residues are

-0.3047 - 0.4341i

-0.3047 + 0.4341i

-1.1719

1.0000 + 0.7758i

1.0000 - 0.7758i

Poles are at

-0.2500 + 0.6614i

-0.2500 - 0.6614i

-0.5000

-0.2500 + 0.4330i

-0.2500 - 0.4330i

Constant value

0.1875

Q8.1程序P8.1设计了什么类型的滤波器？其指标是什么？滤波器的阶数是多少？为了验证仿真，需计算多少个冲激响应样本？仿真是正确的吗？

s[0.1,0.8]，Rp1dB，滤波器类型为IIR带阻滤波器，其指标为p[0.4,0.5]，Rs30dB，滤波器阶数为4。需计算的冲激响应样本为２个。仿真正确。

Q8.3生成在习题Q8.1中产生的传输函数的一个级联实现，并编写出一个程序来仿真它，其精品文档

精品文档

中每个单独的部分用一个直接Ⅱ型实现。验证仿真。

format short

num=input('Numerator coefficients=');

den=input('Denominator coefficients=');

Numfactors=factorize(num);

Denfactors=factorize(den);

disp('Numerator Factor'),disp(Numfactors)

disp('Denominator Factor'),disp(Denfactors)

Numerator coefficients=[0.0571 0 -0.1143 0 0.0571]

Denominator coefficients=[1 -0.5099 1.2862 -0.335 0.4479]

Numerator Factors

Columns 1 through 3

1.0000 1.0211 0

1.0000 0.9793 0

1.0000 -1.0211 0

1.0000 -0.9793 0

Denominator Factors

Columns 1 through 3

1.0000 -0.5976 0.6785

1.0000 0.0877 0.6601

Q8.5生成在习题Q8.1中产生的传输函数的一个并联Ⅰ型实现，并编写出一个程序来仿真它，其中每个单独的部分用一个直接Ⅱ型实现。验证仿真。

num=input('Numerator coefficient vector=');

den=input('Denominator coefficient vector=');

[r1,p1,k1]=residuez(num,den);

[r2,p2,k2]=residue(num,den);

disp('Parallel Form 1')

disp('Residues are');disp(r1);

disp('Poles are at');disp(p1);

disp('Constant value');disp(k1);

disp('Parallel Form 2')

disp('Residues are');disp(r2);

disp('Poles are at');disp(p2);

disp('Constant value');disp(k2);

[b1,a1]=residuez(R1,P1,0);

R1=[r1(1) r1(2) r1(3) r1(4)];

P1=[p1(1) p1(2) p1(3) p1(4)];

disp('b1=');disp(b1);

disp('a1=');disp(a1);

Numerator coefficient vector=[0.0571 0 -0.1143 0 0.0571]

Denominator coefficient vector=[1 -0.5099 1.2862 -0.335 0.4479]

Parallel Form 1

Residues are

-0.0235 + 0.1977i

精品文档

精品文档

-0.0235 - 0.1977i

-0.0117 - 0.2130i

-0.0117 + 0.2130i

Poles are at

0.2988 + 0.7676i

0.2988 - 0.7676i

-0.0439 + 0.8113i

-0.0439 - 0.8113i

Constant value

0.1275

Parallel Form 2

Residues are

-0.1588 + 0.0411i

-0.1588 - 0.0411i

0.1734 - 0.0002i

0.1734 + 0.0002i

Poles are at

0.2988 + 0.7676i

0.2988 - 0.7676i

-0.0439 + 0.8113i

-0.0439 - 0.8113i

Constant value

0.0571

b1=

-0.0470 -0.2895

a1=

1.0000 -0.5976

Q8.9程序P8.３设计了什么类型的滤波器？其指标是什么？滤波器的阶数是多少？形成输入的正弦序列的频率是多少？

IIR低通椭圆滤波器，指标为p0.25，s0.55，Rp0.5dB，Rs50dB，阶数为8，输入的正弦序列的频率为f=0.7kHz。

Q8.10运行程序P8.3并产生两个图形。哪种输入成分会在滤波器输出出现？为什么输出序列的开始部分不是一种理想的正弦序列？修改程序P8.3，以便只过序列X２［ｎ］。产生的输出序列和预料一样吗？证明你的答案。

精品文档

精品文档

Input Sequence420-2-4Time index nOutput Sequence35404550AmplitudeAmplitude420-2-4Time index n35404550

答：产生的输出序列和预料是不一样的。

Q8.14程序P8.４设计了什么类型的滤波器？其指标是什么？滤波器的阶数是多少？为了验证仿真，需计算多少个冲激响应样本？仿真是正确的吗？

FIR低通滤波器，指标为p[0,0.3]，s[0.5,]，阶数为9，为了验证仿真需计算的冲激响应样本为10个。仿真正确。

实验名称：数字滤波器的设计

精品文档

精品文档

实验目的和任务：

（1）用窗口法设计满足指标的FIR数字滤波器。

（2）以模拟低通滤波器为原型设计IIR数字滤波器

（3）选定一个信号滤波问题，设计数字滤波器，验证滤波效果。

实验内容：

阅读Page 91-93相应的函数和程序P7.1,完成Q7.1, Q7.5,Q7.6

阅读Page 94-96相应的函数, 完成Q7.9,Q7.13,Q7.14,Q7.20

sinc函数的功能与使用，可通过help查询：

Matlab-Help-Search-Function Name-输入sinc

实验指导书Page 49有sinc函数使用实例

幅度响应的分析：

通过DTFT定义计算幅度响应

通过freqz函数分析

实验过程与结果分析：…参见实验一格式…

Q7.1用MATLAB确定一个数字无限冲击响应低通滤波器所有四种类型的最低阶数。指标如下：40kHz的抽样率，4kHz的通带边界频率，8kHz的阻带边界频率，0.5dB的通带波纹，40dB的最小阻带衰减。评论你的结果。

解：FT=40kHz，Fp=4kHz，Fs=8kHz，Rp=0.5dB，Rs=40dB。

p2FpFT241030.2

34010Wpp0.2

2Fs28103s0.4

3FT4010Wss0.4

（1）巴特沃兹滤波器：

[N,Wn]=buttord(0.2,0.4,0.5,40)

运行结果：

N =

8

Wn =

0.2469

答：最低阶数是8.

（2）切比雪夫1型滤波器：

[N,Wn]=cheb1ord(0.2,0.4,0.5,40)

运行结果：

N =

精品文档

精品文档

5

Wn =

0.2000

答：最低阶数是5.

（3）切比雪夫2型滤波器：

[N,Wn]=cheb2ord(0.2,0.4,0.5,40)

运行结果：

N =

5

Wn =

0.4000

答：最低阶数是5.

（4）椭圆滤波器：

[N,Wn]=ellipord(0.2,0.4,0.5,40)

运行结果：

N =

4

Wn =

0.2000

答：最低阶数是4.

Q7.5通过运行程序P7.1来设计巴特沃兹带阻滤波器。写出所产生的传输函数的准确表达式。滤波器的指标是什么？你的设计符合指标吗？使用MATLAB，计算并绘制滤波器的未畸变的相位响应及群延迟响应。

MATLAB程序：

Ws = [0.4 0.6]; Wp = [0.2 0.8]; Rp = 0.4; Rs = 50;

[N1, Wn1] = buttord(Wp, Ws, Rp, Rs);

[num,den] = butter(N1,Wn1,'stop');

disp('Numerator Coefficients are ');disp(num);

disp('Denominator Coefficients are ');disp(den);

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

plot(w/pi,g);grid

axis([0 1 -60 5]);

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

title('Gain Response of a Butterworth Bandstop Filter');

运行结果：

Numerator Coefficients are

0.0493 0.0000 0.2465 0.0000 0.4930 0.0000 0.4930 0.0000

0.2465 0.0000 0.0493

Denominator Coefficients are

1.0000 0.0000 -0.0850 0.0000 0.6360 0.0000 -0.0288 0.0000

0.0561 0.0000 -0.0008

0.04930.2465z20.493z40.493z60.2465z80.0493z10表达式：H(z)

10.085z20.636z40.0288z60.0561z80.0008z10精品文档

精品文档

指标：p10.2，p20.8，s10.4，s20.6，Rp0.4dB，Rs50dB。

编写的MATLAB程序如下（计算未畸变的相位响应及群延迟响应）：

% Program Q7_6

Ws = [0.4 0.6]; Wp = [0.2 0.8]; Rp = 0.4; Rs = 50;

% Estimate the Filter Order

[N1, Wn1] = buttord(Wp, Ws, Rp, Rs);

% Design the Filter

[num,den] = butter(N1,Wn1,'stop');

% Find the frequency response; find and plot unwrapped phase

wp = 0:pi/1023:pi;

wg = 0:pi/511:pi;

Hz = freqz(num,den,wp);

Phase = unwrap(angle(Hz));

figure(1);

plot(wp/pi,Phase);

grid;

% axis([0 1 a b]);

xlabel('omega /pi'); ylabel('Unwrapped Phase (rad)');

title('Unwrapped Phase Response of a Butterworth Bandstop Filter');

% Find and plot the group delay

GR = grpdelay(num,den,wg);

figure(2);

plot(wg/pi,GR);

grid;

%axis([0 1 a b]);

xlabel('omega /pi'); ylabel('Group Delay (sec)');

精品文档

精品文档

title('Group Delay of a Butterworth Bandstop Filter');

Q7.6修改程序P7.1来设计符合习题Q7.1所给指标的切比雪夫1型低通滤波器。写出所产生的传输函数的准确表达式。你的设计符合指标吗？使用MATLAB，计算并绘制滤波器的未畸变的相位响应及群延迟响应。

精品文档