数字信号处理上机报告.docx

《数字信号处理上机报告.docx》由会员分享，可在线阅读，更多相关《数字信号处理上机报告.docx（20页珍藏版）》请在三一办公上搜索。

1、数字信号处理上机报告数字信号处理上机报告 实验一：信号、系统及系统响应 一、实验目的 (1) 熟悉连续信号经理想采样前后的频谱变化关系， 加深对时域采样定理的理解。 (2) 熟悉时域离散系统的时域特性。 (3) 利用卷积方法观察分析系统的时域特性。 (4) 掌握序列傅里叶变换的计算机实现方法， 利用序列的傅里叶变换对连续信号、 离散信号及系统响应进行频域分析。 二、实验原理 (1) 时域采样定理 模拟信号经过采样脉冲进行采样，得到离散时间信号。若采用周期脉冲进行等间隔采样，则得到的离散时间信号的频谱是原信号频谱的周期延拓。 假设采样脉冲为：p(t)=原模拟信号是x(t) 则采样输出信号：x(t

2、)n=-d(t-nT) +=x(t)p(t) 根据采样脉冲的周期性，进行傅里叶级数展开，得到p(t)=k=-cek+jkWst1其中：ck=T1c=解得：kTX(jW)=1=T2T2-Tp(t)e-jkWstdt 1+jkWste 则p(t)=Tk=-因此，采样输出信号的频谱为 +-x(t)p(t)e-jWtdt=+-+1x(t)Tk=-+ejkWste-jWtdtk=-+-1x(t)e-j(W-Ws)tdt=Tk=-X(jW-jkWs) 输出信号频谱为原模拟信号频谱周期延拓，延拓周期为T=12p=fsWsfs 为采样周期。 要使采样后信号频谱包含原信号频谱而不发生混叠，要满足fs=2fmax

3、 。 1 fmax 为原信号频谱的最大值。 (2) LTI系统的输入输出关系 LTI系统满足线性和时不变性，如果系统输入x(n) ,输出为y(n) ，系统函数为h(n) ，则满足关系y(n)=h(n)*x(n)=m=-h(m)x(n-m) +用卷积定理表示：Y(ejw)=X(ejw)H(ejw) DFT与DTFT的关系 当有限长序列的长度很大时，直接计算DTFT的计算量很大。而计算DFT可采用FFT，大大减轻计算量，加快程序运行速度。DFT是DTFT在0,2p) 上的等间隔采样，即 x(k)=x(n)e-jwn|n=0N2pw=N=x(ejw)|w=2pN 0kN-1 通过序列补零增大N，则会

4、使x(k)更多地反映x(ejw)的信息。如果N足够大，x(k)近似逼近x(ejw)。通过FFT就可以近似求得x(ejw)。 三、实验内容： 1、分析采样序列的特性 、实验条件 连续模拟信号：x(t)=Ae-at sin(0t)u(t) 其中：A=444.128，a=50*sqrt(2)*pi; 0=50*sqrt(2)*pi; 实验步骤： j取采样频率fs=1 kHz，观察|X(e)|的变化， 并做记录；进一步降低采样频率， 取fs=300 Hz，fs=200 Hz，观察频谱混叠现象，并记录这时的|X(e实验结果： j)|曲线，分析异同。 fs=1KHz 2 fs=300Hz fs=200Hz

5、 3 (4)实验结论： 通过比较不同频率采样的输出信号幅频特性，可以发现当采样频率较大时，输出信号频谱不发生混叠，随着采样频率的降低，输出信号频谱越来越胖，下凹部分越来越尖，当f=300HZ和200HZ时已经发生了混叠。 2、时域离散信号、 系统和系统响应分析 、实验条件 单位脉冲序列： xb(n)=(n) 矩形序列： 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) 、实验步骤 1、观察信号xb(n)和系统hb(n)的时域和频域特性； 2、利用线性卷积求信号xb(n)通

6、过系统hb(n)的响应y(n)， 比较所求响应y(n)和hb(n)的时域及频域特性， 注意它们之间的差别，并绘图。 同时观察系统ha(n)对信号xc(n)的响应特性。 3、卷积定理的验证。 直接调用MATLAB语言中的卷积函数conv得到线性卷积结果y=conv (x, h)。通过对y、x、h进行DTFT运算，验证Y(e)=X(e)H(e) 、实验结果 信号的时域、频域分析和系统输出 jwjwjw 4 系统ha(n)对xc(n)的响应特性 卷积定理的验证 5 (4)实验结论： 由图可知，y(n)和hb(n)的时域图像和幅频特性几乎一致，从而验证了x(n)=x(n)* (n). 通过程序直接计算

7、y(n)=ha(n)*xc(n)和频域计算Y(ejw)=Xc(ejw)Ha(ejw)，对比y(n)的幅度谱可以发现二者几乎一致。从而验证了卷积定理。 思考题： (1) 在分析理想采样序列特性的实验中， 采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同? 它们所对应的模拟频率是否相同? 为什么? 答：采样频率不同时，相应数字频率度量相同，其频谱都是0,2p) 的周期函数。而相同数字频率对应的模拟频率不同，有关系W=wTs 。w 相同，而Ts 不同，对应模拟频率不同。 (2) 在卷积定理验证的实验中， 如果选用不同的频域采样点数M值， 例如， 选M=10和M=20， 分别做

8、序列的傅里叶变换， 求得的结果有无差异？为什么？ 答：有差别。DFT相当于序列频谱的等间隔采样，当取点少时，DFT包含序列频谱的信息会少，与序列频谱的误差会增大。取点多时，DFT会反映更多的频谱信息，误差会小，减轻了栅栏效应。二者因为误差而有差别。 实验程序： clear all;clc; A=444.128; a=50*sqrt(2)*pi; om=50*sqrt(2)*pi; Y=1/1000;%采样周期 t=0:.0001:.05-.0001; 6 n=0:1:max(t)/Y;%从零开始0,0.05) t和n分离 p=length(n); N=1024; %FFT点数 xa=A.*ex

9、p(-1)*a*t).*sin(om*t);%实验序列 xb=dert(p,0);%单位脉冲响应 xc=ones(1,10) zeros(1,p-10);%矩形函数10 ha=xc; hb=xb+2.5.*dert(p,1)+2.5.*dert(p,2)+dert(p,3); %产生采样序列 xa_a=A.*exp(-1)*a*Y*n).*sin(om*Y*n); figure(1); subplot(1,3,1) plot(t,xa,b-);grid on;%绘制图形 title(原序列); subplot(1,3,2) stem(n,xa_a,g-);grid on; title(采样后序

10、列); c=fft(xa_a,N); c=abs(c); c=c c; t=1:2*N; subplot(1,3,3); plot(t,c,r-);grid on; title(采样后序列幅度谱); figure(2) subplot(3,2,1) stem(n,xb,g-); c=fft(xb,N); c=abs(c); t=1:N; subplot(3,2,2) plot(t,c); title(Xb频谱幅度特性) subplot(3,2,3) stem(n,hb,r-); c=fft(hb,N); c=abs(c); subplot(3,2,4) t=1:N; plot(t,c); ti

11、tle(Hb频谱幅度特性); 7 %输出响应yn yn=conv(xb,hb); subplot(3,2,5) n1=0:1:2*p-2; stem(n1,yn); title(输出信号); c=fft(yn,2*N); c=abs(c); t=1:2*N; subplot(3,2,6) plot(t,c); title(输出信号幅度谱); yn=conv(xc,ha); figure(3) subplot(1,3,1) stem(n,xc); title(Xc的波形); subplot(1,3,2) stem(n1,yn); c=fft(yn,2*N); c=abs(c); subplot(

12、1,3,3) plot(t,c); title(输出幅度谱); %验证卷积定理 figure(4) subplot(1,2,1) plot(t,c); title(时域处理结果); xcf=fft(xc,2*N); haf=fft(ha,2*N); xcf=abs(xcf); haf=abs(haf); ynf=xcf.*haf; subplot(1,2,2) plot(t,ynf); title(频域处理结果); hold off clear c haf xcf; 子函数： function y=dert(a,b)%a为n大小,b为位置 y=zeros(1,b) 1 zeros(1,(a-b

13、)-1); 8 end 实验二、用FFT作谱分析 一、实验目的 (1) 进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法， 所以FFT的运算结果必然满足DFT的基本性质)。 (2) 熟悉FFT算法原理和FFT子程序的应用。 (3) 学习用FFT对连续信号和时域离散信号进行谱分析的方法， 了解可能出现的分析误差及其原因， 以便在实际中正确应用FFT。 二、实验内容 利用FFT进行谱分析 实验条件 用到的信号： px4(n)=cosn x1(n)=R4(n)4 n+1,0n3p x5(n)=sinnx2(n)=8-n4n78 x6(n)=cos8pt+cos16pt+c

14、os20pt 0 4-n0n3 x3(n)=n-34n7 0实验步骤 1、对所给的信号逐个进行谱分析。 2、令x(n)=x4(n)+x5(n)， 用FFT计算 8 点和 16 点离散傅里叶变换，得X(k)=FFTx(n) 3、令x(n)=x4(n)+jx5(n)， 重复(2)。 (3)实验结果 1、逐个信号谱分析 9 10 11 2、x(n)=x4(n)+x5(n)的DFT 3、x(n)=x4(n)+jx5(n)的DFT 实验结论 利用FFT可以快速地进行序列的谱分析。 思考题 (1) 在N=8时， x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢? N=8时，序列的幅频特性相

15、同，因为序列x2(n)是1，2，3，4，4，3，2，1正好8个点。 x3(n)是4，3，2，1，1，2，3，4八个点，序列x2(8x3(n)，因此根据DFT循 12 -j4K=1 ，得幅频特性相同。而N=16时，通过对x2(n) 和x3(n)末端环移位特性，并且e补零，两个序列不能通过移位相等，因此幅频特性不同。 (2) 如果周期信号的周期预先不知道， 如何用FFT进行谱分析? 可以预先设定一个较小的N，进行N点FFT，再取2N，进行2N点FFT，如果二个FFT相等，则周期可能是N，直接取此结果作为FFT。，如果不等，取4N进行FFT，并与2N点FFT比较，依次进行，取2N 直到FFT近似相等

16、为止，求的近似的FFT。 实验程序： clear all;clc; %产生所需函数 n=0:1:39; %共40个点，时域取样0,1) p=length(n); t=1/p; %单独分析x6，时域对应0,1),采样40点. N=30; %FFT点数 x1=ones(1,4) zeros(1,p-4); x2=1,2,3,4 4,3,2,1 zeros(1,p-8); x3=4,3,2,1 1,2,3,4 zeros(1,p-8); x4=cos(pi/4).*n); x5=sin(pi/8).*n); x6=cos(8*pi*t).*n)+cos(16*pi*t).*n)+cos(20*pi*

17、t).*n); x=x1;x2;x3;x4;x5;x6; %进行谱分析 k=zeros(6,N); for i=1:6 temp=x(i,:); tmp=fft(temp,N); k(i,:)=abs(tmp); end %绘制波形 %产生频域坐标 ki=0:N-1; %ki=0:N-1*(length(n)/2)/N; for i=1:6 figure(i); subplot(1,2,1); temp=x(i,:); stem(n,temp);grid on; title(原x,int2str(i),序列); subplot(1,2,2); tmp=k(i,:); stem(ki,tmp);

18、grid on; title(FFT序列); end M 13 x=x4+x5; c=fft(x,8); c=abs(c); ki=0:7; figure(7) subplot(1,2,1) stem(ki,c); title(X4+X5的8点DFT); c=fft(x,16); c=abs(c); ki=0:15; subplot(1,2,2) stem(ki,c); title(16点DFT); x=x4+j.*x5; c=fft(x,8); c=abs(c); ki=0:7; figure(8); subplot(1,2,1) stem(ki,c); title(X4+jX5的8点DFT

19、); c=fft(x,16); c=abs(c); ki=0:15; subplot(1,2,2) stem(ki,c); title(16点DFT); clear temp tmp p ki k; 实验三：用窗函数法设计FIR数字滤波器 一、实验目的 (1) 掌握用窗函数法设计FIR数字滤波器的原理和方法。 (2) 熟悉线性相位FIR数字滤波器特性。 (3) 了解各种窗函数对滤波特性的影响。 二、实验原理 1、FIR滤波器 FIR滤波器是指单位脉冲响应是有限长的滤波器。可以具有线性相位，而滤波器也一定是稳定的。如果系统函数表示为H(z)=h(n)zn=0N-1-n则存在N-1个极点均位于z=

20、0.FIR滤波 14 器的频率响应可以表示成H(ejw)=H(w)ejq(w) 。具有线性相位的FIR滤波器有四种情况。 N为奇数，h(n)偶对称。此时H(w) 关于0和2偶对称，偶对称。 N为偶数，h(n)偶对称。此时H(w) 关于0和2偶对称，奇对称。 N为奇数，h(n)奇对称。此时H(w) 关于0和2奇对称，奇对称。 N为偶数，h(n)奇对称。此时H(w) 关于0和2奇对称，偶对称。 2、窗函数设计法 利用窗函数截取无限长的线性相位滤波器，可以得到FIR滤波器，利用此方法，可以先设计满足要求的无限长线性相位滤波器，再选择合适的窗函数截取。时域表现为无限长序列变为有限长，频域则反映为抖动，

21、窗函数不同抖动不同，过渡带宽度不同。 设计FIR滤波器，先根据所需设计理想线性相位滤波器 然后选择合适的窗函数w(n) 。再次确定理想线性滤波器的频率响应H(ejw) ，并求解出单位脉冲响应h(n) 。最后加窗得到h(n)=h(n)w(n)。 三、实验内容 希望对输入低频模拟信号做数字低通滤波器处理，以提取所需要的信号。设系统的采样频率为50kHz，要求通带截止频率为10kHz，阻带截止频率为25kHz，阻带最小衰减为60dB。 1）用窗函数法设计FIR数字低通滤波器，选择合适的窗函数以及滤波器长度N，求出单位脉冲响应h(n)。 2）、画出该FIR数字低通滤波器的幅频响应特性曲线和相频响应特性

22、曲线。 (1)窗函数选择： 根据不同窗函数的阻带最小衰减不同，确定满足条件的窗函数为布莱克曼窗和凯塞窗 (2)设计步骤 先设计理想线性相位低通滤波器。H(e)=H(w)ejwjq(w)H(w)=1 wwc H(w)=0 wcwp 根据要求确定wc和N，确定窗函数的阶次。 求的理想低通滤波器的单位脉冲响应为h(n)=sinwc(n-a)N-1 a= 2p(n-a)直接得到FIR滤波器的单位脉冲响应h(n)=h(n)w(n)。 15 (3)实验结果 1、布莱克曼窗设计 未加窗的理想滤波器单位脉冲图像 由于仅取有限长度，相频特性不理想。 FIR滤波器特性分析 2、凯塞窗设计 FIR滤波器单位脉冲响应

23、。 16 实验结论 通过窗函数设计法，可以得到严格线性相位FIR低通滤波器。 思考题 (1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减， 如何用窗函数法设计线性相位低通滤波器? 写出设计步骤。 答：根据通带截止频率和阻带截止频率，算出截止频率wc= 并设计理想低通线性相位滤波器，得到h(n) 依据阻带最小衰减确定窗函数类型。 1(ws-wp) 和过渡带宽 2 17 根据过渡带宽确定窗大小N，得到窗函数w(n) 。 最后由h(n)=h(n)w(n) 得到FIR低通滤波器的单位脉冲响应。 (2) 如果要求用窗函数法设计带通滤波器， 且给定上、 下边带截止频率为1和2，试求理想带通的单位脉

24、冲响应hd(n)。 答：设计理想带通线性相位滤波器，带宽B=w1-w2 。则频谱可知 H(w)=1 w2ww1 H(w)=0 其他 q(w)=-wN-1 N为常数。 2sin则hd(n)=2w1-w22p(n-a)(n-a)cos(w1+w22n) 其中a=N-1 2实验程序一：布莱克曼窗 clear all;clc; ws=0.5*2*pi; wp=0.2*2*pi; B=ws-wp; N=ceil(11*pi/B); wc=(wp+ws)/2;%截止频率 n=1:N; len=N; a=(N-1)/2; hd=sin(wc.*(n-a)./(pi.*(n-a); if mod(N,2)=1

25、 hd(1,a)=wc/pi; end figure(1) subplot(1,2,1) stem(n,hd); title(原函数时域图像); h=zeros(1,len); beta=4;%凯塞窗参数beta。此项无用 num=5;%窗函数类型 w=windows(N,num,beta); h=hd.*w; %计算DFT 18 FF=1024;%FFT点数 p=1:FF; HFF=zeros(1,FF);%幅度矩阵 HF=HFF;%相位矩阵 pp=fft(hd,FF); subplot(1,2,2) plot(p,angle(pp),r-); title(原函数相频特性); num=5;

26、temp=h; figure(num+1) subplot(1,3,1) stem(n,temp,b-); title(加窗后的函数,num2str(num); temp=fft(temp,FF); HF=abs(angle(temp); HFF=abs(temp); HFF=20.*log(HFF); subplot(1,3,2) plot(p,HFF,g-);grid on title(DFT幅度分析); subplot(1,3,3) plot(p,HF,r-);grid on title(DFT相位分析); 子函数： function y=windows(M,num,beta) M=M+

27、2;%窗函数产生函数,N总长度,M窗长度。 switch num case 1 tmp=boxcar(M);%输入1，产生矩形窗 case 2 tmp=bartlett(M);%输入2 三角窗 case 3 tmp=hanning(M);%汉宁窗 case 4 tmp=hamming(M);%汉明窗 case 5 tmp=blackman(M);%布莱克曼窗 case 6 tmp=kaiser(M,beta);%凯塞窗 otherwise tmp=zeros(1,M); end tmp(1,2,:)=; 19 y=tmp; end 程序二： 凯塞窗 clear all;clc; wp=0.2*

28、2*pi; ws=0.5*2*pi; Rs=60; B=ws-wp; %计算N N=ceil(7.24*pi/B); wc=(wp+ws)/2/pi; h=fir1(N-1,wc,kaiser(N,5.658); figure(1) stem(1:N,h,k-); title(系统频率响应); FF=1024; p=1:FF; temp=fft(h,FF); HFF=20.*log(abs(temp); PFF=angle(temp); figure(2) subplot(2,1,1) plot(p,HFF); title(系统频率响应幅频特性);grid on subplot(2,1,2) plot(p,PFF); title(相位特性);grid on 20