数字信号处理实验指导书(学生版)

一、实验课程编码：105003 二、实验课程名称：数字信号处理

三、实验项目名称： 应用MATLAB分析离散信号频谱 四、实验目的

掌握应用MATLAB分析离散信号频谱的方法，即熟悉应用MATLAB分析离散信号的函数。 五、主要设备

安装有MATLAB软件的电脑 六、实验内容

编写MATLAB程序，实现下面题目：

1. 用快速卷积法计算下面两个序列的线性卷积。

x(n)sin(0.4n)R15(n)，h(n)0.9nR20(n)

cosn2N2．已知序列xn00nN1其它

（1）计算该序列DTFT的表达式Xej，并画出N=10时的Xej曲线； （2）编写MATLAB程序，利用FFT函数，计算N＝10时，序列x[k]的DTFT在m＝2m

N的抽样值。利用hold函数，将抽样点画在Xej的曲线上。

3.理解高密度频谱和高分辨率频谱的概念。 设x(n)cos(0.48n)cos(0.52n)

（1） 取0≤n≤9，求X1(k)

（2） 将（1）中的x(n)补零加长到0≤n≤99，求X2(k) （3） 增加取样值的个数，取0≤n≤99，求X3(k)

4. 用DFT对连续信号做谱分析。

设xa(t)cos(200t)sin(100t)cos(50t)，用DFT分析xa(t)的频谱结构，选择不同的截取长度Tp，观察截断效应，试用加窗的方法减少谱间干扰。

选取的参数：

（1） 频率fs400Hz, T1/fs

（2） 采样信号序列x(n)xa(nT)w(n)，w(n)是窗函数。选取两种窗函数：

矩形窗函数w(n)RN(n)和Hamming窗，后者在程序中调用函数Hamming产生宽度为N的Hamming窗函数向量。

（3） 对x(n)做2048点DFT，作为xa(t)的近似连续频谱Xa(jf)。其中N为

采样点数，NfsTp，Tp为截取时间长度，取三种长度0.04s、2×0.04s、4×0.04s、8×0.04s。

5. 已知一连续信号为xt＝e-3tt，试利用DFT近似分析其频谱。要求频率分辨率为1Hz，确定抽样频率fs、抽样点数N以及持续时间Tp。

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

计算连续信号x(t)的频谱，要考虑几个问题：1）频率混叠：如果x(t)不是有限带宽的，则采样频率fs选择的依据是，使由时域采样所造成的频率混叠小到可以忽略；2）频率分辨率：频率分辨率和信号在时域的持续时间成反比。3）截断效应：如果x(t)是无限长的，就必须进行截断，截断会引起吉布斯效应（波动），也会把窗函数的频谱引入信号频谱，造成混叠。

分析连续信号频谱的一般步骤为：1）先初步选择时间记录长度T0,使得其包括了大部分非零的数据，然后用逐次增大采样频率fs的方法来选择采样频率。最终使得时域采样造成的频率混叠可以近似的忽略不计。2）在确定了采样频率后，进一步选择时间记录长度T0，即逐次将T0加倍，因为此时采样频率相同，两个频谱之间的差别是由截断效应不同引起的，如果两个频谱的误差不大，则这个记录长度T0可以接受。

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

6．验证频域采样定理。

（1） 产生一个三角波序列x(n)，长度为M=40；

（2） 计算N=64时的X(k)=DFT[x(n)]，并图示x(n)和X(k)

（3） 对X(k)在[0,2]上进行32点抽样，得到X1(k)＝X(2k)，k＝0,1,…,31。 （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 17.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、基本序列的离散傅立叶变换计算 复正弦序列：x1(n)ejn8RN(n)，余弦序列：x2(n)cos(n)RN(n)

8分别对以上序列求当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的物理意义

1ej4已知x(n)R4(n)，X(e)DFT[x(n)]1ejj，绘制相应的幅频和相频曲线，

并计算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π , Rp =1dB, ωs=0.54π, RS=15dB 画出所设计的滤波器的幅度响应。

4．利用双线性变换法设计数字高通滤波器，要求用切比雪夫I型滤波器逼近，设计指标为：

p0.6, Rp1dBs0.46, As15dB

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

5．使用双线性变换法设计带通数字滤波器，要求用切比雪夫I型滤波器逼近，设计指标为：

p10.4, p20.5, Rp1dBs10.2, s20.7, As15dB画出所设计的滤波器的幅度响应。 七、实验步骤

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.

[z,p,k] ＝cheb1ap(N,Rp)模拟滤波器 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) 数字域通带转换公式： 变换类型 低通-高通 低通-带通 z1变换公式 变换参数的公式 Z1 11Zcos((cc)/2) cos((cc)/2) z1Z21Z12 212Z1Z1cos((21)/2)，cos((21)/2)1kcot(2)tanc 2212kk1，2 k1k1ω2，ω1为带通滤波器通带的上、下截止频率 低通-带阻 z1Z21Z12 2Z21Z11cos((21)/2)，cos((21)/2)1ktan(2)tanc 22121k，2 1k1kω2，ω1为带阻滤波器通带的上、下截止频率

附录2 例题

例1．设计一个低通模拟巴特沃斯滤波器，使其满足： 通带：11 dB ( 02104 rad/s) 阻带：215 dB ( 21.5104 rad/s)

[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 A = 1.0e+029 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 1.2665

1.2665

例2．设计一个低通模拟切比雪夫I型滤波器，使其满足： 通带截至频率：c2104 rad/s，通带波纹：11 dB 阻带起始频率：s21.5104 rad/s，阻带衰减：215 dB

[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．利用脉冲响应不变法设计一个巴特沃斯数字低通滤波器，使其满足：

p0.2, 11dBs0.3, 215dB

[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型数字低通滤波器，使其满足：

p0.2, 11dBs0.3, 215dB

[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型，使用双线性变换设计数字高通滤波器，使其指标满足：

p0.65, 11dBs0.52, 220dB[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

Wn = 1.0191

b = 0 0 0 0 0.2649 a = 1.0000 0.9710 1.5098 0.7859 0.2972 Bz = 0.0084 0.0335 0.0502 0.0335 0.0084 Az = 1.0000 -2.3741 2.7057 -1.5917 0.4103 bz = 0.0148 -0.0592 0.0888 -0.0592 0.0148 az = 1.0000 1.9963 2.1975 1.2902 0.3548 Rpp = 1.0000 Ass = 23 N = 4 Wn = 0.6500

b = 0.0148 -0.0592 0.0888 -0.0592 0.0148 a = 1.0000 1.9963 2.1975 1.2902 0.3548

例6．利用切比雪夫I型设计一个数字带阻滤波器，使其满足：

p10.25, p20.75, Rp1 dBs10.35, s20.65, As20 dB

[MATLAB程序]： %% 数字带阻滤波器的设计指标 wp1 = 0.25*pi; wp2 = 0.75*pi; ws1 = 0.35*pi; ws2 = 0.65*pi; Rp = 1; As = 20;

%% 1）将数字带阻设计指标转化为数字低通设计指标 %% 设数字低通滤波器的通带截至频率为 cetaP = 0.2*pi; cetaP = 0.2*pi;

beta = cos((wp2+wp1)/2)/cos((wp2-wp1)/2) k = tan((wp2-wp1)/2)*tan(cetaP/2) alpha1 = (-2*beta)/(1+k) alpha2 = (1-k)/(1+k)

cetaS=log((exp(-j*2*ws1)+alpha1*exp(-j*ws1)+alpha2)/(alpha2*exp(-j*2*ws1)+alpha1*exp(-j*ws1)+1))/(-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] = lp2bsz(Bz,Az,alpha1,alpha2) %%（自编函数）

%% 画出所设计数字滤波器的幅频特性，并检测Rp，As是否满足设计指标 [Rpp,Ass] = freqzn(bz, az, [wp1,wp2]/pi, [ws1,ws2]/pi, Rp, As, 'stop') %%（自编函数）

%%%%%% 直接在数字域设计

[N,Wn] = cheb1ord ([wp1,wp2]/pi,[ws1,ws2]/pi,Rp,As) [b,a] = cheby1(N,Rp,Wn,'stop') [运行结果]： beta = 8.6596e-017 k = 0.3249 alpha1 = -1.3072e-016 alpha2 = 0.5095 cetaS = 1.1353 - 0.0000i N = 3 Wn = 0.6498

b = 0 0 0 a = 1.0000 0.6423 0.5230 Bz = 0.0115 0.0344 0.0344 Az = 1.0000 -2.1378 1.7693 bz = 0.1321 -0.0000 0.3964 az = 1.0000 -0.0000 -0.3432 Rpp = 1.0000 Ass = 22 N = 3

Wn = 0.2500 0.7500

0.1348 0.1348 0.0115 -0.5398

-0.0000 -0.0000 0.3964 0.6044 -0.0000 0.0000 0.1321 -0.2041 b = 0.1321 -0.0000 0.3964 -0.0000 0.3964 -0.0000 0.1321 a = 1.0000 0 -0.3432 0.0000 0.6044 -0.0000 -0.2041

例题中用到的自编函数为：

%% Transform lowpass digital filters to highpass %% 2011-11-13 by nlp

function [bz,az] = lp2hpz(b,a,alpha) N = length(a); %% N >=2 bz = zeros(1,N); az = zeros(1,N); bb = [-alpha,-1]; aa = [1,alpha];

for i = 1:N x1 = 1;

for j = 1:(N-i)

x1 = conv(x1,aa); end x2 = 1;

for j = 1 : (i-1)

x2 = conv(x2,bb);

end

az = az + a(i)*conv(x1,x2); bz = bz + b(i)*conv(x1,x2); end

bz = bz/az(1); az = az/az(1);

%% Transform lowpass digital filters to bandpass %% 2011-11-13 by nlp

function [bz,az] = lp2bpz(b,a,alpha1,alpha2) N = length(a); %% N >=2 bz = zeros(1,(N-1)*2+1); az = zeros(1,(N-1)*2+1); aa = [1,alpha1,alpha2]; bb = [-alpha2,-alpha1,-1]; for i = 1:N x1 = 1;

for j = 1:(N-i)

x1 = conv(x1,aa); end x2 = 1;

for j = 1 : (i-1)

x2 = conv(x2,bb); end

az = az + a(i)*conv(x1,x2); bz = bz + b(i)*conv(x1,x2); end

bz = bz/az(1); az = az/az(1);

%% Transform lowpass digital filters to bandstop %% 2011-11-13 by nlp

function [bz,az] = lp2bsz(b,a,alpha1,alpha2) N = length(a); %% N >=2 bz = zeros(1,(N-1)*2+1); az = zeros(1,(N-1)*2+1); aa = [1,alpha1,alpha2]; bb = [alpha2,alpha1,1]; for i = 1:N x1 = 1;

for j = 1:(N-i)

x1 = conv(x1,aa); end x2 = 1;

for j = 1 : (i-1)

x2 = conv(x2,bb); end

az = az + a(i)*conv(x1,x2); bz = bz + b(i)*conv(x1,x2); end

bz = bz/az(1); az = az/az(1);

% Compute Rp(passband ripple) and As(stopband attenuation) of Digital % filter and plot the Megnitude Response of filter in the new figure % window. %

% [Rpp,Ass] = freqzn(num,den,wp,ws,Rp,As,ftype)

% returns the Rp(passband ripple) and As(stopband attenuation) of Digital filter %

% given numerator and denominator coefficients in vectors num and den.

% wp and ws are the passband and stopband edge frequencies,normalized form % 0 to 1(where 1 corresponds to pi radians/sample).For example, % lowpass: wp = 0.1 ws = 0.2 % highpass: wp = 0.2 ws = 0.1

% bandpass: wp = [0.2 0.7], ws = [0.1 0.8] % bandstop: wp = [0.1 0.8], ws = [0.2 0.7] % ftype: 'low' or 'high' or 'band' or 'stop' %

% 2011-11-20 by nlp

function [Rpp,Ass] = freqzn(num,den,wp,ws,Rp,As,ftype)

b = num; a = den; if As <= 35

As_low = As+10; else

As_low = 100; end

[H,w] = freqz(b,a,1000,'whole'); H = H(1:501)'; w = w(1:501)';

dB = 20*log10((abs(H)+eps)/max(abs(H))); rip = 10^(-Rp/20); att = 10^(-As/20); dw = 2/1000; figure;

switch ftype case 'low'

if As_low < 100

subplot(1,2,1),plot(w/pi,abs(H)); title('Lowpass Magnitude Response'); xlabel('Digital frequency in pi units'); ylabel('abs(H)'); axis([0,1,0,1.1]);

line([0,1],[rip,rip],'LineStyle',':'); line([0,1],[att,att],'LineStyle',':');

line([wp(1),wp(1)],[0,1.1],'LineStyle',':'); line([ws(1),ws(1)],[0,1.1],'LineStyle',':'); set(gca,'XTick',[0,wp(1),ws(1),1]); set(gca,'YTick',[0,att,rip,1]); subplot(1,2,2), end

plot(w/pi,dB);

title('Lowpass Magnitude Response in dB'); xlabel('Digital frequency in pi units'); ylabel('dB');

axis([0,1,-As_low,5]);

line([0,1],[-Rp,-Rp],'LineStyle',':'); line([0,1],[-As,-As],'LineStyle',':');

line([wp(1),wp(1)],[5,-As_low],'LineStyle',':'); line([ws(1),ws(1)],[5,-As_low],'LineStyle',':'); set(gca,'XTick',[0,wp(1),ws(1),1]); set(gca,'YTick',[-As_low,-As,-Rp,0]);

Rpp = -min(dB(1 : round(wp(1)/dw) +1));

Ass = -round(max(dB(round(ws(1)/dw) +1 : 501))); case 'high'

if As_low < 100

subplot(1,2,1),plot(w/pi,abs(H));

title('Highpass Magnitude Response'); xlabel('Digital frequency in pi units'); ylabel('abs(H)'); axis([0,1,0,1.1]);

line([0,1],[rip,rip],'LineStyle',':'); line([0,1],[att,att],'LineStyle',':');

line([wp(1),wp(1)],[0,1.1],'LineStyle',':'); line([ws(1),ws(1)],[0,1.1],'LineStyle',':'); set(gca,'XTick',[0,ws(1),wp(1),1]); set(gca,'YTick',[0,att,rip,1]);

subplot(1,2,2), end

plot(w/pi,dB);

title('Highpass Magnitude Response in dB'); xlabel('Digital frequency in pi units'); ylabel('dB');

axis([0,1,-As_low,5]);

line([0,1],[-Rp,-Rp],'LineStyle',':'); line([0,1],[-As,-As],'LineStyle',':');

line([wp(1),wp(1)],[5,-As_low],'LineStyle',':'); line([ws(1),ws(1)],[5,-As_low],'LineStyle',':'); set(gca,'XTick',[0,ws(1),wp(1),1]);

set(gca,'YTick',[-As_low,-As,-Rp,0]);

Rpp = -min(dB(round(wp(1)/dw) +1 : 501)); Ass = -round(max(dB(1 : round(ws(1)/dw) +1))); case 'band'

if As_low < 100

subplot(2,1,1),plot(w/pi,abs(H));

title('Bandpass Magnitude Response'); xlabel('Digital frequency in pi units'); ylabel('abs(H)'); axis([0,1,0,1.1]);

line([0,1],[rip,rip],'LineStyle',':'); line([0,1],[att,att],'LineStyle',':');

line([wp(1),wp(1)],[0,1.1],'LineStyle',':'); line([ws(1),ws(1)],[0,1.1],'LineStyle',':'); line([wp(2),wp(2)],[0,1.1],'LineStyle',':'); line([ws(2),ws(2)],[0,1.1],'LineStyle',':');

set(gca,'XTick',[0,ws(1),wp(1),wp(2),ws(2),1]); set(gca,'YTick',[0,att,rip,1]); subplot(2,1,2), end

plot(w/pi,dB);

title('Bandpass Magnitude Response in dB'); xlabel('Digital frequency in pi units'); ylabel('dB');

axis([0,1,-As_low,5]);

line([0,1],[-Rp,-Rp],'LineStyle',':'); line([0,1],[-As,-As],'LineStyle',':');

line([wp(1),wp(1)],[5,-As_low],'LineStyle',':'); line([ws(1),ws(1)],[5,-As_low],'LineStyle',':'); line([wp(2),wp(2)],[5,-As_low],'LineStyle',':'); line([ws(2),ws(2)],[5,-As_low],'LineStyle',':');

set(gca,'XTick',[0,ws(1),wp(1),wp(2),ws(2),1]); set(gca,'YTick',[-As_low,-As,-Rp,0]);

Rpp = -min(dB(round(wp(1)/dw) +1 : round(wp(2)/dw) +1)); Ass1 = -round(max(dB(1 : round(ws(1)/dw) +1))); Ass2 = -round(max(dB(round(ws(2)/dw) +1 : 501))); Ass = min(Ass1,Ass2); case 'stop'

if As_low < 100

subplot(2,1,1),plot(w/pi,abs(H));

title('Bandstop Magnitude Response'); xlabel('Digital frequency in pi units'); ylabel('abs(H)'); axis([0,1,0,1.1]);

line([0,1],[rip,rip],'LineStyle',':'); line([0,1],[att,att],'LineStyle',':');

line([wp(1),wp(1)],[0,1.1],'LineStyle',':'); line([ws(1),ws(1)],[0,1.1],'LineStyle',':'); line([wp(2),wp(2)],[0,1.1],'LineStyle',':'); line([ws(2),ws(2)],[0,1.1],'LineStyle',':');

set(gca,'XTick',[0,wp(1),ws(1),ws(2),wp(2),1]); set(gca,'YTick',[0,att,rip,1]); subplot(2,1,2), end

plot(w/pi,dB);

title('Bandstop Magnitude Response in dB'); xlabel('Digital frequency in pi units'); ylabel('dB');

axis([0,1,-As_low,5]);

line([0,1],[-Rp,-Rp],'LineStyle',':'); line([0,1],[-As,-As],'LineStyle',':');

line([wp(1),wp(1)],[5,-As_low],'LineStyle',':'); line([ws(1),ws(1)],[5,-As_low],'LineStyle',':'); line([wp(2),wp(2)],[5,-As_low],'LineStyle',':'); line([ws(2),ws(2)],[5,-As_low],'LineStyle',':'); set(gca,'XTick',[0,wp(1),ws(1),ws(2),wp(2),1]); set(gca,'YTick',[-As_low,-As,-Rp,0]);

Rpp1 = -min(dB(1:round(wp(1)/dw)+1)); Rpp2 = -min(dB(round(wp(2)/dw)+1:501)); Rpp = max(Rpp1,Rpp2);

Ass = -round(max(dB(round(ws(1)/dw)+1:round(ws(2)/dw)+1))); end

“数字信号处理”实验指导书（三）

一、实验课程编码：103044 二、实验课程名称：数字信号处理A

三、实验项目名称：应用MATLAB设计FIR数字滤波器及数字滤波器结构的MATLAB实现 四、实验目的

掌握应用MATLAB实现IIR、FIR数字滤波器结构的方法，即熟悉应用MATLAB实现IIR、FIR数字滤波器结构的相关函数。

掌握应用MATLAB设计FIR数字滤波器的方法，即熟悉应用MATLAB设计FIR滤波器相关的函数。 五、主要设备

安装有MATLAB软件的电脑 六、实验内容

编写MATLAB程序，实现以下各题：

12z1z21．已知某系统的系统函数为H(z)，求其级联型、并联型2411.76z6.25z结构。

4(1z1)(11.4z1z2)2．已知某系统的系统函数为H(z)，求其直接112(10.5z)(10.9z0.81z)型结构。

14.7512.90z124.5026.82z13．已知某系统的系统函数为H(z)，求其直

713211zz1z1z28322接型结构。

4．分别用N＝21和N＝51的Blackman窗对截频ωc＝0.4π的理想低通滤波器截断设计FIR滤波器并分别求出两个滤波器的幅度响应和最小阻带衰减As。从中

能得出什么结论。

5．用Hamming设计一个满足下列指标的线性相位FIR高通滤波器：

ωs＝0.4πrad，ωp＝0.6πrad， As= 45dB, Rp =1dB

画出所设计滤波器的增益响应和求出实际的As和 Rp。

6．分别用Blackman窗和Hamming窗设计满足下列指标的线性相位的FIR低通滤波器：

ωp＝0.4πrad，ωs＝0.6πrad， As= 45dB, Rp =0.5dB

画出所设计滤波器的幅频响应。简单评述两种窗的设计结果。 7．设计一个FIR数字带阻滤波器，使其满足：

p10.3, Rp10.5 dBs20.6, As240 dB，

s10.4, As140 dBp20.7, Rp20.5 dB七、实验步骤

1、熟悉与离散信号频谱分析相关的MATLAB函数（参考附录1）； 2、通过运行附录2中提供的例题，熟悉用MATLAB设计FIR数字滤波器的基本方法；

3、根据“六、实验内容”中各个题目的要求，编写MATLAB程序代码，调试程序，分析并保存结果。 八、实验结果

对实验练习题编写MATLAB程序并运行，在计算机上输出仿真结果。 附录1 主要的相关MATLAB函数

1．部分转换滤波器结构的文件 （1）tf2sos.m

调用格式：[sos,g] = tf2sos(b,a); （2）residuez.m

调用格式：[r,p,k] = residuez(num,den); [num,den] = residuez(r,p,k);

（3）conv.m

调用格式：c = conv(a,b); （4）poly.m

调用格式：p = poly(r); （5）filter.m

调用格式：y = filter(b,a,x); 2．部分窗函数文件

（1）rectwin.m / boxcar.m (矩形窗）

调用格式：w＝blackman ( N )

文件（2）～（5）的调用格式与文件（1）的调用格式类似。 （2）triang.m（三角窗 / 巴特利特窗） （3）hanning.m（汉宁窗） （4）hamming.m（海明窗） （5）blackman .m（布莱克曼窗） （6）kaiser.m（凯泽窗）

调用格式：w＝ kaiser ( N, beta )

3．部分FIR滤波器设计的文件 （1）fir1.m 窗函数设计法

调用格式：h＝fir1 (M, Wn, ‘type’, window) （2）fir2.m 频率采样设计法

调用格式：h＝fir2（M，F, A） （3）firpm.m 最优等波动设计法

调用格式：h＝firpm（M，F, A）

附录2 例题

13z111z227z318z4例1．将系统函数H(z)转换为级联

116z112z22z34z4z5型结构。

[MATLAB程序]： b = [1,-3, 11, 27, 18, 0]; a = [1, 16, 12, 2, -4, -1]; [sos, g] = tf2sos(b, a)

[运行结果]： sos =

1.0000 0 0 1.0000 15.2214 0 1.0000 -4.7366 18.2390 1.0000 -0.2477 -0.1246 1.0000 1.7366 0.9869 1.0000 1.0263 0.5272 g = 1

所以该系统函数的级联型结构为：

114.7366z118.239z211.7366z10.9869z2H(z)112115.2214z10.2477z0.1246z11.0263z10.5272z2

13z111z227z318z4例2．将系统函数H(z)转换为并联

116z112z22z34z4z5型结构。[MATLAB程序]： b = [1,-3, 11, 27, 18, 0]; a = [1, 16, 12, 2, -4, -1]; [r,p,k] = residuez (b, a)

[运行结果] r =

1.3030 0.3356 + 0.7215i 0.3356 - 0.7215i 2.2385 -3.2127 p =

-15.2214 -0.5132 + 0.5137i -0.5132 - 0.5137i 0.4979 -0.2503 k = 0

说明，系统没有常数项，有五个极点，其中第2、3个为共轭极点，应该把2、3极点组合起来构成一个二阶基本阶，然后再把4、5极点组合起来构成一个二阶基本阶。

[MATLAB程序]： [b1,a1] = residuez(r(2:3),p(2:3),[]) [b2,a2] = residuez(r(4:5),p(4:5),[])

[运行结果]

b1 = 0.6712 -0.3969 0

a1 = 1.0000 1.0263 0.5272 b2 = -0.9741 2.1600 0 a2 = 1.0000 -0.2477 -0.1246 所以该系统函数的并联型结构为：

1.3030.67120.3969z10.97412.16z1H(z)112115.2214z11.0263z0.5272z10.2477z10.1246z2

例3．使用窗函数设计法设计一个FIR数字低通滤波器，使其满足：

sample21.5104(rad/sec)p21.5103(rad/sec)s2310(rad/sec)3

As50 dB[MATLAB程序]： omega_s = 2*pi*1.5*10^4;

omega_p = 2*pi*1.5*10^3; omega_st = 2*pi*3*10^3; As = 50;

wp = 2*pi*omega_p/omega_s; ws = 2*pi*omega_st/omega_s; deltaw = ws - wp;

%% 由阻带衰减 >50dB，选择hamming窗；hamming窗的带宽为 6.6*pi/N N0 = ceil(6.6*pi/deltaw);

N = N0 + mod(N0+1,2) %% 保证N为奇数 wc = (wp + ws)/2;

%% 得到理想单位脉冲响应 hd n = 0:N-1; n0 = (N-1)/2; m = n -n0 + eps; hd = sin(wc*m)./(pi*m);

%% 得到N点的hamming窗 wdham = hamming(N)'; %% 得到脉冲响应 h = hd.*wdham;

%% 画出幅频特性(dB)，并检测As是否满足设计指标 [Rpp,Ass] = freqzn(h,1,wp/pi,ws/pi,0.5,As,'low') %% 如果Ass <50 则N = N+2 重新计算 [运行结果]： N = 33 As = 46

不满足设计指标，所以N=N+2后重新计算： N = 35 As = 55

例4．设计一个FIR数字高通滤波器，使其满足：

p0.6, Rp0.5dBs0.4, As60dB

[MATLAB程序]： wp = 0.6*pi; ws = 0.4*pi; Rp = 0.5; As = 60; deltaw = wp - ws;

%% 由阻带衰减 >60dB，选择blackman窗；blackman窗的带宽为 11*pi/N N0 = ceil(11*pi/deltaw);

N = N0 + mod(N0+1,2) %% 保证N为奇数 wc = (wp + ws)/2;

%% 得到高通滤波器的理想单位脉冲响应 hd n = 0:N-1; n0 = (N-1)/2; m = n -n0 + eps;

hd = sin(pi*m)./(pi*m) - sin(wc*m)./(pi*m); %% 高通滤波器的理想单位脉冲响应 %% 得到N点的blackman窗 wdham = blackman(N)'; %% 得到脉冲响应 h = hd.*wdham;

%% 画出幅频特性(dB)，并检测As是否满足设计指标 [Rpp,Ass] = freqzn(h,1,wp/pi,ws/pi,Rp,As,'high') %% 如果As不满足设计指标 则N = N+2 重新计算 hh = fir1(N-1,wc/pi,'high',blackman(N)); [运行结果]： N = 55 Rp = 0.0039 As = 71

例5．设计一个FIR数字带通滤波器，使其满足：

s10.2, As160 dBp20.65, Rp21 dB，

p10.35, Rp11 dBs20.8, As260 dB [MATLAB程序]： wp1 = 0.35*pi; wp2 = 0.65*pi; ws1 = 0.2*pi; ws2 = 0.8*pi; Rp = 1; As = 60;

deltaw = min((wp1 - ws1),(ws2-wp2));

%% 由阻带衰减 >60dB，选择blackman窗；blackman窗的带宽为 11*pi/N N0 = ceil(11*pi/deltaw);

N = N0 + mod(N0+1,2) %% 保证N为奇数 %% 得到带通滤波器的理想单位脉冲响应 hd wc1 = (wp1 + ws1)/2; wc2 = (wp2 + ws2)/2; n = 0:N-1; n0 = (N-1)/2;

m = n -n0 + eps;

hd = sin(wc2*m)./(pi*m) - sin(wc1*m)./(pi*m); %% 带通 %% 得到N点的blackman窗 wdham = blackman(N)'; %% 得到脉冲响应 h = hd.*wdham;

%% 画出幅频特性(dB)，并检测As是否满足设计指标 [Rpp,Ass] = freqzn(h,1,[wp1,wp2]/pi,[ws1,ws2]/pi,Rp,As,'band') %% 如果Ass <60 则N = N+2 重新计算 hh = fir1(N-1,[wc1,wc2]/pi,'band',blackman(N)); [运行结果]： N = 75 Rpp = 0.0030 Ass = 75

例题中用到的自编函数为：

% Compute Rp(passband ripple) and As(stopband attenuation) of Digital % filter and plot the Megnitude Response of filter in the new figure

% window. %

% [Rpp,Ass] = freqzn(num,den,wp,ws,Rp,As,ftype)

% returns the Rp(passband ripple) and As(stopband attenuation) of Digital filter %

% given numerator and denominator coefficients in vectors num and den.

% wp and ws are the passband and stopband edge frequencies,normalized form % 0 to 1(where 1 corresponds to pi radians/sample).For example, % lowpass: wp = 0.1 ws = 0.2 % highpass: wp = 0.2 ws = 0.1

% bandpass: wp = [0.2 0.7], ws = [0.1 0.8] % bandstop: wp = [0.1 0.8], ws = [0.2 0.7] % ftype: 'low' or 'high' or 'band' or 'stop' %

% 2011-11-20 by nlp

function [Rpp,Ass] = freqzn(num,den,wp,ws,Rp,As,ftype)

b = num; a = den; if As <= 35

As_low = As+10; else

As_low = 100; end

[H,w] = freqz(b,a,1000,'whole'); H = H(1:501)'; w = w(1:501)';

dB = 20*log10((abs(H)+eps)/max(abs(H))); rip = 10^(-Rp/20); att = 10^(-As/20); dw = 2/1000; figure;

switch ftype case 'low'

if As_low < 100

subplot(1,2,1),plot(w/pi,abs(H)); title('Lowpass Magnitude Response'); xlabel('Digital frequency in pi units'); ylabel('abs(H)'); axis([0,1,0,1.1]);

line([0,1],[rip,rip],'LineStyle',':'); line([0,1],[att,att],'LineStyle',':');

line([wp(1),wp(1)],[0,1.1],'LineStyle',':'); line([ws(1),ws(1)],[0,1.1],'LineStyle',':');

set(gca,'XTick',[0,wp(1),ws(1),1]); set(gca,'YTick',[0,att,rip,1]); subplot(1,2,2), end

plot(w/pi,dB);

title('Lowpass Magnitude Response in dB'); xlabel('Digital frequency in pi units'); ylabel('dB');

axis([0,1,-As_low,5]);

line([0,1],[-Rp,-Rp],'LineStyle',':'); line([0,1],[-As,-As],'LineStyle',':');

line([wp(1),wp(1)],[5,-As_low],'LineStyle',':'); line([ws(1),ws(1)],[5,-As_low],'LineStyle',':'); set(gca,'XTick',[0,wp(1),ws(1),1]); set(gca,'YTick',[-As_low,-As,-Rp,0]);

Rpp = -min(dB(1 : round(wp(1)/dw) +1));

Ass = -round(max(dB(round(ws(1)/dw) +1 : 501))); case 'high'

if As_low < 100

subplot(1,2,1),plot(w/pi,abs(H));

title('Highpass Magnitude Response'); xlabel('Digital frequency in pi units'); ylabel('abs(H)'); axis([0,1,0,1.1]);

line([0,1],[rip,rip],'LineStyle',':'); line([0,1],[att,att],'LineStyle',':');

line([wp(1),wp(1)],[0,1.1],'LineStyle',':'); line([ws(1),ws(1)],[0,1.1],'LineStyle',':'); set(gca,'XTick',[0,ws(1),wp(1),1]); set(gca,'YTick',[0,att,rip,1]); subplot(1,2,2), end

plot(w/pi,dB);

title('Highpass Magnitude Response in dB'); xlabel('Digital frequency in pi units'); ylabel('dB');

axis([0,1,-As_low,5]);

line([0,1],[-Rp,-Rp],'LineStyle',':'); line([0,1],[-As,-As],'LineStyle',':');

line([wp(1),wp(1)],[5,-As_low],'LineStyle',':'); line([ws(1),ws(1)],[5,-As_low],'LineStyle',':'); set(gca,'XTick',[0,ws(1),wp(1),1]);

set(gca,'YTick',[-As_low,-As,-Rp,0]);

Rpp = -min(dB(round(wp(1)/dw) +1 : 501)); Ass = -round(max(dB(1 : round(ws(1)/dw) +1))); case 'band'

if As_low < 100

subplot(2,1,1),plot(w/pi,abs(H));

title('Bandpass Magnitude Response'); xlabel('Digital frequency in pi units'); ylabel('abs(H)'); axis([0,1,0,1.1]);

line([0,1],[rip,rip],'LineStyle',':'); line([0,1],[att,att],'LineStyle',':');

line([wp(1),wp(1)],[0,1.1],'LineStyle',':'); line([ws(1),ws(1)],[0,1.1],'LineStyle',':'); line([wp(2),wp(2)],[0,1.1],'LineStyle',':'); line([ws(2),ws(2)],[0,1.1],'LineStyle',':');

set(gca,'XTick',[0,ws(1),wp(1),wp(2),ws(2),1]); set(gca,'YTick',[0,att,rip,1]); subplot(2,1,2), end

plot(w/pi,dB);

title('Bandpass Magnitude Response in dB'); xlabel('Digital frequency in pi units'); ylabel('dB');

axis([0,1,-As_low,5]);

line([0,1],[-Rp,-Rp],'LineStyle',':'); line([0,1],[-As,-As],'LineStyle',':');

line([wp(1),wp(1)],[5,-As_low],'LineStyle',':'); line([ws(1),ws(1)],[5,-As_low],'LineStyle',':'); line([wp(2),wp(2)],[5,-As_low],'LineStyle',':'); line([ws(2),ws(2)],[5,-As_low],'LineStyle',':'); set(gca,'XTick',[0,ws(1),wp(1),wp(2),ws(2),1]); set(gca,'YTick',[-As_low,-As,-Rp,0]);

Rpp = -min(dB(round(wp(1)/dw) +1 : round(wp(2)/dw) +1)); Ass1 = -round(max(dB(1 : round(ws(1)/dw) +1))); Ass2 = -round(max(dB(round(ws(2)/dw) +1 : 501))); Ass = min(Ass1,Ass2); case 'stop'

if As_low < 100

subplot(2,1,1),plot(w/pi,abs(H));

title('Bandstop Magnitude Response'); xlabel('Digital frequency in pi units'); ylabel('abs(H)');

axis([0,1,0,1.1]);

line([0,1],[rip,rip],'LineStyle',':'); line([0,1],[att,att],'LineStyle',':');

line([wp(1),wp(1)],[0,1.1],'LineStyle',':'); line([ws(1),ws(1)],[0,1.1],'LineStyle',':'); line([wp(2),wp(2)],[0,1.1],'LineStyle',':'); line([ws(2),ws(2)],[0,1.1],'LineStyle',':');

set(gca,'XTick',[0,wp(1),ws(1),ws(2),wp(2),1]); set(gca,'YTick',[0,att,rip,1]); subplot(2,1,2), end

plot(w/pi,dB);

title('Bandstop Magnitude Response in dB'); xlabel('Digital frequency in pi units'); ylabel('dB');

axis([0,1,-As_low,5]);

line([0,1],[-Rp,-Rp],'LineStyle',':'); line([0,1],[-As,-As],'LineStyle',':');

line([wp(1),wp(1)],[5,-As_low],'LineStyle',':'); line([ws(1),ws(1)],[5,-As_low],'LineStyle',':'); line([wp(2),wp(2)],[5,-As_low],'LineStyle',':'); line([ws(2),ws(2)],[5,-As_low],'LineStyle',':'); set(gca,'XTick',[0,wp(1),ws(1),ws(2),wp(2),1]); set(gca,'YTick',[-As_low,-As,-Rp,0]);

Rpp1 = -min(dB(1:round(wp(1)/dw)+1)); Rpp2 = -min(dB(round(wp(2)/dw)+1:501)); Rpp = max(Rpp1,Rpp2);

Ass = -round(max(dB(round(ws(1)/dw)+1:round(ws(2)/dw)+1))); end

“数字信号处理”实验指导书（四）

一、实验课程编码：103044

二、实验课程名称：数字信号处理A

三、实验项目名称：应用MATLAB分析音频信号的时频特性、并进行低通滤波 四、实验目的

综合应用数字信号处理理论知识，分析实际信号的频谱特性；应用设计得到的滤波器对信号进行滤波；理解“采样频率”；理解离散频率和数字角频率、模拟角频率、模拟频率之间的关系。 五、主要设备

安装有MATLAB软件的电脑；耳机 六、实验内容

1. 读音频文件：

 用wavread()函数获取音频数据，返回参数应当包含采样频率Fs。  实验中的音频是双声道音频，包括左右声道。左右声道很接近，频谱

分析时，可以只分析一个声道的音频的时频特性。

2. 画出音频信号时域波形

 要求：横坐标为时间（单位为秒） 3. 播放音频文件

 目的：理解音频主观听觉质量和音频频谱特性之间的关系 4. 对音频信号进行频谱分析：

1）直接对整个音频文件进行DFT变换，画出幅度谱，要求横坐标为模拟频率（kHz）。

2）对音频信号分帧，加窗，再进行DFT变换。要求：  每1024个样本为一帧。  加hanning窗

 每向前推移一帧有50%的重叠

 对每一帧进行DFT变换，以观察当前时刻的频谱特性  画出频谱图（幅度谱），横坐标为时间（单位为秒），纵坐标为频率

（kHz）

提示：第2种方法中画出的频谱图是一个二维图，参见附录2

5. 对实验中给出的有“毛刺”的音频信号进行低通滤波

 结合第4步频谱图，确定低通滤波器的通带截止频率和阻带起始频

率，通带波纹和阻带波纹可结合滤波结果自行调节设定。  采用Butterworth滤波器设计。

 对受损音频信号进行低通滤波，要求根据差分方程编程实现。

 用wavwrite ()函数，把滤波后的信号保存为wav格式的音频文件。

七、实验步骤

1、熟悉相关的MATLAB函数（参考附录1）；

2、根据“六、实验内容”中的要求，编写MATLAB程序代码，调试程序，分析并保存结果。

八、实验结果

编写MATLAB程序实现上述功能，在计算机上输出仿真结果。

附录1 主要的相关MATLAB函数

wavread( ) % 读wav格式的音频文件 wavwrite () % 写wav格式的音频文件

length() % 获取一维数组的长度

sound( ) % 播放音频文件 fft()

plot( ) % 画一维曲线 imagesc() % 画二维图形

xlabel( ) % 标明横坐标内容 ylabel( ) % 标明纵坐标内容

hanning() ％ hanning窗

buttord() % 获取Butterworth滤波器的阶数和截至频率

butter() % 获取Butterworth滤波器的分子多项式和分母多项式

的系数

附录2 时频曲线示例