数字信号处理实验报告

数字信号处理 实验名称:  采样,系统性质及滤波 系统频率响应和样本处理算法实现 加窗和离散傅氏变换 数字滤波器设计 姓名:                      学号:                      专业:                    实验一: 采样,系统性质及滤波 一 、观察采样引起的混叠；熟悉Matlab的使用和基本编程（信号产生，矢量运算，结果显示） 设模拟信号为x(t)=cos(5πt)+4sin(2πt)*sin(3πt),t的单位为毫秒(ms)。 1、设采样频率为3kHz, 通过手工计算确定与x(t) 混叠的采样重建信号 。 2、画出x(t)和 在0 (ms)范围内的连续波形，（因数字计算机无法真正画出连续波形，可用较密的离散点的连线来近似。） 3、分别用o和x在两信号波形上标记出3kHz采样点。 解：(1)经过手工计算得重建新号 =cos(πt) (2)(3)程序代码及波形如下： time_period=6; fs1=50; T1=1/fs1; n1=0:fix(time_period/T1); x=cos(5*pi*n1*T1)+4*sin(2*pi*n1*T1).*sin(3*pi*n1*T1); xa=cos(pi*n1*T1); fs=3; T=1/fs; n=0:fix(time_period/T); x_sample=cos(5*pi*n*T)+4*sin(2*pi*n*T).*sin(3*pi*n*T); xa_sample=cos(pi*n*T); figure,plot(n1*T1,x,'r',n1*T1,xa,'b',n*T,x_sample,'ro'), hold on,stem(n*T,xa_sample,'b:x') legend('x(t)','xa(t)','x(nT)','xa(nT)'),xlabel('t(ms)') 思考编程 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 ： 1.两信号x(t)和xa(t)波形是否相同？采样后的两序列x(nT)和xa(nT) 是否相同？反映了什么现象？ 答：x(t)和 波形不同，采样后两序列波形相同。反映了采样重构时发生混迭现象。 2．参考程序第一段的语句：x=cos(5*pi*n1*T1)+4*sin(2*pi*n1*T1).*sin(3*pi*n1*T1) 中用计算符“*”代替“.*”，结果如何？如果进一步将参与计算的两正弦矢量sin(2*pi*n1*T1)和sin(3*pi*n1*T1)分别进行转置（提示：矢量y的转置为y’），再进行“*”，结果又如何？ 答：前者程序报错。后者计算正确。 4. 改用绿色画出。 答：将程序中的figure,plot(n1*T1,x,'r',n1*T1,xa,'b',n*T,x_sample,'ro'),改为： figure,plot(n1*T1,x,'r',n1*T1,xa,'g',n*T,x_sample,'ro'), 波形如下图 二、判别离散时间系统的时不变性。 设输入序列为x(n),系统y(n)=x(2n)实现对x(n)的抽取。 (1) 设 ,n=1,2,….,500。取延迟量D（例如D=30）。记 ，画出x(n), 的序列波形。 (2) 编程求出系统对x(n)的响应y(n)以及对 的响应 (3) 画出y(n-D), 的波形。 解：程序及波形图如下 D=30; N=500; n=1:N; x=sin(2*pi/100*n); for n=1:N+D, if (n-D)<=0, xD(n)=0; else xD(n)=x(n-D); end end figure,subplot(2,1,1), plot(1:N,x,'r:',1:length(xD),xD,'b') legend('x(n)','xD(n)'),xlabel('n') for n=1:fix(N/2) y(n)=x(2*n); end for n=1:length(y)+D, if (n-D)<=0, y_delay(n)=0; else y_delay(n)=y(n-D); end end for n=1:fix(length(xD)/2) yD(n)=xD(2*n); end subplot(2,1,2), plot(1:length(y),y,'r:',1:length(y_delay),y_delay,'g.:',1:length(yD),yD,'b.') legend('y(n)','y(n-D)','yD(n)'),xlabel('n') axis([0 530 -1 1]); figure,subplot(4,1,1),plot(1:N,x,'r:'),legend('x(n)'),xlabel('n') subplot(4,1,2),plot(1:length(y),y,'k:'),legend('y(n)'),xlabel('n') subplot(4,1,3),plot(1:length(y_delay),y_delay,'g:'),legend('y(n-D)'),xlabel('n') subplot(4,1,4),plot(1:length(yD),yD,'b:'),legend('yD(n)'),xlabel('n') 思考编程题: 1、该系统是否为时不变的？ 答：该系统不是时不变系统。 2、利用subplot命令，自上而下用四个小窗口分别画出x(n), y(n), y(n-D)的波形和yD(n)的波形。 答：利用subplot画出的函数图形如下图所示。 三、利用卷积计算出信号通过FIR滤波器的输出，并观察输出响应的Input-on暂态，input-off暂态和稳态部分。(计算卷积可用conv命令) 考虑下面两个滤波器，第一个的单位脉冲响应为 ，另一个的单位脉冲响应为h= [1,-5,10,-10,5,-1];输入为周期方波，在一个周期内 。分别画出两个滤波器的输出y(n)(0 )的波形，并与书上p144例4.1.8的两幅图比较是否一致。 解：实验程序及波形如下 clear h1=0.25*0.75.^(0:14); h2=1/5*[1 -5 10 -10 5 -1]; N=200; n=0:N-1; x1=[ones(1,25) zeros(1,25)]; % one period of 'x(n)' x=[x1 x1 x1 x1]; y1=conv(x,h1); y2=conv(x,h2); figure,subplot(2,1,1),plot(n,x,'r:',n,y1(1:N),'b'), axis([0 200 -0.5 2.5]),grid on, legend('input','output'),xlabel('n') subplot(2,1,2),plot(n,x,'r:',n,y2(1:N),'b'), axis([0 200 -1.5 2.5]),grid on, legend('input','output'),xlabel('n') 思考编程题: 1．两个滤波器的DC gain分别是多少？ 答：第一个： 第二个： 2.响应的input-on暂态、input-off暂态和稳态段范围分别是多少？ 答：第一个周期：input-on暂态： input-off暂态： 稳态: 实验二:系统频率响应和样本处理算法实现 一、观察序列频谱，观察信号通过系统后波形与频谱的变化 已知输入信号x(n)= (n)+2 ,其中， ， ，n=0,1,2….N-1,N可取5000点。 (1) 画出 ， ，x(n)的前100点波形 (2) 画出x(n)的DTFT频谱X(w)( ) (3) 某LTI系统h(n)=u(n)-u(n-8),画出系统的幅度频响｜H(w)｜ (4) 求系统h(n)对x(n)的响应y(n),画出y(n)的波形并与 的波形比较(各画100点)；画出y(n)的幅度谱|Y(w)|，并与|X(w)H(w)|比较 (1) 解：（1） ， ，x(n)的前100点波形如下图所示 N=5000; x1=cos(pi/8*(0:N-1)); x2=cos(pi/4*(0:N-1)); x=5*x1+2*x2; m=100; figure,plot(1:m,x1(1:m),1:m,x2(1:m),'g',1:m,x(1:m),'r.-'),legend('x1','x2','x') M=5000; w=pi/M*(0:M-1); j=sqrt(-1); Xw=zeros(1,M); for k=1:M, Xw(k)= x*(exp(-j*w(k)*(0:N-1)')); end figure,plot(w/pi,abs(Xw)),xlabel('\omega(\pi)'),ylabel('|X(\omega)|') h=ones(1,8); for k=1:M, Hw(k)=h*(exp(-j*w(k)*(0:7)')); end figure,plot(w/pi,abs(Hw)),xlabel('\omega(\pi)'),ylabel('|H(\omega)|') y=conv(x,h); H1=h*(exp(-j*pi/8*(0:7)')); H2=h*(exp(-j*pi/4*(0:7)')); m=100; yy=5*abs(H1)*cos(pi/8*(0:m-1)+angle(H1))+2*abs(H2)*cos(pi/4*(0:m-1)+angle(H2)); figure,plot(y(1:m),'r.-'),hold on,plot(yy,'g'),legend('y','yy') Ny=length(y); for k=1:M, Yw(k)=y*(exp(-j*w(k)*(0:Ny-1)')); end figure,plot(w/pi,abs(Yw),'r.-') ,hold on,plot(w/pi,abs(Xw.*Hw)), xlabel('\omega(\pi)'),legend('|Y(\omega)|','|X(\omega)|*|H(\omega)|') (1) (2) x(n)的DTFT频谱X(w)( )如下图所示 (3) 系统的幅度频响｜H(w)｜如下图所示 (4)y(n)的波形与幅度谱如下所示 思考编程题: 1.滤波器 是什么频响类型的滤波器? 答： 低通 2.你从以上实验中观察到什么？与课本的什么重要结论相吻合？ （1） （2）  对正弦信号： 二、系统函数 ，根据正准型结构（canonical form）编写样本处理算法。内部状态的初始值设为零，输入信号x(n) 采用逐个样本手动输入的方式(用input命令)，求输出信号y(n)。 解：样本处理算法及输出信号波形如下所示 w1=0; for i=1:10, x=input('input x ='); w0=0.8*w1+x; y=w0-0.5*w1; w1=w0; y end 实验三：加窗和离散傅氏变换 一、观察窗函数的影响 信号为x(t)=cos(2 )+cos(2 )+cos(2 )， =2KHz, =2.5KHz, =3KHz,采样频率 =10KHz。 a) 写出x(n)(- )的频谱 。 b) 分别画出窗长度L=10,L=20,L=40,L=100 的矩形窗频谱和Hamming窗频谱。 c) 时域采样点分别取L=10,L=20,L=40,L=100,分别画出x(n)的加矩形窗及加Hamming窗时DTFT频谱 答：实验程序如下： figure,L=10; windows_spectrum(L,[2 2 1 2 2 2]) L=20; windows_spectrum(L,[2 2 3 2 2 4]) figure,L=40; windows_spectrum(L,[2 2 1 2 2 2]) L=100; windows_spectrum(L,[2 2 3 2 2 4]) L=100; x=cos(0.4*pi*(0:(L-1)))+cos(0.5*pi*(0:(L-1)))+cos(0.6*pi*(0:(L-1))); L=10; figure,windowed_spectrum(x(1:L),[2 2 1 2 2 2]) L=20; windowed_spectrum(x(1:L),[2 2 3 2 2 4]) L=40; figure,windowed_spectrum(x(1:L),[2 2 1 2 2 2]) L=100; windowed_spectrum(x(1:L),[2 2 3 2 2 4]) 子程序1： function windowed_spectrum(x,a) L=length(x); w=hamming(L); W1=abs(fft(x,1024)); W2=abs(fft(x.*w',1024)); omega=(0:1023)*2/1024; subplot(a(1),a(2),a(3)),plot(omega,W1),xlabel('\omega(\pi)'),axis([0 2 0 50]) title(['X(w) with rectangular window, L=',num2str(L)]),grid on subplot(a(4),a(5),a(6)),plot(omega,W2),xlabel('\omega(\pi)'),axis([0 2 0 50]) title(['X(w) with Hamming window, L=',num2str(L)]),grid on 子程序2： function windows_spectrum(L,a) rect=rectwin(L); % rectangular window hamm=hamming(L); % Hamming window w=2*pi/1024*(-511:512); % discreted frequency W1=fft(rect,1024); W1=abs(fftshift(W1)); W2=fft(hamm,1024); W2=abs(fftshift(W2)); subplot(a(1),a(2),a(3)),hold on,plot(1:L,rect,'bo',1:L,hamm,'r*'),stem(1:L,rect),axis([1 L -0.5 2]), xlabel('n'),title(['L=',num2str(L),' (waveform)']),legend('rectangular','Hamming') subplot(a(4),a(5),a(6)),plot(w/pi,abs(W1),w/pi,abs(W2),'r'),xlabel('\omega(\pi)'),title(['L=',num2str(L), ' (Magnitude spectrum)']) legend('rectangular','Hamming') 实验所得图形如下： 思考编程题: 1．观察窗长L的变化对窗函数频谱的主瓣宽度、旁瓣密集度、相对旁瓣水平的影响。 答：L越长，主瓣宽度越窄、旁瓣宽度越密集、相对旁瓣水平影响越小。 2. 你能否从信号频谱上分辨出信号的三个频率分量？若能分辨出，它们的位置和相对大小是否准确？ 答：能，不准确 3.L的大小对频率的物理分辨率（physical frequency resolution）有何影响？ 答：L越大，物理分辨率越高。 二、 理解频率的物理分辨率和计算分辨率的区别 信号同上，加矩形窗。时域采样点数分别取L=10，L=20, L=40, L=100 画出以上各种时长情况下，频域采样点数分别为N=32，N=64时的DFT（在同一个图上用虚线画出相应的DTFT频谱，用于比较）。 答：实验程序如下： L=10; figure,windowed_dft_spectrum(x(1:L),[2 2 1 2 2 2]) L=20; windowed_dft_spectrum(x(1:L),[2 2 3 2 2 4]) L=40; figure,windowed_dft_spectrum(x(1:L),[2 2 1 2 2 2]) L=100; windowed_dft_spectrum(x(1:L),[2 2 3 2 2 4]) 子程序１： function []=windowed_dft_spectrum(x,a); L=length(x); W1=abs(fft(x,1024)); omega=(0:1023)*2/1024; N=32; kk=2/N*(0:N-1); if N>=L, dft=abs(fft(x,N)); else dft=zeros(N,1); for k=1:N, dft(k)=sum(x(:)'*exp(-j*2*pi/N*(k-1)*(0:L-1)')); end end subplot(a(1),a(2),a(3)),plot(omega,W1,':'),hold on,stem(kk,abs(dft),'r.'), xlabel('\omega(\pi)'),axis([0 2 0 50]),title(['L=',num2str(L),', N=',num2str(N)]) N=64; kk=2/N*(0:N-1); if N>=L, dft=abs(fft(x,N)); else dft=zeros(N,1); for k=1:N, dft(k)=sum(x(:)'*exp(-j*2*pi/N*(k-1)*(0:L-1)')); end end subplot(a(4),a(5),a(6)),plot(omega,W1,':'),hold on,stem(kk,abs(dft),'r.') xlabel('\omega(\pi)'),axis([0 2 0 50]),title(['L=',num2str(L),', N=',num2str(N)]) 实验所得图形如下： 思考编程题: 1.离散频谱DFT和连续频谱DTFT有什么关系？ 答： DTFT定义： 2.L一定的情况下，能否通过增加N改善频率的物理分辨率？N的作用是什么？ 答：L一定时，不能通过增加N改善物理分辨率，N的作用是改善计算分辨率。 实验四：数字滤波器设计 一、窗口法设计FIR数字滤波器 (a) 用矩形窗设计长度分别为N=11、41、81、121的低通FIR滤波器，要求截止频率为ωc =0.3π画出滤波器的单位冲激响应和幅度频响h(n)、 | H(ω)|曲线。 (b）用汉明窗再次设计同样的滤波器。 答：实验程序如下： N=121; FIR_design(N,wc,win,[2 2 3 2 2 4],p) p=2; N=11; figure,FIR_design(N,wc,win,[1 1 1 1 1 1],p) win=2; p=1; N=11; figure,FIR_design(N,wc,win,[2 2 1 2 2 2],p) N=41; FIR_design(N,wc,win,[2 2 3 2 2 4],p) N=81; figure,FIR_design(N,wc,win,[2 2 1 2 2 2],p) N=121; FIR_design(N,wc,win,[2 2 3 2 2 4],p) 子程序： function FIR_design(N,wc,win,a,p) % LP FIR design with window method n=0:N-1; M=(N-1)/2; warning off MATLAB:divideByZero h=sin(wc*(n-M))./(pi*(n-M)); h(M+1)=wc/pi; if win==2, h =h(:).*hamming(N); end [H w]=freqz(h,1,16400); if p==1, subplot(a(1),a(2),a(3)),stem(0:N-1,h,'.'),xlabel('n'),ylabel('h(n)'),title(['N=',num2str(N),', h(n)']) subplot(a(4),a(5),a(6)),plot(w/pi,abs(H)),xlabel('\omega(\pi)'),ylabel('|H(\omega)|') ,title(['N=',num2str(N),', |H(\omega)| ']),grid else subplot(a(1),a(2),a(3)),plot(w/pi,angle(H)/pi,w/pi,abs(H),'r:'),legend('ang(H(\omega)','|H(\omega)|') xlabel('\omega(\pi)'),ylabel('ang(H(\omega))(\pi)'),title(['N=',num2str(N),', phase response']) end 实验所得图形如下： 思考编程题: 1.理想滤波器的频响是怎样的？ 答： Lowpass : Highpass: Bandpass: Bansstop: 2.当N增大时，FIR滤波器在ωc附近的最大纹波幅度是否降低？其余纹波的幅度是否减小？纹波的密度怎样变化？过渡带宽度怎样变化？（如有必要可增大N值观察）。 答：当N增大时，FIR滤波片在ωc附近的最大纹波幅度没有降低，其余纹波的幅度减少，纹波的密度增大，过渡带宽宽度减小。 3.在N=11时，画出滤波器的相频曲线。它是否是线性的？ 答：是线性的 4.用汉明窗设计出的滤波器与用矩形窗相比有什么特点？ 答：用汉明窗设计的滤波器没有纹波，但过渡带宽大。 二、以Butterworth 模拟低通滤波器为原型，设计IIR数字滤波器。 （a）截止频率ωc =0.3π。设计11阶IIR数字低通滤波器，画出幅频、相频曲线。 （b）截止频率ωc =0.3π。设计11阶IIR数字高通滤波器，画出幅频、相频曲线。 答：实验程序如下： N=11; Wn=wc/pi; [b,a] = BUTTER(N,Wn); % lowpass [H,w]=freqz(b,a,8192); figure,subplot(2,1,1),plot(w/pi,abs(H)),xlabel('\omega(\pi)'),ylabel('|H(\omega)|') ,title(['order=',num2str(N),', Lowpass IIR DF: magnitude response']),grid subplot(2,1,2),plot(w/pi,angle(H)/pi), xlabel('\omega(\pi)'),ylabel('ang(H(\omega))(\pi)'),title('phase response') N=11; [b,a] = BUTTER(N,Wn,'high'); % highpass [H,w]=freqz(b,a,8192); figure,subplot(2,1,1),plot(w/pi,abs(H)),xlabel('\omega(\pi)'),ylabel('|H(\omega)|') ,title(['order=',num2str(N),', HIghpass IIR DF: magnitude response']),grid subplot(2,1,2),plot(w/pi,angle(H)/pi), xlabel('\omega(\pi)'),ylabel('ang(H(\omega))(\pi)'),title('phase response') 实验所得图形如下： 思考编程题: 1.所设计的IIR滤波器与FIR滤波器的频率特性有何区别？ 答：IIR滤波器不是线性相位，但计算量低，时延小。