《现代信号处理》课程设计

1、目 录一、 课程设计目的-2二、 课程设计任务-2三、 详细设计过程-4四、 设计结果及分析-13五、 心得体会-24一、课程设计目的：1全面复习课程所学理论知识，巩固所学知识重点和难点，将理论与实践很好地结合起来。2提高综合运用所学知识独立分析和解决问题的能力；3熟练使用一种高级语言进行编程实现。二、课程设计内容（一）信号分析1、编制信号生成程序，产生下述各序列，绘出它们的时域波形 1） 单位抽样序列 2）矩形序列 3） 三角波序列4) 反三角波序列 5）gaussian序列6） 正弦序列 取7）衰减正弦序列 对连续信号 进行采样，可得到测试序列 。令a=50,采样周期t=1ms,即fs=1

2、000hz,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、的现象，绘出相应的时域序列和幅频特性曲线。3）对于正弦序列，取数据长度n分别等于8，16，32，分别作n点fft，观察它们的的时域和幅频特性，说明它们的差别 ，简要说明原因。4)、观察衰减正弦序列的时域和幅频特性，绘出幅频特性曲线，改变采样频率fs,使fs=300hz， 观察此时的频谱的形状和谱峰出现位置？说明产生现象的原因。3设有一连续时间信号s(t),其由20hz、220hz和750hz的正弦信号叠加而成，分析确定采样频率及数据分析长度，计算并绘出信号的频谱，指出各个频率份量。 （二）、数字滤波器的设计与实现 设计和实现iir或fir数字滤波器： 要求：输入数字滤波器的滤波指标,包括通带截

4、止频率,通带最大衰减,阻带截止频率,阻带最小衰减。设计出相应的数字滤波器。显示得到的滤波器的系数和阶数，求出系统的单位脉冲响应h(n)。并画出设计得到的滤波器的幅频特性图(要有坐标标度)和相频特性。 1对于iir数字滤波器：使用双线性变换法。可选用bw型、cb型和椭圆型等滤波器类型来进行设计2对于fir数字滤波器：要对多种窗口（三种及以上）、不同阶数设计结果进行分析比较。（三）、综合应用对于给出的任意一段语音信号，进行分析处理。1、 读出语音信号，并显示语音信号的时域波形，然后对读出的信号进行fft变换，得到信号的频谱特性。2、 分析其频谱特性，应用设计的滤波器将主要的频率分量分离出来。3、

5、回放此语音信号，感觉滤波前后的声音有什么变化三、详细设计过程1、新建一个空白的图形gui界面2、添加axes及pop-up menu。可以用鼠标拖拽，将其放在合适的位置。pop-up menu下拉菜单可以根据任务自行设置名称。3、在图形界面生成的源程序中编程，以实现课程设计的各项任务。源程序如下：function varargout = wo(varargin)% wo m-file for wo.fig% wo, by itself, creates a new wo or raises the existing% singleton*.% h = wo returns the handle

6、 to a new wo or the handle to% the existing singleton*.% wo(callback,hobject,eventdata,handles,.) calls the local% function named callback in wo.m with the given input arguments.% wo(property,value,.) creates a new wo or raises the% existing singleton*. starting from the left, property value pairs a

7、re% applied to the gui before wo_openingfunction gets called. an% unrecognized property name or invalid value makes property application% stop. all inputs are passed to wo_openingfcn via varargin.% *see gui options on guides tools menu. choose gui allows only one% instance to run (singleton).% see a

8、lso: guide, guidata, guihandles % copyright 2002-2003 the mathworks, inc. % edit the above text to modify the response to help wo % last modified by guide v2.5 15-sep-2010 08:37:54 % begin initialization code - do not editgui_singleton = 1;gui_state = struct(gui_name, mfilename, . gui_singleton, gui

9、_singleton, . gui_openingfcn, wo_openingfcn, . gui_outputfcn, wo_outputfcn, . gui_layoutfcn, , . gui_callback, );if nargin & ischar(varargin1) gui_state.gui_callback = str2func(varargin1);end if nargout varargout1:nargout = gui_mainfcn(gui_state, varargin:);else gui_mainfcn(gui_state, varargin:);end

10、% end initialization code - do not edit % - executes just before wo is made visible.function wo_openingfcn(hobject, eventdata, handles, varargin)% this function has no output args, see outputfcn.% hobject handle to figure% eventdata reserved - to be defined in a future version of matlab% handles str

11、ucture with handles and user data (see guidata)% varargin command line arguments to wo (see varargin) % choose default command line output for wohandles.output = hobject; % update handles structureguidata(hobject, handles); % uiwait makes wo wait for user response (see uiresume)% uiwait(handles.figu

12、re1); % - outputs from this function are returned to the command line.function varargout = wo_outputfcn(hobject, eventdata, handles) % varargout cell array for returning output args (see varargout);% hobject handle to figure% eventdata reserved - to be defined in a future version of matlab% handles

13、structure with handles and user data (see guidata) % get default command line output from handles structurevarargout1 = handles.output; % - executes on selection change in popupmenu1.function popupmenu1_callback(hobject, eventdata, handles)% hobject handle to popupmenu1 (see gcbo)% eventdata reserve

14、d - to be defined in a future version of matlab% handles structure with handles and user data (see guidata) n = get(handles.popupmenu1, value);switch n case 1 n=1:50;x=zeros(1,50);x(1)=1;subplot(421); stem(x,.);xlabel(n);ylabel(x1(n); n=1:50;x=rectpuls(n-25,20);subplot(422);stem(x,.);xlabel(n);ylabe

15、l(x2(n); n1=0:3;x1=n1+1; n2=4:7;x2=8-n2;n=n1 n2; x=x1 x2;subplot(423);stem(x,.);xlabel(n);ylabel(x3(n); n1=0:3;x1=4-n1;n2=4:7;x2=n2-3;n=n1 n2;x=x1 x2;subplot(424);stem(x,.);xlabel(n);ylabel(x4(n); n=0:50;p=8;q=2;x=exp(-1*(n-p).2/q);subplot(425);stem(x,.);xlabel(n);ylabel(x5(n); fs=64;n=16;n=0:2*pi*f

16、s/n;x=sin(n*pi/4);subplot(426);stem(x,.);xlabel(n);ylabel(x6(n); n=0:15; a=0.1; f=0.0625; x=50*exp(-a*n).*sin(2*pi*f*n);subplot(427);stem(x,.);xlabel(n);ylabel(x7(n); case 2 n1=0:3;x1=n1+1; n2=4:7;x2=8-n2;n=n1 n2; x=x1 x2; y=fft(x,8);subplot(421);stem(y,.);xlabel(k);ylabel(x3(k);grid; magy=abs(y);k=

17、0:7;w=2*pi/8*k;subplot(422); stem(w/pi,magy,.);xlabel(n=8);ylabel(|x3(k)|);grid; x=x zeros(1,8);y=fft(x,16);subplot(425);stem(y,.);xlabel(k);ylabel(x3(k);grid; magy=abs(y);k=0:15;w=2*pi/16*k;subplot(426); stem(w/pi,magy,.);xlabel(n=16);ylabel(|x3(k)|);grid; n1=0:3;x1=4-n1;n2=4:7;x2=n2-3;n=n1 n2;x=x1

18、 x2; y=fft(x,8);subplot(423);stem(y,.);xlabel(k);ylabel(x4(k);grid; magy=abs(y);k=0:7;w=2*pi/8*k;subplot(424); stem(w/pi,magy,.);xlabel(n=8);ylabel(|x4(k)|);grid; x=x zeros(1,8);y=fft(x,16);subplot(427);stem(y,.);xlabel(k);ylabel(x4(k);grid; magy=abs(y);k=0:15;w=2*pi/16*k;subplot(428); stem(w/pi,mag

19、y,.);xlabel(n=16);ylabel(|x4(k)|);grid; case 3 figure(1);n=0:15;p=8;q=2;x=exp(-1*(n-p).2/q);subplot(321);stem(x,.);title(高斯，p=8,q=2); y=fft(x);magy=abs(y);subplot(322);stem(magy,.);title(幅频特性，p=8,q=2); n=0:15;p=8;q=4;x=exp(-1*(n-p).2/q);subplot(323);stem(x,.);title(高斯，p=8,q=4); y=fft(x);magy=abs(y);

20、subplot(324);stem(magy,.);title(幅频特性，p=8,q=2); n=0:15;p=8;q=8;x=exp(-1*(n-p).2/q);subplot(325);stem(x,.);title(高斯，p=8,q=8); y=fft(x);magy=abs(y);subplot(326);stem(magy,.);title(幅频特性，p=8,q=2); figure(2);n=0:15;p=8;q=8;x=exp(-1*(n-p).2/q);subplot(321);stem(x,.);title(高斯，p=8,q=8); y=fft(x);magy=abs(y);

21、subplot(322);stem(magy,.);title(幅频特性，p=8,q=8); n=0:15;p=13;q=8;x=exp(-1*(n-p).2/q);subplot(323);stem(x,.);title(高斯，p=13,q=8); y=fft(x);magy=abs(y);subplot(324);stem(magy,.);title(幅频特性，p=13,q=8); n=0:15;p=14;q=8;x=exp(-1*(n-p).2/q);subplot(325);stem(x,.);title(高斯，p=14,q=8); y=fft(x);magy=abs(y);subpl

22、ot(326);stem(magy,.);title(幅频特性，p=14,q=8); case 4 fs=64;n=8;n=0:2*pi*fs/n;x=sin(n*pi/4);subplot(321);stem(x,.);title(正弦，n=8);y=abs(fft(x,8);subplot(322);stem(y,.);title(正弦幅频特性，n=8); fs=64;n=16;n=0:2*pi*fs/n;x=sin(n*pi/4);subplot(323);stem(x,.);title(正弦，n=16);y=abs(fft(x,16);subplot(324);stem(y,.);ti

23、tle(正弦幅频特性，n=16); fs=64;n=20;n=0:2*pi*fs/n;x=sin(n*pi/4);subplot(325);stem(x,.);title(正弦，n=20);y=abs(fft(x,20);subplot(326);stem(y,.);title(正弦幅频特性，n=20); case 5 n=0:15; a=0.1; f=0.0625; x=50*exp(-a*n).*sin(2*pi*f*n);subplot(221);stem(x,.);title(fs=1000hz); y=abs(fft(x);subplot(222);stem(y,.);title(f

24、s=1000hz,幅频特性); n=0:15; a=1/3; f=5/24; x=50*exp(-a*n).*sin(2*pi*f*n);subplot(223);stem(x,.);title(fs=300hz); y=abs(fft(x);subplot(224);stem(y,.);title(fs=300hz，幅频特性); case 6 fs=2000;n=200;n=0:n-1;t=n/fs;x=sin(pi*40*t)+sin(pi*440*t)+sin(pi*1500*t);subplot(2,1,1);plot(x); y=fft(x,n);mag=abs(y);f=(0:le

25、ngth(y)-1).*fs/length(y);subplot(212);plot(f(1:n/2),mag(1:n/2);xlabel(f/hz); case 7 wp=input(please enter wp:); ws=input(please enter ws:); rp=input(please enter rp:); rs=input(please enter rs:); k=input(请选择滤波器类型 1、cb型 2、bw型 3、椭圆型:); switch k case 1 %cb型 n,wn= cheb1ord(wp,ws,rp,rs) %阶数n和截止频率wn b,a=c

26、heby1(n,rp,wn) %系统函数系数向量b和a figure(1); freqs(b,a);h,w = freqz(b,a); magh = abs(h);phah = angle(h); num=b; den=a; figure(2);subplot(311);plot(w/pi,magh);grid xlabel(frequency unit:pi);ylabel(magnitude); title(幅频特性); subplot(312);plot(w/pi,phah/pi);grid xlabel(frequency unit:pi);ylabel(phase unit:pi);

27、 title(相频特性); subplot(313); impulse(num,den); title(系统的单位脉冲响应h(n); case 2 %bw型 n,wn=buttord(wp,ws,rp,rs) %阶数n和截止频率wn b,a=butter(n,wn) %系统函数系数向量b和a figure(1); freqs(b,a);h,w = freqz(b,a); magh = abs(h);phah = angle(h); num=b; den=a; figure(2);subplot(311);plot(w/pi,magh);grid xlabel(frequency unit:pi

28、);ylabel(magnitude); title(幅频特性); subplot(3,1,2);plot(w/pi,phah/pi);grid xlabel(frequency unit:pi);ylabel(phase unit:pi); title(相频特性); subplot(313); impulse(num,den); title(系统的单位脉冲响应h(n); case 3 %椭圆型 n,wn=ellipord(wp,ws,rp,rs) %阶数n和截止频率wn b,a=ellip(n,rp,rs,wn) %系统函数系数向量b和a figure(1); freqs(b,a); h,w

29、 = freqz(b,a); magh = abs(h);phah = angle(h); num=b; den=a; figure(2);subplot(311);plot(w/pi,magh);grid xlabel(frequency unit:pi);ylabel(magnitude); title(幅频特性); subplot(312);plot(w/pi,phah/pi);grid xlabel(frequency unit:pi);ylabel(phase unit:pi); title(相频特性); subplot(313); impulse(num,den); title(系

30、统的单位脉冲响应h(n); end case 8 n=21;wc=1/4;n=0:20; hn=fir1(n-1,wc,boxcar(n); fh=fft(hn,1024); fh=20*log10(abs(fh); wk=0:1023;wk=2*wk/1024; subplot(421);stem(n,hn);grid; title(矩形窗设计结果：); subplot(422);plot(wk,fh);grid; title(矩形窗);xlabel(w/pi); wind=hanning(n); hn=hn.*wind; fh=fft(hn,1024); fh=20*log10(abs(f

31、h); subplot(423);stem(n,hn);grid; title(汉宁窗设计结果：); subplot(424);plot(wk,fh);grid; title(汉宁窗);xlabel(w/pi); axis(0 2 -130 0); wind=hamming(n); hn=hn.*wind; fh=fft(hn,1024); fh=20*log10(abs(fh); subplot(425);stem(n,hn);grid; title(哈明窗设计结果：); subplot(426);plot(wk,fh);grid; title(哈明窗);xlabel(w/pi); axis

32、(0 2 -160 0); wind=blackman(n); hn=hn.*wind; fh=fft(hn,1024); fh=20*log10(abs(fh); subplot(427);stem(n,hn);grid; title(布莱克曼窗设计结果：); subplot(428);plot(wk,fh);grid; title(布莱克曼窗);xlabel(w/pi); case 9 fb=1000;fc=1200;as=100;ap=1;fs=22050;x2=wavread(f:sha we.wav );ts=1/fs;r1=10;wp=2*pi*fb/fs;ws=2*pi*fc/f

33、s;wp1=2/ts*tan(wp/2); %将模拟指标转换成数字指标ws1=2/ts*tan(ws/2);n,wn=buttord(wp1,ws1,ap,r1,s); %选择滤波器的最小阶数z,p,k=buttap(n); %创建butterworth模拟滤波器bap,aap=zp2tf(z,p,k);b,a=lp2lp(bap,aap,wn); bz,az=bilinear(b,a,fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换h,w=freqz(bz,az); %绘制频率响应曲线figure(1)plot(w*fs/(2*pi),abs(h);gridxlabel(频率hz)

34、ylabel(频率响应幅度)title(iir低通滤波器)f1=filter(bz,az,x2);figure(2)subplot(2,1,1)plot(x2) %画出滤波前的时域图title(iir低通滤波器滤波前的时域波形);subplot(2,1,2)plot(f1); %画出滤波后的时域图title(iir低通滤波器滤波后的时域波形);sound(f1,44100); %播放滤波后的信号f0=fft(f1,1024);f=fs*(0:511)/1024;figure(3)y2=fft(x2,1024);subplot(2,1,1);plot(f,abs(y2(1:512); grid

35、%画出滤波前的频谱图title(iir低通滤波器滤波前的频谱)xlabel(频率/hz);ylabel(幅值);subplot(2,1,2)f1=plot(f,abs(f0(1:512); grid %画出滤波后的频谱图title(iir低通滤波器滤波后的频谱)xlabel(频率/hz);ylabel(幅值);end% hints: contents = get(hobject,string) returns popupmenu1 contents as cell array% contentsget(hobject,value) returns selected item from popu

36、pmenu1 % - executes during object creation, after setting all properties.function popupmenu1_createfcn(hobject, eventdata, handles)% hobject handle to popupmenu1 (see gcbo)% eventdata reserved - to be defined in a future version of matlab% handles empty - handles not created until after all createfc

37、ns called % hint: popupmenu controls usually have a white background on windows.% see ispc and computer.if ispc & isequal(get(hobject,backgroundcolor), get(0,defaultuicontrolbackgroundcolor) set(hobject,backgroundcolor,white);end% - executes during object creation, after setting all properties.funct

38、ion axes1_createfcn(hobject, eventdata, handles)% hobject handle to axes1 (see gcbo)% eventdata reserved - to be defined in a future version of matlab% handles empty - handles not created until after all createfcns called % hint: place code in openingfcn to populate axes1四、设计结果及分析（一）、信号分析部分1、单位抽样序列

39、、矩形序列 、三角波序列、反三角波序列、gaussian序列、正弦序列、衰减正弦序列的波形如下图所示：2、对三角波序列和反三角波序列，作n=8点的fft；在和的尾部补零，作n=16点的fft。如下图所示：由上图可知：当求n=8点的fft时，所得的三角波和反三角波的幅频特性相差不大，基本看不出差别。当取n=16时的幅频特性差别明显。对n=16点的幅频特性进行一个间隔一点的采样即得到n=8点的幅频特性。由所学知识可知求fft的公式是（k=0,1.n-1）其中幅度()n的取值越大所得的包络越接近y的dtft的幅度谱。3、固定信号中的参数p=8，令q分别等于2，4，8；固定q=8，令p分别等于8，13

40、，14；如下图所示：由上图可知：当固定p=8,而q在取不同值时它们的时域和幅频度不一样。q的值越大，时域图形的零点越少，图形越“胖”，值越分散。q的值越大，频域图形的零点越多，但其他点的幅值相对增大，发生频谱泄露以及混叠。这是因为值越分散在经过窗函数时，n的值就越大，泄漏也就越大。由上图可知：当固定q=8,而p在取不同值时它们的时域和幅频度不一样。p的值越大，时域图形右移，形状不变。p的值越大，频域图形的频谱泄露就越明显并伴随有频谱的混叠。这是因为谱线右移，经过窗函数时n的取值会相应增大，从而使泄露以及混叠现象越明显。4、对于正弦序列，取数据长度n分别等于8，16，32，分别作n点fft。如下

41、图所示：由上图可知：当n越大时，所求得的fft的幅值就越密，包络曲线就越接近原函数的包络。这是因为n与它的周期紧密相连，序列fft的幅值是关于n/2对称的。5、观察衰减正弦序列的时域和幅频特性，绘出幅频特性曲线，改变采样频率fs,使fs=300hz。如下图所示：采样我们知道是对信号进行时间上的量化，即采样后频谱发生变化。结合我们所学的理论的知识和上图体现出不管是序列的时域还是幅频特性它们的图形都与采样周期相关联，时域的衰减长度及快慢与采样频率联系，fs越大其的衰减长度就越长衰减的也越慢，当fs是300是序列衰减长度是fs=1000的300/1000倍即衰减的速度是fs=1000的300/100

42、0倍。同时它的幅频特性也发生了变化，fs=300hz的一个完整的幅频长度是fs=1000hz的长度发的3/10倍，从上图我们可以看出在实际中我们不一定只能根据信号的宽度来选择滤波器，也可根据所得滤波器去去调整信号宽度减少信号的泄露。6、设有一连续时间信号s(t),其由20hz、220hz和750hz的正弦信号叠加而成，分析确定采样频率及数据分析长度，计算并绘出信号的频谱，指出各个频率份量。如下图所示：根据连续信号的叠加性质，将几个简单信号叠加即可。根据采样定理fs要大于等于2倍的fm，在此可以取fs=2000hz。7、对于iir数字滤波器：使用双线性变换法。可选用bw型、cb型和椭圆型等滤波器

43、类型来进行设计。cb型bw型 椭圆型 由上图可知：在同等条件下（输入的wp=.2,ws=.3,rp=1,rs=34）选择不同的类型滤波器得出来的各项参数也不一样。各项参数如下所示：cb型n = 6wn =0.2000b =1.0e-003 * 0.0464 0.2782 0.6956 0.9274 0.6956 0.2782 0.0464a = 1.0000 -4.8694 10.3807 -12.3367 8.5969 -3.3268 0.5586bw型n =11wn = 0.2183b=1.0e-003*0.0010 0.0111 0.0556 0.1668 0.3336 0.4670 0.4670 0.3336 0.1668 0.0556 0.0111a =1.0000 -6.1863