南京邮电大学
实 验 报 告
实验名称：离散时间信号与系统的时、频域表示
离散傅立叶变换和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一个单位。
nnQ1.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^{j\\omega})') xlabel('\\omega /\\pi'); ylabel('Amplitude'); subplot(2,1,2)
plot(w/pi,imag(h));grid
title('Imaginary part of H(e^{j\\omega})') xlabel('\\omega /\\pi'); ylabel('Amplitude'); pause
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'); 修改程序后的运行结果为:
答：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^{j\\omega})') xlabel('\\omega /\\pi'); ylabel('Amplitude'); subplot(2,1,2)
plot(w/pi,imag(h));grid
title('Imaginary part of H(e^{j\\omega})') xlabel('\\omega /\\pi'); ylabel('Amplitude'); pause
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');
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
input('Hit
-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
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^{j\\omega})|') xlabel('\\omega /\\pi'); ylabel('Amplitude'); subplot(2,1,2)
plot(w/pi,angle(y));grid
title('Phase Spectrum arg[Y(e^{j\\omega})]') 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^{j\\omega})|') xlabel('\\omega /\\pi'); ylabel('Amplitude'); subplot(2,1,2)
plot(w/pi,angle(y));grid
title('Phase Spectrum arg[Y(e^{j\\omega})]') 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^{j\\omega})|') xlabel('\\omega /\\pi'); ylabel('Amplitude'); subplot(2,1,2)
plot(w/pi,angle(y));grid
title('Phase Spectrum arg[Y(e^{j\\omega})]') 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-4051015202530Time index nOutput Sequence35404550AmplitudeAmplitude420-2-4051015202530Time 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，计算并绘制滤波器的未畸变的相位响应及群延迟响应。
修改P7.1程序：
[N, Wn] = cheb1ord(0.2,0.4,0.5,40);Rp=0.5;
[num,den] = cheby1(N,Rp,Wn);[g, w] = Gain(num,den); % Plot the gain response figure(1);
plot(w/pi,g);grid; axis([0 1 -60 5]);
xlabel('\\omega /\\pi'); ylabel('Gain in dB');
title('Gain Response of a Type 1 Chebyshev Lowpass Filter'); figure(2);
w2 = 0:pi/511:pi;
Hz = freqz(num,den,w2); Phase = unwrap(angle(Hz)); plot(w2/pi,Phase);grid;
xlabel('\\omega /\\pi'); ylabel('Unwrapped Phase (rad)'); title('Unwrapped Phase Response of a Type 1 Chebyshev Lowpass Filter'); % Find and plot the group delay figure(3);
GR = grpdelay(num,den,w2); plot(w2/pi,GR);grid;
xlabel('\\omega /\\pi'); ylabel('Group Delay (sec)');
title('Group Delay of a Type 1 Chebyshev Lowpass Filter');
Gain Response of a Type 1 Chebyshev Lowpass Filter0-10-20Gain in dB-30-40-50-6000.10.20.30.40.5? /?0.60.70.80.91