数字信号处理西电上机实验.docx
- 文档编号:8180191
- 上传时间:2023-01-29
- 格式:DOCX
- 页数:28
- 大小:190.51KB
数字信号处理西电上机实验.docx
《数字信号处理西电上机实验.docx》由会员分享,可在线阅读,更多相关《数字信号处理西电上机实验.docx(28页珍藏版)》请在冰豆网上搜索。
数字信号处理西电上机实验
数字信号处理实验报告
实验一:
信号、系统及系统响应
一、实验目的:
(1)熟悉连续信号经理想采样前后的频谱变化关系,加深对时域采样定理的理解。
(2)熟悉时域离散系统的时域特性。
(3)利用卷积方法观察分析系统的时域特性。
(4)掌握序列傅里叶变换的计算机实现方法,利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
二、实验原理与方法:
(1)时域采样。
(2)LTI系统的输入输出关系。
三、实验内容、步骤
(1)认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容,阅读本实验原理与方法。
(2)编制实验用主程序及相应子程序。
①信号产生子程序,用于产生实验中要用到的下列信号序列:
a.xa(t)=A*e^-at*sin(Ω0t)u(t)
A=444.128;a=50*sqrt
(2)*pi;
b.单位脉冲序列:
xb(n)=δ(n)
c.矩形序列:
xc(n)=RN(n),N=10
②系统单位脉冲响应序列产生子程序。
本实验要用到两种FIR系统。
a.ha(n)=R10(n);
b.hb(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)
③有限长序列线性卷积子程序
用于完成两个给定长度的序列的卷积。
可以直接调用MATLAB语言中的卷积函数conv。
conv用于两个有限长度序列的卷积,它假定两个序列都从n=0开始。
调用格式如下:
y=conv(x,h)
四、实验内容
调通并运行实验程序,完成下述实验内容:
①分析采样序列的特性。
a.取采样频率fs=1kHz,即T=1ms。
b.改变采样频率,fs=300Hz,观察|X(ejω)|的变化,并做记录(打印曲线);进一步降低采样频率,fs=200Hz,观察频谱混叠是否明显存在,说明原因,并记录(打印)这时的|X(ejω)|曲线。
②时域离散信号、系统和系统响应分析。
a.观察信号xb(n)和系统hb(n)的时域和频域特性;利用线性卷积求信号xb(n)通过系统hb(n)的响应y(n),比较所求响应y(n)和hb(n)的时域及频域特性,注意它们之间有无差别,绘图说明,并用所学理论解释所得结果。
b.观察系统ha(n)对信号xc(n)的响应特性。
③卷积定理的验证。
五、思考题
(1)在分析理想采样序列特性的实验中,采样频率不同时,相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同?
它们所对应的模拟频率是否相同?
为什么?
(2)在卷积定理验证的实验中,如果选用不同的频域采样点数M值,例如,选M=10和M=20,分别做序列的傅里叶变换,求得的结果有无差异?
为什么?
六、参考资料
数字信号处理教程——Matlab释义与实现第二版,陈怀琛著,电子工业出版社,2008年
七、实验程序及结果
本实验使用自定义函数的方法产生信号:
门函数:
function[y,n]=gate(np,ns,nf)
ifns>np|ns>nf|np>nf
error('ÊäÈëλÖòÎÊý²»Âú×ãns<=np<=nf');
else
n=[ns:
nf];
y=(n>=0&n end 冲击函数: function[y,n]=impseq(np,ns,nf) ifns>np|ns>nf|np>nf error('ÊäÈëλÖòÎÊý²»Âú×ãns<=np<=nf'); else n=[ns: nf]; y=[(n-np)==0]; end 时域函数: function[y,ts]=sig(t0,tp,t1); A=444.128;a=50.*sqrt (2).*pi;w0=50.*sqrt (2).*pi; ts=t0: tp: t1; y=A.*exp(-a.*ts).*sin(w0.*ts).*(ts>=0); ;hb=impseq(0,0,10)+2.5.*impseq(1,0,10)+2.5.*impseq(2,0,10)+impseq(3,0,10); ;[ha,N]=gate(10,0,10); 第一部分主程序: [y1,t1]=sig(0,0.001,1); %t1=0: 0.001: 1; %y1=0.5.*sin(2.*pi.*15.*t1); figure; plot(t1,y1); axis([00.15-10150]); title('f=1000ʱÓòÐźÅ'); figure; f=-500: 500; yy1=fftshift(fft(y1)); plot(f,abs(yy1)); title('f=1000ƵÆ×'); [y2,t2]=sig(0,1/300,1); figure; plot(t2,y2); axis([00.15-10150]); title('f=300ʱÓòÐźÅ'); figure; yy2=fftshift(fft(y2)); plot(-150: 150,abs(yy2)); title('f=300ƵÆ×'); [y3,t3]=sig(0,0.02,1); figure; plot(t3,y3); title('f=200ʱÓòÐźÅ'); axis([00.15-85]); figure; yy3=fftshift(fft(y3)); plot(-25: 25,abs(yy3)); title('f=200ƵÆ×'); 第二部分主程序: xb=impseq(0,0,100); figure; stem(0: length(xb)-1,xb); axis([-1,3,0,2]); figure; y1=fft(xb); plot(0: length(y1)-1,y1); hb=impseq(0,0,100)+2.5.*impseq(1,0,100)+2.5.*impseq(2,0,100)+impseq(3,0,100); figure; stem(0: length(hb)-1,hb); axis([-1,5,0,3]) figure y2=fftshift(fft(hb)); plot(-(length(y2)-1)/2: (length(y2)-1)/2,abs(y2)); figure; y3=conv(double(xb),double(hb)); plot(0: length(y3)-1,y3); y4=fftshift(fft(y3)); figure; plot(-(length(y4)-1)/2: (length(y4)-1)/2,abs(y4)); xc=gate(10,0,100); ha=xc; y5=fftshift(fft(conv(double(xc),double(ha)))); figure; plot(-(length(y5)-1)/2: (length(y5)-1)/2,abs(y5)); 实验结果: 第一部分: 采样频率为200Hz时时域恢复200Hz时频谱 采样频率为300Hz时时域恢复300Hz时频谱 采样频率为1000Hz是时域恢复1000Hz时频谱 第二部分: Xb频域分析xb时域分析 hb频域分析hb时域分析 信号通过xb频域信号通过xc频域 卷积定理的验证: %%%%%卷积定理的验证 y1=[0,1,2,3,4,5,6,7,8]; y2=[8,7,6,5,4,3,2,1,0]; y3=conv(y1,y2); f1=fft(y1,20); f2=fft(y2,20); f3=fft(y3,20); f=f1.*f2; subplot(2,1,1); stem(0: length(f)-1,f,'.'); subplot(2,1,2); stem(0: length(f3)-1,f3,'.'); 时域卷积后求频谱和频域相乘 可见,f=conv(y1,y2)的频谱和y1,y2的频谱相乘后结果相同。 即满足卷积定理 八、回答问题 (1)在分析理想采样序列特性的实验中,采样频率不同时,相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同? 它们所对应的模拟频率是否相同? 为什么? 答: 数字频率度量不相同,但他们所对应的模拟频率相同。 由w=Ω*Ts公式得,采样间隔变化时模拟频率对应的数字频率会有相应的变化,故其度量会有所变化。 而且采样频率的大小直接关系到能否将能否将原始信号恢复出来。 (2)在卷积定理验证的实验中,如果选用不同的频域采样点数M值,例如,选M=10和M=20,分别做序列的傅里叶变换,求得的结果有无差异? 答: 有差异,所到的结果点数不同。 九、总结 本实验主要是后续实验的基础。 涉及内容也比较浅显。 单位冲击序列与hb卷积后得到的频谱与hb原频谱相同。 原因很简单,是因为单位冲击序列卷积任何函数仍然是原函数。 实验二: 用FFT作谱分析 一、实验目的 (1)进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法,所以FFT的运算结果必然满足DFT的基本性质)。 (2)熟悉FFT算法原理和FFT子程序的应用。 (3)学习用FFT对连续信号和时域离散信号进行谱分析的方法,了解可能出现的分析误差及其原因,以便在实际中正确应用FFT。 二、实验步骤 (1)复习DFT的定义、性质和用DFT作谱分析的有关内容。 (2)复习FFT算法原理与编程思想,并对照DIT-FFT运算流图和程序框图,读懂本实验提供的FFT子程序。 (3)编制信号产生子程序,产生以下典型信号供谱分析用: (4)编写主程序。 下图给出了主程序框图,供参考。 本实验提供FFT子程序和通用绘图子程序。 (5)按实验内容要求,上机实验,并写出实验报告。 三、实验内容 (1)对2中所给出的信号逐个进行谱分析。 (2)令x(n)=x4(n)+x5(n),用FFT计算8点和16点离散傅里叶变换, X(k)=DFT[x(n)] (3)令x(n)=x4(n)+jx5(n),重复 (2)。 四、思考题 (1)在N=8时,x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢? (2)如果周期信号的周期预先不知道,如何用FFT进行谱分析? 五、参考资料 数字信号处理教程——Matlab释义与实现第二版,陈怀琛著,电子工业出版社,2008年 数字信号处理第三版,高西泉等著,西安电子科技大学出版社,第三版 六、实验程序及结果 x1=[1,1,1,1,0,0,0,0,0,0]; x2=[1,2,3,4,4,3,2,1,0,0]; x3=[4,3,2,1,1,2,3,4,0,0]; n=0: 1000; x4=cos(n); x5=sin(n); t=0: 0.001: 1; x6=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t); %figure; %title('Öð¸öÆ×·ÖÎö') subplot(3,2,1); stem(0: length(x1)-1,x1); subplot(3,2,2); stem(0: length(x2)-1,x2); subplot(3,2,3); stem(0: length(x3)-1,x3); subplot(3,2,4); stem(0: 20,x4(1: 21)); subplot(3,2,5); stem(0: 20,x5(1: 21)); subplot(3,2,6); stem(0: 20,x6(1: 21)); figure; stem(0: 200,x6(1: 201)); y1=abs(fft(x1,8)); y2=abs(fft(x2,8)); y3=abs(fft(x3,8)); y4=abs(fft(x4,8)); y5=abs(fft(x5,8)); y6=abs(fft(x6,8)); figure; subplot(3,2,1); stem(0: length(y1)-1,y1); subplot(3,2,2); stem(0: length(y2)-1,y2); subplot(3,2,3); stem(0: length(y3)-1,y3); subplot(3,2,4); stem(0: length(y1)-1,y1); subplot(3,2,5); stem(0: length(y5)-1,y5); subplot(3,2,6); stem(0: length(y6)-1,y6); %%%figure; %%%stem(0: 15,y6(1: 16)); %%%%%%%%%%%%%%% y1=abs(fft(x1,16)); y2=abs(fft(x2,16)); y3=abs(fft(x3,16)); y4=abs(fft(x4,16)); y5=abs(fft(x5,16)); y6=abs(fft(x6,16)); figure; subplot(3,2,1); stem(0: length(y1)-1,y1); subplot(3,2,2); stem(0: length(y2)-1,y2); subplot(3,2,3); stem(0: length(y3)-1,y3); subplot(3,2,4); stem(0: length(y6)-1,y4); subplot(3,2,5); stem(0: length(y6)-1,y5); subplot(3,2,6); stem(0: length(y6)-1,y6); %figure; %stem(0: 15,y6(1: 16)); x7=x4+j*x5; y71=abs(fft(x7,8)); y72=abs(fft(x7,16)); figure; subplot(1,2,1); stem(0: length(y71)-1,y71); subplot(1,2,2); stem(0: length(y72)-1,y72); 程序结果如下所示: 每个信号时域分析 因为第六个信号不易在上图中展示(包含信息不到一个周期)所以单独画出图形如下: x6的时域波形 逐个进行8点离散傅里叶变换 逐个进行16点离散傅里叶变换 x(n)8点和16点离散傅里叶变换比较 七、回答问题 (1)在N=8时,x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢? 答: N=8时幅频特性一样,N=16时幅频特性不一样。 Fn所能分辨到频率为为Fs/N,如果采样频率Fs为1024Hz,采样点数为1024点,则可以分辨到1Hz。 1024Hz的采样率采样1024点,刚好是1秒,也就是说,采样1秒时间的信号并做FFT,则结果可以分析到1Hz,如果采样2秒时间的信号并做FFT,则结果可以分析到0.5Hz。 如果要提高频率分辨力,则必须增加采样点数,也即采样时间。 频率分辨率和采样时间是倒数关系。 总结: 采样点N的不同,分辨率不同,在N=8,俩函数不同的部分没有分辨出来,因此特性相同,而N=16时已经足以分辨出不同的频率,因此频谱不同。 (2)如果周期信号的周期预先不知道,如何用FFT进行谱分析? 答: 设一个定长的m值,先取2m,看2m/m的误差是否大,如大的话再取4m,看4m/2m的误差是否大,如不大,4m(4倍的m值)则可近似原来点的谱分析。 八、总结 本实验主要掌握fft函数的用法,包括采样点的意义。 这个意义已经在七、回答问题中描述,不再赘述。 实验三: 用窗函数法设计FIR数字滤波器 一、实验目的 (1)掌握用窗函数法设计FIR数字滤波器的原理和方法。 (2)熟悉线性相位FIR数字滤波器特性。 (3)了解各种窗函数对滤波特性的影响。 二、实验内容及步骤 (1)复习用窗函数法设计FIR数字滤波器一节内容,阅读本实验原理,掌握设计步骤。 (2)编写程序。 ①编写能产生矩型窗、升余弦窗、改进升余弦窗和二阶升余弦窗的窗函数子程序。 ②编写主程序。 其中幅度特性要求用dB表示。 窗函数法设计滤波器主程序框图 三、实验内容 四、思考题 (1)如果给定通带截止频率和阻带截止频率以及阻带最小衰减,如何用窗函数法设计线性相位低通滤波器? 写出设计步骤。 (2)如果要求用窗函数法设计带通滤波器,且给定上、下边带截止频率为ω1和ω2,试求理想带通的单位脉冲响应hd(n)。 五、参考资料 数字信号处理教程——Matlab释义与实现第二版,陈怀琛著,电子工业出版社,2008年 数字信号处理第三版,高西泉等著,西安电子科技大学出版社,第三版 XX知道 六、实验程序及结果 myfreqz(b,a)函数,是自编的改进freqz()的函数。 function[db,mag,pha,grd,w]=myfreqz(b,a); %advancedfreqz %[db,mag,pha,grd,w]=myfreqz(b,a); [H,w]=freqz(b,a,1000,'whole'); H=(H(1: 1: 501))';w=(w(1: 1: 501))'; mag=abs(H); db=20*log10((mag+eps)/max(mag)); pha=angle(H); grd=grpdelay(b,a,w); 以下是主程序部分 wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp; N0=ceil(6.6*pi/deltaw); N=N0+mod(N0+1,2); wdham=(hamming(N))'; wc=(ws+wp)/2; %hd=ideallp(wc,N); tao=(N-1)/2; n=[0: (N-1)]; m=n-tao+eps; hd=sin(wc*m)./(pi*m); h=hd.*wdham; [db,mag,pha,grd,w]=myfreqz(h,[1]); dw=2*pi/1000; Rp=-(min(db(1: wp/dw+1))); As=-round(max(db(ws/dw+1: 501))); figure; subplot(2,2,1); stem(0: length(hd)-1,hd,'.'); subplot(2,2,2); stem(0: length(wdham)-1,wdham,'.'); subplot(2,2,3); stem(0: length(h)-1,h,'.'); subplot(2,2,4); plot(1/length(db): 1/length(db): 1,db); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% wn1=boxcar(N); figure; subplot(2,2,1); stem(0: N-1,wn1,'.'); title('boxcar'); wn2=bartlett(N); subplot(2,2,2); stem(0: N-1,wn2,'.');title('bartlett'); wn3=hanning(N); subplot(2,2,3); stem(0: N-1,wn3,'.');title('hanning'); %wn4=hamming(N); wn4=blackman(N); subplot(2,2,4); stem(0: N-1,wn4,'.');title('blackman'); %wn=(N,beta); 得到的结果如下图: 其中,图一是理想脉冲响应,第一行图二是汉明窗;图三是实际脉冲响应,图四是幅度响应。 一下是几个典型窗函数的图像。 ,汉明窗上例已画出,不再画。 每个窗的名字已在每个子图上方标出。 典型窗函数 七、回答问题 (1)给定通带截止频率和阻带截止频率以及阻带最小衰减,用窗函数法设计线性相位低通滤波器的设计步骤: 答: 技术指标Wp=0.2*pi,Ws=0.4*pi,Ap=0.25dB,As=50dB 方法一 选择海明窗 clearall; Wp=0.2*pi; Ws=0.4*pi; tr_wide=Ws-Wp;%过渡带宽度 N=ceil(6.6*pi/tr_wide)+1;%滤波器长度 n=0: 1: N-1; Wc=(Wp+Ws)/2;%理想低通滤波器的截止频率 hd=ideal_lp1(Wc,N);%理想滤波器的单位冲击响应 w_ham=(hamming(N))';%海明窗 h=hd.*w_ham;%实际海明窗的响应 [db,mag,pha,w]=freqz_m2(h,[1]);%计算实际滤波器的幅度响应 delta_w=2*pi/1000; Ap=-(min(db(1: 1: Wp/delta_w+1)))%实际通带纹波 As=-round(max(db(Ws/delta_w+1: 1: 501)))%实际阻带纹波 subplot(221) stem(n,hd) title('理想单位脉冲响应hd(n)') subplot(222) stem(n,w_ham) title('海明窗') subplot(223) stem(n,h) title('实际单位脉冲响应hd(n)') subplot(224) plot(wi/pi,db) title('幅度响应(dB)') axis([0,1,-100,10]) 方法二 Window=blackman(16); b=fir1(15,0.3*pi,'low',Window); freqz(b,128) (2)用窗函数法设计带通滤波器,且给定上、下边带截止频率为ω1和ω2,理想带通的单位脉冲响应hd(n)的求解过程: 答: 由hd(n)=1/2π∫w1woe-jwaejnwdw+1/2π∫w2w3e-jwaejnwdw (其中w0=-w0-wc,w1=-w0+wc,w2=w0-wc,w3=w0+wc) 计算整理后可得: hd(n)=2/((n-a)*π)*sin[(n-a)wc]*cos[(n-a)w0] =2wc/π*sa[(n-a)wc]*cos[(n-a)w0] 实验四: 用双线性变换法设计IIR数字滤波器 一、实验目的 (1)掌握用双线性变换法设计IIR数字滤波器的原理和方法。 (2)熟悉IIR数字滤波器特性。 二、实验内容及步骤 (1)复习用双线性变换法设计IIR数字滤波器一节内容,阅读本实验原理,掌握设计步骤。 (2)编写程序。 ①设计产生巴特沃斯滤波器原型,采用双线性变换法设计巴特沃斯数字低通滤波器。 ②编写主程序。 其中幅度特性要求用dB表示。 三、上机实验内容 (1)实验所需设计的IIR数字低通滤波器的指标为wp=0.2*pi,ws=0.3*pi,Rp=1dB,As=15dB,滤波器
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 上机 实验