数字信号处理实验指导书.docx
- 文档编号:4382652
- 上传时间:2022-12-01
- 格式:DOCX
- 页数:48
- 大小:388.36KB
数字信号处理实验指导书.docx
《数字信号处理实验指导书.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验指导书.docx(48页珍藏版)》请在冰豆网上搜索。
数字信号处理实验指导书
“数字信号处理”实验指导书
(一)
一、实验课程编码:
105003
二、实验课程名称:
数字信号处理
三、实验项目名称:
应用MATLAB分析离散信号频谱
四、实验目的
掌握应用MATLAB分析离散信号频谱的方法,即熟悉应用MATLAB分析离散信号的函数。
五、主要设备
安装有MATLAB软件的电脑
六、实验内容
编写MATLAB程序,实现下面题目:
1.用快速卷积法计算下面两个序列的线性卷积。
,
2.已知序列
(1)计算该序列DTFT的表达式
,并画出N=10时的
曲线;
(2)计算N=10时,序列x[n]的10点DTF,并画出曲线。
(3)利用hold函数,将DTFT和DFT画在一起。
3.理解高密度频谱和高分辨率频谱的概念。
设
(1)取0≤n≤9,求X1(k)
(2)将
(1)中的序列补零加长到0≤n≤99,求X2(k);
(3)增加取样值的个数,取0≤n≤99,求X3(k);
(4)将(3)中的序列补零加长到0 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 实验 指导书