4．用DFT对连续信号做谱分析。
设
，用DFT分析
的频谱，选择不同的截取长度Tp，观察截断效应，并用加窗的方法减少谱间干扰。
选取的参数：
（1）采样频率fs=400Hz；
（2）加窗时选取两种窗函数：
矩形窗函数和Hamming窗函数。
（3）
为截取时间长度，选取四种截取长度0.04s、2×0.04s、4×0.04s、8×0.04s。
对四种截取结果分别加矩形窗和Hamming窗后做2048点DFT。
5.已知连续信号为
，利用DFT近似分析其频谱。
要求频率分辨率为1Hz，确定抽样频率
、持续时间Tp和抽样点数N。
说明：
连续信号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=41；
（2）计算N=64时的X（k）=DFT[x（n）]，并画出x（n）和X（k）
（3）对X（k）在
上进行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.mz
调用格式：
〔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.000017.0000-2.00007.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.000018.000021.00006.00009.0000
例2、基本序列的离散傅立叶变换计算
复正弦序列：
，余弦序列：
分别对以上序列求当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的物理意义
已知
，
，绘制相应的幅频和相频曲线，并计算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=ones（1,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）|'）;holdon;
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）]'）;holdon;
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）|'）;holdon;
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）]'）;holdon;
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型滤波器逼近，设计指标为：
画出所设计的滤波器的幅度响应。
5．使用双线性变换法设计带通数字滤波器，要求用切比雪夫I型滤波器逼近，设计指标为：
画出所设计的滤波器的幅度响应。
七、实验步骤
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为带通滤波器通带的上、下截止频率
低通-带阻
，
，
ω2，ω1为带阻滤波器通带的上、下截止频率
附录2例题
例1．设计一个低通模拟巴特沃斯滤波器，使其满足：
通带：
阻带：
[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（'MagnitudeResponse'）;
xlabel（'Analogfrequencyinpiunits'）;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（'MagnitudeindB'）;
xlabel（'Analogfrequencyinpiunits'）;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
0000001.2665
A=1.0e+029*
0.00000.00000.00000.00000.00000.00011.2665
例2．设计一个低通模拟切比雪夫I型滤波器，使其满足：
通带截至频率：
，通带波纹：
阻带起始频率：
，阻带衰减：
[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（'MagnitudeResponse'）;
xlabel（'Analogfrequencyinpiunits'）;
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（'MagnitudeindB'）;
xlabel（'Analogfrequencyinpiunits'）;
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*
00003.8286
a=1.0e+018*
0.00000.00000.00000.00024.2958
例3．利用脉冲响应不变法设计一个巴特沃斯数字低通滤波器，使其满足：
[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=0000000.1266
a=1.00002.73803.74843.25331.88240.69050.1266
bz=0.00000.00070.01050.01670.00420.00010
az=1.0000-3.34435.0183-4.21902.0725-0.56000.0647
Rpp=0.9202
Ass=15
例4．利用双线性变换法设计一个切比雪夫I型数字低通滤波器，使其满足：
[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,A
升级会员 