现代信号处理课程设计报告

现代信号处理课程设计报告 目录一、 课程设计目的2二、 课程设计内容2三、 课程设计流程及分析43.1熟悉MATLAB基本操作. 43.2 MATLAB数字信号处理图形界面. 53.3用MATLAB绘出信号波形. .5 3.4对信号频谱分析. . . . . . . . .73.5设计滤波器并对信号滤波.11四、 设计经验总结27五、 参考资料27一、课程设计目的1全面复习课程所学理论知识，巩固所学知识重点和难点，将理论与实践很好地结合起来。2提高综合运用所学知识独立分析和解决问题的能力；3熟练使用一种高级语言进行编程实现。二、课程设计内容（一）信号分析1、编制信号生成程序，产生下述各序列，绘出它们的时域波形 1） 单位抽样序列 2）矩形序列 3） 三角波序列4) 反三角波序列 5）Gaussian序列6） 正弦序列 取7）衰减正弦序列 对连续信号 进行采样，可得到测试序列 。令A=50,采样周期T=1ms,即fs=1000Hz,f0=62.5,a=100。2. 对上述信号完成下列信号分析 1）对三角波序列和反三角波序列，作N=8点的FFT，观察比较它们的幅频特性，说明它们有什么异同？绘出两序列及其它们的幅频特性曲线。 在和的尾部补零，作N=16点的FFT，观察它们的幅频特性发生了什么变化？分析说明原因。 2)、观察高斯序列，固定信号中的参数p=8，令q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，令p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。3）对于正弦序列x6(n)，取数据长度N分别等于8，16，32，分别作N点FFT，观察它们的的时域和幅频特性，说明它们的差别 ，简要说明原因。4)、观察衰减正弦序列的时域和幅频特性，绘出幅频特性曲线，改变采样频率fs,使fs=300Hz， 观察此时的频谱的形状和谱峰出现位置？说明产生现象的原因。3设有一连续时间信号s(t),其由20Hz、220Hz和750Hz的正弦信号叠加而成，分析确定采样频率及数据分析长度，计算并绘出信号的频谱，指出各个频率份量。 （二）、数字滤波器的设计与实现 设计和实现IIR或FIR数字滤波器： 要求：输入数字滤波器的滤波指标,包括通带截止频率,通带最大衰减,阻带截止频率,阻带最小衰减。设计出相应的数字滤波器。显示得到的滤波器的系数和阶数，求出系统的单位脉冲响应h(n)。并画出设计得到的滤波器的幅频特性图(要有坐标标度)和相频特性。 1对于IIR数字滤波器：使用双线性变换法。可选用BW型、CB型和椭圆型等滤波器类型来进行设计2对于FIR数字滤波器：要对多种窗口（三种及以上）、不同阶数设计结果进行分析比较。（三）、综合应用对于给出的任意一段语音信号，进行分析处理。1、 读出语音信号，并显示语音信号的时域波形，然后对读出的信号进行FFT变换，得到信号的频谱特性。2、 分析其频谱特性，应用设计的滤波器将主要的频率分量分离出来。3、 回放此语音信号，感觉滤波前后的声音有什么变化三、课程设计流程3.1熟悉MATLAB基本操作，基本语法，以及GUI界面设计中各控件的调用方式。MATLAB是一款功能超强的数学软件，应用于各个行业。而其基本操作却很大众化，操作起来很人性化。其中所附带的GUI面向对象的友好编程方式更便于学术交流和设计开发。 3.2 MATLAB数字信号处理图形界面3.3用MATLAB绘出信号波形global xypopup_sel_index = get(handles.popupmenu2, Value);switch popup_sel_index case 1 n0=0;n1=-10;n2=10;n=n1:n2;x=(n-n0)=0;stem(n,x,.);xlabel(n);ylabel(x(n);title(单位抽样序列); case 2 n=-1:10;n0=2;width=5;x=rectpuls(n-n0,width);stem(n,x,.);xlabel(n);ylabel(x(n);axis(-1 10 0 2);title(矩形序列); case 3 n1=0:3;n2=4:7;x1=n1+1;x2=8-n2;n=(-10:0) n1 n2 (7:17);x=zeros(1,11) x1 x2 zeros(1,11);stem(n,x,.);xlabel(n);ylabel(x(n);title(三角波序列); case 4 n1=0:3;n2=4:7;x1=4-n1;x2=n2-3;n=(-10:0) n1 n2 (7:17);x=zeros(1,11) x1 x2 zeros(1,11);stem(n,x,.);xlabel(n);ylabel(x(n);title(反三角波序列); case 5 n=0:15;p=8;q=2;x=exp(-(n-p).2/q);stem(n,x,.);xlabel(n);ylabel(x(n);title(Gaussian序列); case 6 t=1/64;n=0:t:16;phase=0;A=1;arg=2*pi*n-phase;x=A*sin(arg);stem(n,x,.);xlabel(n);ylabel(x(n);title(正弦序列); case 7 t=0.001;n=0:t:0.2;x1=sin(2*62.5*n*pi);x2=exp(-100*n);x=x1.*x2;stem(n,x,.);xlabel(n);ylabel(x(n);title(衰减正弦序列); case 8f1=20;f2=220;f3=750;n=0:0.5/f3:1;x=sin(2*pi*f1*n)+sin(2*pi*f2*n)+sin(2*pi*f3*n);stem(n,x,.);xlabel(n);ylabel(x(n);title(正弦叠加序列); endxy=x;3.4对信号频谱分析三角波的8点傅里叶变换： 三角波的32点傅里叶变换：对信号N点离散傅里叶变换相当于序列傅里叶变换在频域里的N点采样,所以采样点数越高越能反映信号的频谱。Gaussian序列P=8，Gaussian序列的15点傅里叶变换P=8，Gaussian序列的15点傅里叶变换由图可知，如果N=15，当P=13时，会发生明显的泄漏现象，混叠现象也随之出现。当采样频率fs=300时，衰减正弦序列的频谱图当采样频率fs=3000时，衰减正弦序列的频谱图采样频率为fs=300时，由于采样频率太小，产生频谱混叠显现。当采样频率为1500时,正弦叠加序列的频谱图当采样频率为1500时,正弦叠加序列的频谱图由于正弦叠加序列的频率由2*20、2*220、 2* 750三种频率组成，所以采样频率应大于或等于1500；由于三个正弦序列的周期都为1，所以截取0到1进行傅里叶分析。3.5设计滤波器并对信号滤波 global yglobal fsviir=get(handles.IIR,Value);vfir=get(handles.FIR,Value);vbw=get(handles.BW,Value);vcb=get(handles.CB,Value);vclass1=get(handles.class1,Value);vclass2=get(handles.class2,Value);vclass3=get(handles.class3,Value);vclass4=get(handles.class4,Value);vpopupmenu1=get(handles.popupmenu1, Value);vp1=str2num(get(handles.p1,String);vs1=str2num(get(handles.s1,String);vp2=str2num(get(handles.p2,String);vs2=str2num(get(handles.s2,String);vap=str2num(get(handles.ap,String);vas=str2num(get(handles.as,String);fs1=str2num(get(handles.Fs,String);vs=get(handles.checkbox1,Value);if(viir) if(vs) T=1/fs; else T=1/fs1; end if(vbw) if(vclass1) wp=2/T*tan(vp1*pi/2); ws=2/T*tan(vs1*pi/2); N,wc=buttord(wp,ws,vap,vas,s); b,a=butter(N,wc,s); elseif(vclass2) wp=2/T*tan(vp1*pi/2); ws=2/T*tan(vs1*pi/2); N,wc=buttord(wp,ws,vap,vas,s); b,a=butter(N,wc,high,s); elseif(vclass3) wp=2/T*tan(vp1*pi/2),2/T*tan(vp2*pi/2); ws=2/T*tan(vs1*pi/2),2/T*tan(vs2*pi/2); N,wc=buttord(wp,ws,vap,vas,s); b,a=butter(N,wc,s); else wp=2/T*tan(vp1*pi/2),2/T*tan(vp2*pi/2); ws=2/T*tan(vs1*pi/2),2/T*tan(vs2*pi/2); N,wc=buttord(wp,ws,vap,vas,s); b,a=butter(N,wc,stop,s); end else if(vclass1) wp=2/T*tan(vp1*pi/2); ws=2/T*tan(vs1*pi/2); N,wc=cheb1ord(wp,ws,vap,vas,s); b,a=cheby1(N,vap,wc,s); elseif(vclass2) wp=2/T*tan(vp1*pi/2); ws=2/T*tan(vs1*pi/2); N,wc=cheb1ord(wp,ws,vap,vas,s); b,a=cheby1(N,vap,wc,high,s); elseif(vclass3) wp=2/T*tan(vp1*pi/2),2/T*tan(vp2*pi/2); ws=2/T*tan(vs1*pi/2),2/T*tan(vs2*pi/2); N,wc=cheb1ord(wp,ws,vap,vas,s); b,a=cheby1(N,vap,wc,s); else wp=2/T*tan(vp1*pi/2),2/T*tan(vp2*pi/2); ws=2/T*tan(vs1*pi/2),2/T*tan(vs2*pi/2); N,wc=cheb1ord(wp,ws,vap,vas,s); b,a=cheby1(N,vap,wc,stop,s); end end if(vs) bz,az=bilinear(b,a,fs); else bz,az=bilinear(b,a,fs1); end if(vs) x=filter(bz,az,y); X=fft(x,4096); plot(0:length(x)-1,x);title(滤波后的信号波形); figure; f=(0:4095)/4096; plot(f*2,abs(X); xlabel(w/);ylabel(X|k|); title(滤波后的频谱); sound(x,fs); else wk=0:pi/512:pi; hz,w=freqz(bz,az,wk); hz1=ifft(hz); stem(w,hz1,.);figure; L=length(hz); hx=angle(hz); plot(2*(0:L-1),hx,.); xlabel(w/);ylabel(相位);figure; plot(wk/pi,abs(hz); xlabel(w/);ylabel(幅度X|k|);figure; plot(0:L-1)/L,20*log10(abs(hz);grid on; axis(0 1 -100 5); xlabel(w/);ylabel(幅度(db); end N wc b a else if(vclass1|vclass2) DB=abs(vs1-vp1); wc=(vp1+vs1)/2; switch vpopupmenu1 case 1 if(vasvp1&vclass1) N=N0; hn=fir1(N,wc,boxcar(N+1); elseif(vs1vp1&vclass2) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,high,boxcar(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(矩形窗的阻带最小衰减小于或等于21,Error); end case 2 if(vasvp1&vclass1) N=N0; hn=fir1(N,wc,bartlett(N+1); elseif(vs1vp1&vclass2) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,high,bartlett(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(三角窗的阻带最小衰减小于或等于25,Error); end case 3 if(vasvp1&vclass1) N=N0; hn=fir1(N,wc,hanning(N+1); elseif(vs1vp1&vclass2) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,high,hanning(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(汉宁窗的阻带最小衰减小于或等于44,Error); end case 4 if(vasvp1&vclass1) N=N0; hn=fir1(N,wc,hamming(N+1); elseif(vs1vp1&vclass2) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,high,hamming(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(哈明窗的阻带最小衰减小于或等于53,Error); end case 5 if(vasvp1&vclass1) N=N0; hn=fir1(N,wc,blackman(N+1); elseif(vs1vp1&vclass2) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,high,blackman(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(布莱克曼窗的阻带最小衰减小于或等于74,Error); end case 6 if(vas50) beta=0.112*(vas-8.7); elseif(vas=21) beta=0.5842*(vas-21)0.4+0.07886*(vas-21); else beta=0; end N0=ceil(vas-8)/2.285/DB); if(vs1vp1&vclass1) N=N0; hn=fir1(N,wc,kaiser(N+1,beta); elseif(vs1vp1&vclass2) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,high,kaiser(N+1,beta); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(凯塞窗的阻带最小衰减小于或等于80,Error); end end elseif(vclass3|vclass4) DB=abs(vs1-vp1); wc=(vs1+vp1)/2,(vs2+vp2)/2; switch vpopupmenu1 case 1 if(vas=21) N0=ceil(1.8/DB); if(vs1vp1&vclass4) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,stop,boxcar(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(矩形窗的阻带最小衰减小于或等于21,Error); end case 2 if(vas=25) N0=ceil(6.1/DB); if(vs1vp1&vclass4) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,stop,bartlett(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(三角窗的阻带最小衰减小于或等于25,Error); end case 3 if(vas=44) N0=ceil(6.2/DB); if(vs1vp1&vclass4) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,stop,hanning(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(汉宁窗的阻带最小衰减小于或等于44,Error); end case 4 if(vas=53) N0=ceil(6.6/DB); if(vs1vp1&vclass4) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,stop,hamming(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(哈明窗的阻带最小衰减小于或等于53,Error); end case 5 if(vas=74) N0=ceil(11/DB); if(vs1vp1&vclass4) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,stop,blackman(N+1); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(布莱克曼窗的阻带最小衰减小于或等于74,Error); end case 6 if(vas50) beta=0.112*(vas-8.7); elseif(vas=21) beta=0.5842*(vas-21)0.4+0.07886*(vas-21); else beta=0; end N0=ceil(vas-8)/2.285/DB); if(vs1vp1&vclass4) if(mod(N0+1,2) N=N0; else N=N0+1; end hn=fir1(N,wc,stop,kaiser(N+1,beta); else errordlg(输入错误，将通带截止频率与阻带截止频率互换,Error); end else errordlg(凯塞窗的阻带最小衰减小于或等于80,Error); end end end if(vs) x=fftfilt(hn,y); X=fft(x,4096); plot(0:length(x)-1,x); xlabel(n);ylabel(h(n);figure; f=(0:4095)/4096; plot(f*2,abs(X); xlabel(w/);ylabel(X|k|); title(滤波后的频谱); sound(x,fs); else stem(0:N,hn,.); xlabel(n);ylabel(h(n);figure; y=fft(hn,4096*2); M=1:4096; mag=abs(y(M); w=(0:4095)/4096; plot(w,mag);grid on; xlabel(w/);ylabel(幅度X|k|);figure; plot(w,20*log10(mag);grid on; xlabel(w/);ylabel(幅度db);figure; phay=angle(y(M); plot(2*(0:length(phay)-1)/length(phay),phay); xlabel(w/);ylabel(相位(w); end N end用双线性变换法设计IIR滤波器和用窗函数法设计FIR滤波器低通IIR滤波器（通带截止频率0.2，阻带截止频率0.3，通带最大衰减系数1，阻带最小衰减系数40）序列图：频谱图：损耗函数：N = 6wc = 2.5994e+004b = 1.0e+025 * Columns 1 through 3 0 0 0 Columns 4 through 6 0 0 0 Column 71.8943a =1.0e+025 * Columns 1 through 3 0.0000 0.0000 0.0000 Columns 4 through 6 0.0000 0.0000 0.0004 Column 72.1255带通IIR滤波器（通带截止频率0.30.4，阻带截止频率0.20.5，通带最大衰减系数1，阻带最小衰减系数40）序列图：频谱图：损耗函数：N = 4wc = 1.0e+004 * 4.0762 5.8123b = 1.0e+016 * Columns 1 through 7 0 0 0 0 2.2318 0 0 Columns 8 through 9 0 0a = 1.0e+037 * Columns 1 through 7 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 8 through 9 0.0000 3.1508带通FIR滤波器（通带截止频率0.30.4，阻带截止频率0.20.5，通带最大衰减系数1，阻带最小衰减系数40）（下一页）序列图：频谱图：耗函数：N = 67对语音信号频谱分析并滤波对已有函数频谱分析并滤波四、设计经验总结Malta无疑是一款功能超强的数学软件，基于它可实现多领域的科学研究和开发。在本次课程设计中,收获最多的是如何看待问题以及如何解决问题。一个很现实的问题就是搞科研开发，遇到问题是必不可免，关键是如何看待这样的问题，持怎样的心态，持怎样的目的。再就是问题是遇到了，要如何解决以及用怎样最优的方案去处理，都得下心思去思考和摸索。再一个很现实的问题是，要学会学习，学会自学。毕竟社会日新月异，科技蓬勃发展，新领域新知识不断涌现，怎么样才能在社会保持有效的立足，怎样才能翻弄潮流，面对社会变迁做到运筹帷幄之中呢？在我看来，关键要保持良好心态，保持一种积极永不妥协永不放弃的信念，不断进取不断开拓。在学好基本课内的专业知识之外，应尽量多补充课外的和本专业相协调的专业知识。强势未来就业利剑。五、参考资料MATLAB数字信号处理与运用数字信号处理28