[工学]DSP实验报告

1、实 验 报 告 实验名称 实验一 采样、系统性质及滤波 实验二_系统频率响应和样本处理算法实现 实验三 加窗和离散傅氏变换 实验四 数字滤波器设计 课程名称 数字信号处理a（双语） 班级学号_ _姓 名_ _ 开课时间 实验一 采样、系统性质及滤波一、 实验目的和任务（1）熟悉matlab的主要操作命令。（2）学会简单的矩阵输入和数据读写。（3）掌握简单的绘图命令。（4）用matlab编程并学会创建函数。（5）观察采样引起的混叠。（6）判别离散时间系统的时不变性。（7）卷积计算二、 实验内容与结果a、观察采样引起的混叠设模拟信号为，t的单位为毫秒(ms)。1. 设采样频率为3khz，确定与混叠

2、的采样重建信号。2. 画出和在范围内的连续波形。（因数字计算机无法真正画出连续波形，可用较密的离散点的连线来近似。）3. 分别用和在两信号波形上标记出3khz采样点。两信号波形是否相同？采样后的两序列是否相同？ 实验程序：clear% estimate x(t) and xa(t) with a much higher sampling freq. fs1time_period=6; % unit: msfs1=50; % unit: khzt1=1/fs1; % unit: msn1=0:fix(time_period/t1);x=cos(5*pi*n1*t1)+4*sin(2*pi*n1*

3、t1).*sin(3*pi*n1*t1);xa=cos(pi*n1*t1);% obtain x(nt) and xa(nt) with given sampling freq. fsfs=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),x

4、a(t),x(nt),xa(nt),xlabel(t(ms) 运行结果：b、判别离散时间系统的时不变性。设输入序列为，系统实现对的抽取。1. 设。取延迟量d（例如d30）。记，画出、的序列波形。2. 编程求出系统对的响应以及对的响应3. 画出、的波形。 实验程序：clear% plot x(n) and x(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); endendfigure,subplot(2,1,1),plot(1:n,x,r:,1:length(xd

5、),xd,b),legend(x(n),xd(n),xlabel(n)% plot y(n) and yd(n)for n=1:fix(n/2) y(n)=x(2*n);endfor n=1:length(y)+d, if (n-d)=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); endendsubplot(a(1),a(2),a(3),plot(omega,w1,:),hold on,stem(kk,abs(dft),r.),xlabel(om

6、ega(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); endendsubplot(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=,num2s

7、tr(l), n=,num2str(n)子程序：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 onsubplot(a(4),a(5),a(6),plot(o

8、mega,w2),xlabel(omega(pi),axis(0 2 0 50)title(x(w) with hamming window, l=,num2str(l),grid on子程序：function windows_spectrum(l,a)rect=rectwin(l); % rectangular windowhamm=hamming(l); % hamming windoww=2*pi/1024*(-511:512); % discreted frequencyw1=fft(rect,1024);w1=abs(fftshift(w1);w2=fft(hamm,1024);w2

9、=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)legen

10、d(rectangular,hamming)三、实验分析观察窗长l的变化对窗函数频谱的主瓣宽度、旁瓣密集度、相对旁瓣水平的影响。 答：l越长，主瓣宽度越窄、旁瓣宽度越密集、相对旁瓣水平影响越小。你能否从信号频谱上分辨出信号的三个频率分量？若能分辨出，它们的位置和相对大小是否准确？ 答：能，不准确。矩形窗与汉明窗 答：矩形窗有旁瓣，主瓣宽度是同等条件下汉明窗的二分之一，主瓣高度高于汉明窗。 ，其中，矩形窗c取1，汉明窗c取2。离散频谱dft和连续频谱dtft有什么关系？ 答：dtft定义： l的大小对频率的物理分辨率（physical frequency resolution）有何影响？ l一定

11、的情况下，能否通过增加n改善频率的物理分辨率？n的作用是什么？ 答：l越大，物理分辨率越高。 l一定时，不能通过增加n改善物理分辨率，n的作用是改善计算分辨率。实验四 数字滤波器设计一、实验目的和任务（1）熟悉matlab的主要操作命令。（2）学会简单的矩阵输入和数据读写。（3）掌握简单的绘图命令。（4）用matlab编程并学会创建函数。（5）观察采样引起的混叠。（6）判别离散时间系统的时不变性。（7）卷积计算二、实验内容与结果a、窗口法设计fir数字滤波器1、用矩形窗设计长度分别为n=11、41、81、121的低通fir滤波器，要求截止频率为。画出滤波器的单位冲激响应和幅度频响曲线。 运行结

12、果：2、用汉明窗再次设计同样的滤波器。 运行结果：b、以butterworth 模拟低通滤波器为原型，设计iir数字滤波器。1、截止频率。设计11阶iir数字低通滤波器，画出幅频、相频曲线。2、截止频率。设计11阶iir数字高通滤波器，画出幅频、相频曲线。 运行结果： 实验程序：%=% design fir df with window method% =wc=0.3*pi;% - rectangular windowwin=1; p=1; % observe magnitude responsen=11; figure,fir_design(n,wc,win,2 2 1 2 2 2,p)n=

13、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)p=2; % observe phase responsen=11; figure,fir_design(n,wc,win,1 1 1 1 1 1,p)% - hamming windowwin=2; p=1;n=11; figure,fir_design(n,wc,win,2 2 1 2 2 2,p) n=41; fir_design(n,wc

14、,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)%=%design iir df %=n=11;wn=wc/pi;b,a = butter(n,wn); % lowpassh,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 ii

15、r df: magnitude response),gridsubplot(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); % highpassh,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

16、 iir df: magnitude response),gridsubplot(2,1,2),plot(w/pi,angle(h)/pi), xlabel(omega(pi),ylabel(ang(h(omega)(pi),title(phase response) 子程序：function fir_design(n,wc,win,a,p) % lp fir design with window methodn=0:n-1;m=(n-1)/2;warning off matlab:dividebyzeroh=sin(wc*(n-m)./(pi*(n-m); h(m+1)=wc/pi; if

17、win=2, h =h(:).*hamming(n); endh 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)| ),gridelse subplot(a(1),a(2),a(3),plot(w/pi,

18、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 三、实验分析理想滤波器的频响是怎样的？ 答：lowpass : highpass: bandpass: bansstop: 当n增大时，fir滤波器在c附近的最大纹波幅度是否降低？其余纹波的幅度是否减小？纹波的密度怎样变化？过渡带宽度怎样变化？（如有必要可增大n值观察）。 答：当n增大时，fir滤波片在c附近的最大纹波幅度没有降低，其余纹波的幅度减少，纹波的密度增大，过渡带宽宽度增加。在n=11时，画出滤波器的相频曲线。它是否是线性的？ 答：是线性的用汉明窗设计出的滤波器与用矩形窗相比有什么特点？ 答：用汉明窗设计的滤波器没有纹波，但过渡带宽大。所设计的iir滤波器与fir滤波器的频率特性有何区别？ 答：iir滤波器不是线性相位，但计算量低，时延小。实验总结通过这几次实验，让我增强了使用matlab的能力，同时增强了我对fft、dtft等相关概念的理解，同时对