数字信号处理实验指导书（学生版）

数字信号处理实验指导书（学⽣版）
“数字信号处理”实验指导书（⼀）
⼀、实验课程编码：105003 ⼆、实验课程名称：数字信号处理
三、实验项⽬名称：应⽤MATLAB 分析离散信号频谱四、实验⽬的
掌握应⽤MATLAB 分析离散信号频谱的⽅法，即熟悉应⽤MATLAB 分析离散信号的函数。

五、主要设备
安装有MATLAB 软件的电脑六、实验内容
编写MATLAB 程序，实现下⾯题⽬：
1. ⽤快速卷积法计算下⾯两个序列的线性卷积。

)()4.0(s )(15n R n in n x =，)(9.0)(20n R n h n =
2．已知序列[]()
cos 0120
n n N N
x n π?≤≤-?=?
其它
（1）计算该序列DTFT 的表达式()j X e ω，并画出N=10时的()j X e ω曲线；（2）编写MATLAB 程序，利⽤FFT 函数，计算N ＝10时，序列x [k ]的DTFT 在2m m
N
πω＝的抽样值。

利⽤hold 函数，将抽样点画在()j X e ω的曲线上。

3.理解⾼密度频谱和⾼分辨率频谱的概念。

设)52.0cos()48.0(co )(n n s n x ππ+=
（1）取0≤n ≤9，求)(1k X
（2）将（1）中的)(x n 补零加长到0≤n ≤99，求)(2k X （3）增加取样值的个数，取0≤n ≤99，求)(3k X
4. ⽤DFT 对连续信号做谱分析。

设)50cos()100sin()200cos()(t t t t x a πππ++=，⽤DFT 分析)(t x a 的频谱结构，选择不同的截取长度Tp ，观察截断效应，试⽤加窗的⽅法减少谱间⼲扰。

选取的参数：
（1）频率s s f T Hz f /1 ,400==
（2）采样信号序列)()()(n w nT x n x a =，)(n w 是窗函数。

选取两种窗函数：
矩形窗函数)()(n R n w N =和Hamming 窗，后者在程序中调⽤函数Hamming 产⽣宽度为N 的Hamming 窗函数向量。

（3）对)(n x 做2048点DFT ，作为)(t x a 的近似连续频谱)(jf X a 。

其中N 为
采样点数，p s T f N =，p T 为截取时间长度，取三种长度0.04s 、2×0.04s 、4×0.04s 、8×0.04s 。

5. 已知⼀连续信号为()()x t t ε-3t ＝e ，试利⽤DFT 近似分析其频谱。

要求频率分辨率为1Hz ，确定抽样频率s f 、抽样点数N 以及持续时间T p 。

说明：连续信号x (t )的频谱X (j Ω)可以由其离散信号x (n )的DFT 近似求得： X (j Ω) ≈ T ．FFT[x (n )]，其中T 为采样周期。

计算连续信号x (t )的频谱，要考虑⼏个问题：1）频率混叠：如果x (t )不是有限带宽的，则采样频率f s 选择的依据是，使由时域采样所造成的频率混叠⼩到可以忽略；2）频率分辨率：频率分辨率和信号在时域的持续时间成反⽐。

3）截断效应：如果x (t )是⽆限长的，就必须进⾏截断，截断会引起吉布斯效应（波动），也会把窗函数的频谱引⼊信号频谱，造成混叠。

分析连续信号频谱的⼀般步骤为：1）先初步选择时间记录长度T 0,使得其包括了⼤部分⾮零的数据，然后⽤逐次增⼤采样频率fs 的⽅法来选择采样频率。

最终使得时域采样造成的频率混叠可以近似的忽略不计。

2）在确定了采样频率后，进⼀步选择时间
记录长度T 0，即逐次将T 0加倍，因为此时采样频率相同，两个频谱之间的差别是由截断效应不同引起的，如果两个频谱的误差不⼤，则这个记录长度T 0可以接受。

（本题直接给出了频率分辨率的要求，也就是记录长度T 0可以直接确定（不必进⾏步骤2），只需进⾏步骤1逐次选择采样频率。

）
6．验证频域采样定理。

（1）产⽣⼀个三⾓波序列x(n)，长度为M=40；
（2）计算N=64时的X(k)=DFT[x(n)]，并图⽰x(n)和X(k)
2,0[ 上进⾏32点抽样，得到X1(k)＝X(2k)，k＝0,1, (31)
（3）对X(k)在]
（4）求X1(k)的32点IDFT，即x1(n)＝IDFT[X1(k)]。

（5）绘制x1((n))32的波形图，观察x1((n))32和x(n)的关系，并加以说明。

七、实验步骤
1、熟悉与离散信号频谱分析相关的MATLAB函数（参考附录1）；
2、通过运⾏附录2中提供的例题，熟悉⽤MATLAB分析离散信号频谱的基本⽅法；
3、根据“六、实验内容”中各个题⽬的要求，编写MATLAB程序代码，调试程序，分析并保存结果。

⼋、实验结果
对实验练习题编写MATLAB程序并运⾏，在计算机上输出仿真结果。

附录1 主要的相关MATLAB函数
1.fft.m和ifft.m
调⽤格式：〔X〕＝fft（x）
〔x〕＝ifft（X）
〔X〕＝fft（x,N）
〔x〕＝ifft（X,N）
2.czt.m
调⽤格式：〔y〕＝czt（x,m,w,s）
3.fftshift.m
调⽤格式：〔y〕＝fftshift（x）
附录2 例题
例1利⽤DFT的性质，编写MATLAB程序，计算下列序列的6点圆周卷积。

（1）x [n]= {1,-3,4,2,0,-2}, h[n]= {3,0,1,-1,2,1}
（2）x[n]=cos(πn/2), h[n]=3n, n=0,1,2,3,4,5 [MATLAB 程序]：
N=6;
xn=[1,-3,4,2,0,-2]; hn=[3,0,1,-1,2,1];
Xk=fft(xn,N); %计算N 点的DFT[x(n)] Hk=fft(hn,N); %计算N 点的DFT[h(n)] Yk=Xk.*Hk; %DFT[x(n)].*DFT[h(n)]
y=ifft(Yk,N) %计算N 点的IDFT[Y(k)],即为x(n)和h(n)的圆周卷积 [运⾏结果]： y =
6.0000 -3.0000 1
7.0000 -2.0000 7.0000 -13.0000 [MATLAB 程序]：
N=6; n=0:N-1; xn=cos(pi*n/2); hn=3*n;
Xk=fft(xn,N); %计算N 点的DFT[x(n)] Hk=fft(hn,N); %计算N 点的DFT[h(n)] Yk=Xk.*Hk; %DFT[x(n)].*DFT[h(n)] y=ifft(Yk,N) %计算N 点的IDFT[Y(k)],即为x(n)和h(n)的圆周卷积 [运⾏结果]： y =
-6.0000 -3.0000 18.0000 21.0000 6.0000 9.0000
例2、基本序列的离散傅⽴叶变换计算复正弦序列：)(e
)(8
1n R n x N n j π
=，余弦序列：)()8
cos()(2n R n n x N π
=
分别对以上序列求当N=16和N=8时的DFT，并绘出幅频特性曲线，对其结果进⾏分析。

[MATLAB程序]：
%基本序列的离散傅⽴叶变换计算
N=16;N1=8;
n=0:N-1;
k=0:N1-1;
x1n=exp(j*pi*n/8); %产⽣x1(n)
X1k=fft(x1n,N); %计算N点的DFT[x1(n)]
X2k=fft(x1n,N1); %计算N1点的DFT[x1(n)]
x2n=cos(pi*n/8); %产⽣x2(n)
X3k=fft(x2n,N); %计算N点的DFT[x2(n)]
X4k=fft(x2n,N1); %计算N1点的DFT[x2(n)]
subplot(2,2,1);
stem(n,abs(X1k),'.');
axis([0,20,0,20]);
ylabel('|X1(k)|');
title('16点的DFT[x1(n)]');
subplot(2,2,3);
stem(k,abs(X2k),'.');
axis([0,20,0,20]);
ylabel('|X1(k)|');
title('8点的DFT[x1(n)]');
subplot(2,2,2);
stem(n,abs(X3k),'.');
axis([0,20,0,20]);
ylabel('|X2(k)|');
title('16点的DFT[x2(n)]');
subplot(2,2,4);
stem(k,abs(X4k),'.');
axis([0,20,0,20]);
ylabel('|X2(k)|');
title('8点的DFT[x2(n)]');
[运⾏结果]：
例3、验证N 点DFT 的物理意义
已知)()(4n R n x =，ω
ω
ω
j j j e e n x e X ----=
=11)]([DFT )(4，绘制相应的幅频和相频曲线，
并计算N=8和N=16时的DFT 。

[MATLAB 程序]：
%验证N 点DFT 的物理意义
N1=8; N2=16; n=0:N1-1; k1=0:N1-1; k2=0:N2-1;
w=2*pi*(0:2047)/2048;
Xw=(1-exp(-j*4*w))./(1-exp(-j*w));%对x(n)的频谱函数采样2048个点可以近似看作是连续的频谱
xn=[(n>=0&n<4)]; % 产⽣x(n)，即将n 中第0到第3个数变为1，其余为零 X1k=fft(xn,N1); %计算N1点的DFT X2k=fft(xn,N2); %计算N2点的DFT
subplot(3,2,1);plot(w/pi,abs(Xw));xlabel('w/pi');
subplot(3,2,2);plot(w/pi,angle(Xw));axis([0,2,-pi,pi]);line([0,2],[0,0]);xlabel('w/pi');
subplot(3,2,3);stem(k1,abs(X1k),'.');xlabel('k(w=2pik/N1)');ylabel('|X1(k)|');hold on; plot(N1/2*w/pi,abs(Xw)); %在频谱上叠加连续频谱的幅度曲线
subplot(3,2,4);stem(k1,angle(X1k),'.');axis([0,N1,-pi,pi]);line([0,N1],[0,0]); xlabel('k(w=2pik/N1)');ylabel('Arg[X1(k)]');hold on; plot(N1/2*w/pi,angle(Xw)); %在频谱上叠加连续频谱的相位曲线
subplot(3,2,5);stem(k2,abs(X2k),'.');axis([0,N2,0,4]); xlabel('k(w=2pik/N2)');ylabel('|X2(k)|');hold on; plot(N2/2*w/pi,abs(Xw)); subplot(3,2,6);stem(k2,angle(X2k),'.');axis([0,N2,-pi,pi]);line([0,N2],[0,0]); xlabel('k(w=2pik/N2)');ylabel('Arg[X2(k)]');hold on; plot(N2/2*w/pi,angle(Xw));
[运⾏结果]：
“数字信号处理”实验指导书（⼆）
⼀、实验课程编码：103044 ⼆、实验课程名称：数字信号处理A
三、实验项⽬名称：应⽤MATLAB 设计IIR 数字滤波器四、实验⽬的
掌握应⽤MATLAB 设计IIR 数字滤波器的⽅法，即熟悉应⽤MATLAB 设计IIR 滤波器相关的函数。

五、主要设备
安装有MATLAB 软件的电脑六、实验内容
编写MATLAB 程序，实现下列题⽬：
1．⽤冲激响应不变法设计巴特沃思数字低通滤波器，采样频率为10kHz ，通带截⾄频率1.5kHz ，通带最⼤衰减3dB ，阻带截⾄频率3kHz ，阻带最⼩衰减为12dB 。

画出所设计的滤波器的幅度响应。

2．⽤双线性法设计巴特沃思低通数字滤波器，采样频率10kHz ，通带截⾄频率2.5kHz ，通带最⼤衰减2dB ，阻带截⾄频率3.5kHz ，阻带最⼩衰减15dB 。

画出所设计的滤波器的幅度响应。

3．使⽤双线性变换法设计低通数字滤波器，要求⽤切⽐雪夫I 型滤波器逼近，设计指标为：
ωp =0.4π , R p =1dB , ωs =0.54π, R S =15dB 画出所设计的滤波器的幅度响应。

4．利⽤双线性变换法设计数字⾼通滤波器，要求⽤切⽐雪夫I 型滤波器逼近，设计指标为：
dB
15 ,46.0dB 1 ,6.0====s s p p A R πωπω
画出所设计的滤波器的幅度响应。

5．使⽤双线性变换法设计带通数字滤波器，要求⽤切⽐雪夫I 型滤波器逼近，设计指标为：
dB 15 ,7.0 ,2.0dB 1 ,5.0 ,4.02121======s s s p p p A R πωπωπωπω
画出所设计的滤波器的幅度响应。

七、实验步骤
1、熟悉与设计滤波器相关的MATLAB 函数（参考附录1）；
2、通过运⾏附录2中提供的例题，熟悉⽤MATLAB 设计IIR 数字滤波器的基本⽅法；
3、根据“六、实验内容”中各个题⽬的要求，编写MATLAB 程序代码，调试程序，分析并保存结果。

⼋、实验结果对实验练习题编写MATLAB 程序并运⾏，在计算机上输出仿真结果。

附录1 主要的相关MATLAB 函数 1．buttord.m 调⽤格式：[N,Wn]＝buttord ( Wp ，Ws ，Rp ，Rs ) [N,Wn]＝buttord ( Wp ，Ws ，Rp ，Rs ，’s ’) 2．buttap.m
调⽤格式：[z,p,k] ＝buttap （N ） 3．butter.m
调⽤格式：[B,A] = butter(N,Wn)
[B,A] = butter(N,Wn,'s')
4．cheb1ord.m
调⽤格式：[n,Wn] = cheb1ord(Wp, Ws, Rp, Rs)；
[n,Wn] = cheb1ord(Wp, Ws, Rp, Rs, 's')
5．cheb1ap.m 6．cheby1.m
调⽤格式：[B,A] = cheby1(N,Rp,Wn)
[B,A] = cheby1(N,Rp,Wn,’s ’)
7．cheb2ord.m
8．cheb2ap.m
9．cheby2.m
10．ellipord.m
11．ellip.m
⽂件4～11的调⽤格式与⽂件1~3的调⽤格式类似。

12．bilinear.m
调⽤格式：[Bz,Az]＝bilinear（B，A，Fs）13．impinvar.m
调⽤格式：[Bz,Az]＝impinvar（B，A，Fs）14．Grpdelay.m
调⽤格式：[Gd,W] = grpdelay (B,A, N)
15．lp2hpz.m （⾃编函数）
调⽤格式：[bz,az] = lp2hpz(b,a,alpha)
16．lp2bpz.m （⾃编函数）
调⽤格式：[bz,az] = lp2bpz(b,a,alpha1,alpha2)
17．lp2bsz.m （⾃编函数）
调⽤格式：[bz,az] = lp2bsz(b,a,alpha1,alpha2) 18．freqzn.m（⾃编函数）
调⽤格式：[Rp,As] = freqzn(num,den,wp,ws,Rp,As,ftype) 数字域通带转换公式：
附录2 例题
例1．设计⼀个低通模拟巴特沃斯滤波器，使其满⾜：通带：rad/s) 1020 ( dB 141?≤≤≤πΩδ阻带：rad/s) 105.12 ( dB 1542??≥≥πΩδ
[MATLAB 程序]： omegaP = 20000*pi; omegaS = 30000*pi; Rp = 1; As = 15;
[N, Wn] = buttord(omegaP, omegaS ,Rp,As,'s') %%得到模拟滤波器的阶次和3dB 截⽌频率
[B,A] = butter (N,Wn,'s') %%由阶次和3dB 截⽌频率得到模拟滤波器的系统函数
%% 画出模拟滤波器的幅频特性 [H,w] = freqs (B,A); figure;
subplot(1,2,1), plot(w/pi,abs(H)); %%画出幅频特性（绝对幅度）
title('Magnitude Response');
xlabel('Analog frequency in pi units'); ylabel('abs(H)');
axis([0,50000,0,1.1]); %% 只画出0~50000部分的图形
line([20000,20000],[0,1.1],'LineStyle',':');
line([30000,30000],[0,1.1],'LineStyle',':');
rip = 10^(-Rp/20); %% 由通带波动得到通带容限
att = 10^(-As/20); %% 由阻带衰减得到阻带容限
line([0,50000],[rip,rip],'LineStyle',':');
line([0,50000],[att,att],'LineStyle',':');
set(gca, 'XTick', [0,20000,30000,50000]);
set(gca, 'YTick', [0,att,rip,1]);
subplot(1,2,2),plot(w/pi,20*log10(abs(H))); %%画出幅频特性（dB）title('Magnitude in dB');
xlabel('Analog frequency in pi units'); ylabel('dB');
axis([0,50000,-20,2]);
line([20000,20000],[-20,2],'LineStyle',':');
line([30000,30000],[-20,2],'LineStyle',':');
line([0,50000],[-As,-As],'LineStyle',':');
line([0,50000],[-Rp,-Rp],'LineStyle',':');
set(gca, 'XTick',[0,20000,30000,50000]);
set(gca, 'YTick',[-20,-As,-Rp,0]);
[运⾏结果]：
N = 6
Wn = 7.0865e+004
B = 1.0e+029
0 0 0 0 0 0 1.2665
A = 1.0e+029 *
0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 1.2665
例2．设计⼀个低通模拟切⽐雪夫I 型滤波器，使其满⾜：通带截⾄频率：rad/s 1024?=πΩc ，通带波纹：dB 11=δ阻带起始频率：rad/s 105.124??=πΩs ，阻带衰减：dB 152=δ
[MATLAB 程序]： omegaP = 2*pi*10000; omegaS = 3*pi*10000; Rp=1; As=15;
[N,Wn] = cheb1ord (omegaP, omegaS, Rp, As, 's') %%得到模拟滤波器的阶次 [b,a] = cheby1(N,Rp,Wn,'s') %%得到模拟滤波器的系统函数 %%%%%% 画出幅频特性 [H,w] = freqs (b,a); figure;
subplot(1,2,1),plot(w/pi,abs(H)); title('Magnitude Response'); xlabel('Analog frequency in pi units'); ylabel('abs(H)');
axis([0,50000,0,1.1]);
line([20000,20000],[0,1.1],'LineStyle',':');
line([30000,30000],[0,1.1],'LineStyle',':');
rip = 10^(-Rp/20);
att = 10^(-As/20);
line([0,50000],[rip,rip],'LineStyle',':');
line([0,50000],[att,att],'LineStyle',':');
set(gca,'XTickMode','manual','XTick',[0,20000,30000,50000]); set(gca,'YTickMode','manual','YTick',[0,att,rip,1]);
subplot(1,2,2),plot(w/pi,20*log10(abs(H)));
title('Magnitude in dB');
xlabel('Analog frequency in pi units');
ylabel('dB');
axis([0,50000,-20,2]);
line([20000,20000],[-20,2],'LineStyle',':');
line([30000,30000],[-20,2],'LineStyle',':');
line([0,50000],[-As,-As],'LineStyle',':');
line([0,50000],[-Rp,-Rp],'LineStyle',':');
set(gca,'XTickMode','manual','XTick',[0,20000,30000,50000]); set(gca,'YTickMode','manual','YTick',[-20,-As,-Rp,0]);
[运⾏结果]：
N = 4
Wn = 6.2832e+004
b = 1.0e+018*
0 0 0 0 3.8286
a = 1.0e+018*
0.0000 0.0000 0.0000 0.0002 4.2958
例3．利⽤脉冲响应不变法设计⼀个巴特沃斯数字低通滤波器，使其满⾜：
dB
15 ,3.0dB 1 ,2.021====δπωδπωs p
[MATLAB 程序]： %% 数字滤波器的设计指标 wp = 0.2*pi; ws = 0.3*pi; Rp = 1; As = 15;
%% 将数字滤波器设计指标转换为模拟滤波器的设计指标 %% 假设采样周期T=1; T = 1;
omegaP = wp/T; omegaS = ws/T;
%% 得到模拟滤波器的阶次
[N,Wn] = buttord (omegaP, omegaS, Rp, As, 's') %% 得到模拟滤波器的系统函数 [b,a] = butter (N,Wn,'s')
%% 使⽤冲激响应不变法将模拟滤波器转换为数字滤波器 [bz,az] = impinvar(b,a,1/T);
%% 画出所设计数字滤波器的幅频特性，并检测Rp ，As 是否满⾜设计指标 [Rpp,Ass] = freqzn(bz, az, wp/pi, ws/pi, Rp, As, 'low') %%（⾃编函数）
[运⾏结果]： N = 6 Wn = 0.7087
b = 0 0 0 0 0 0 0.1266 a = 1.0000 2.7380 3.7484 3.2533 1.8824 0.6905 0.1266 bz = 0.0000 0.0007 0.0105 0.0167 0.0042
0.0001 0 az = 1.0000 -3.3443 5.0183 -4.2190 2.0725 -0.5600 0.0647 Rpp = 0.9202 Ass = 15
例4．利⽤双线性变换法设计⼀个切⽐雪夫I 型数字低通滤波器，使其满⾜：
dB
15 ,3.0dB 1 ,2.021====δπωδπωs p
[MATLAB程序]：
%% 数字滤波器的设计指标
wp = 0.2*pi;
ws = 0.3*pi;
Rp = 1;
As = 15;
%% 将数字滤波器设计指标转换为模拟滤波器的设计指标
%% 当使⽤双线性变化法时要预畸变!
%% 假设采样周期T=1;
T = 1;
omegaP = (2/T)*tan(wp/2);
omegaS = (2/T)*tan(ws/2);
%% 得到模拟滤波器的阶次
[N,Wn] = cheb1ord (omegaP,omegaS,Rp,As,'s')
%% 设计模拟滤波器
[b,a] = cheby1(N,Rp,Wn,'s')
%% 使⽤双线性变换法将模拟滤波器转换为数字滤波器
[bz,az] = bilinear(b,a,1/T)
%% 画出所设计数字滤波器的幅频特性，并检测Rp，As是否满⾜设计指标[Rpp,Ass] = freqzn(bz, az, wp/pi, ws/pi, Rp, As, 'low') %%（⾃编函数）%%%%%% 直接在数字域设计
[N,Wn] = cheb1ord (wp/pi, ws/pi,Rp,As)
[b,a] = cheby1(N,Rp,Wn,'low')
[运⾏结果]：
N = 4
Wn = 0.6498
b = 0 0 0 0 0.0438
a = 1.0000 0.6192 0.6140 0.2038 0.0492
bz = 0.0018 0.0073 0.0110 0.0073 0.0018
az = 1.0000 -3.0543 3.8290 -2.2925 0.5507 Rpp = 0.9997 Ass = 24 N = 4 Wn = 0.2000
b = 0.0018 0.0073 0.0110 0.0073 0.0018 a = 1.0000 -3.0543 3.8290 -2.2925
0.5507
例5．利⽤切⽐雪夫I 型，使⽤双线性变换设计数字⾼通滤波器，使其指标满⾜：
dB ,.dB ,.20520165021=δπ=θ=δπ=θs p
[MATLAB 程序]： %% 数字⾼通滤波器的设计指标 wp = 0.65*pi; ws = 0.52*pi; Rp = 1; As = 20;
%% 1）将数字⾼通设计指标转化为数字低通设计指标
%% 设数字低通滤波器的通带截⾄频率为cetaP = 0.3*pi;
cetaP = 0.3*pi;
alpha = -cos((wp+cetaP)/2)/cos((wp-cetaP)/2)
cetaS = log(-(exp(j*ws)+alpha)/(1+alpha*exp(j*ws)))/(-j)
%% 2）将数字低通滤波器设计指标转换为模拟低通滤波器的设计指标
%% 使⽤双线性转换，故需要预畸变!
%% 假设采样周期T=1;
T = 1;
omegaP = (2/T)*tan(cetaP/2);
omegaS = (2/T)*tan(cetaS/2);
%% 3）设计模拟低通滤波器
%% 得到模拟滤波器的阶次
[N,Wn] = cheb1ord (omegaP,omegaS,Rp,As,'s')
%% 得到模拟滤波器的系统函数
[b,a] = cheby1(N,Rp,Wn,'s')
%% 4）使⽤双线性变化法将模拟滤波器转换为数字滤波器
[Bz,Az] = bilinear(b,a,1/T)
%% 5）将数字低通滤波器转化为数字⾼通滤波器
[bz,az] = lp2hpz(Bz,Az,alpha) %%（该函数为⾃编函数见附录）
%% 画出所设计数字滤波器的幅频特性，并检测Rp，As是否满⾜设计指标[Rpp,Ass] = freqzn(bz, az, wp/pi, ws/pi, Rp, As, 'high') %%（⾃编函数）%%%%%% 直接在数字域设计
[N,Wn] = cheb1ord (wp/pi,ws/pi,Rp,As) [b,a] = cheby1(N,Rp,Wn,'high')
[运⾏结果]：
alpha = -0.0920
cetaS = 1.3258
N = 4。