数字信号处理(胡广书例题作业程序)
- 格式:doc
- 大小:637.51 KB
- 文档页数:48
!!"#$%&’!"#$()*+,-./!!!"!!!"!""!!"#$!!""#!"$"#%$#"#%"%##"#$#"$%&%&’(!""9:!!""+;<&=>?@A+(%!!"BC !!""D&EF+GHIJ !!""%!&"K &"!""#!!!"%""&B9:&"!""+;<%!$"K &!!""#&!!"$!"&B9:&!!""+;<%!’"L !!""G H $M N O A P Q &R S T &U &&!""&B 9:&&!""+;<%!%"VL !!""ST &PGH $MNOAUW &$!""&B9:&$!""+;<%!’!""!!""+;<X;"’"’"YJ %;!"("("!!"!!""#%!!""$%!!"%""$%!!"%!"$%!!"%&"$%!!"%$"$)!!"$""$%!!"$!"$$!!"$&"$!!!"$$"!!"#$%&’()*+,-!!!!!&"&"!""#!!!"%""Z[4\!!""GH"MNO]^&P_Q‘-!UW+&&;<X;"’"’!YJ%;!"("(!!$"&!!""#&!!"$!"Z[4\!!""a7!MNO]^&=_Q‘-&UW+&;< X;"’"’&YJ%;!"("(&!’"L!!""GH$MNOAU!(!""#!!"%$"&PL!(!""Q&RSTU&&!""# !(!%""#!!%"%$"&&&!""+;<X;"’"’$YJ%;!"("($!%"&$!""bIc%’&$!""#!!%"$$"&&;<X;"’"’’YJ%!"#!!"#"!d"’"!:+!!""’!""9:!!%""+;<%"!./01#$2./01345’67(8"!;!"("(’!&"ef !+!""#"!(!!""%!!%"")&=9:!+!""+;<%!$"BC !*!""&!+!""IJ !!""&=ghLiM4\jk%iMldm4\niM odm4\+pq %!’!""!!%""+;<X;"’!’"YJ %;!"(!("!!"!*!""#"!(!!""$!!%"")#"$)%$#"#%"%"$)"#"#$%"###"$%&%&’(&;<X;"’!’!YJ %!!"#$%&’()*+,-#!!&"!+!""#"!(!!""%!!%"")#"$!%$#"#%""%!"#"#$#"$%&%&’(&;<X;"’!’&YJ %;!"(!(&!$"drstuvwx4\!!""&bQL&jk%iMldm4\!*!""yiMo dm4\!+!""zy &{Zuvwx#$jk+i|}~p & !!""#!*!""$!+!"" *!*!""#"!(!!""$!!%"")!+!""#"!(!!""%!!%"$%&") &!*!""y !+!""j !*!""#!*!%""&!+!""#%!+!%""+dm ‘%!"$!#"!"!!" ‘ ’!""&!""#!!""$!!"%""$!!"%!"%!!"&!""#&!%""%!&"&!""#!!"!"%!$"&!""#!!!""%!’"&!""#!!"",-.!"""%!%"&!""#)!!""$*&&*)&*% -%B iM‘ Z * += 0 %!’!""d‘ &!""#!!""$!!"%""$!!"%!"&!" !"!""y !!!""& Y! j p &"!./01#$2./01345’67(8$! K!!""##!"!""$$!!!""‘ d!!""+ F&!""#+(!!"")##!"!""$$!!!""$#!"!"%""$$!!!"%""$#!"!"%!"$$!!!"%!"##(!"!""$!"!"%""$!"!"%!")$$(!!!""$!!!"%""$!!!"%!")&!""##&"!""$$&!!""¡¢‘ !""Z +%r&!""#+(!!"")#!!""$!!"%""$!!"%!"£¤‘ d!!"%,"+ F&,!""Z&,!""#+(!!"%,")#!!"%,"$!!"%,%""$!!"%,%!"¥&!"%,"#!!"%,"$!!"%,%""$!!"%,%!"¦§&!"%,"#+(!!"%,")#&,!""¡¢‘ !"" %!!"d‘ &!""#&!%""&!" !"!""y!!!""& Y! jp &&"!""#+(!"!"")#!"!%""&!!""#+(!!!"")#!!!%""K!!""##!"!""$$!!!""£¤‘ d!!""+ F&!""#+(!!"")##!"!%""$$!!!%""##&"!""$$&!!"" ? +¨©ªZ&"!""y&!!""+«¬&‘ !!"Z +%r&!""#+(!!"")#!!%""£¤‘ d!!"%,"+ F&,!""Z&,!""#+(!!"%,")#!(%!"%,")¥&!"%,"#!(%!"%,")!!"#$%&’()*+,-%!YQ‘ !!" %!&"d‘ &!""#!!"!"&!" !"!""y !!!""& Y! jp & &"!""#+(!"!"")#!"!"!"&!!""#+(!!!"")#!!!"!"K!!""##!"!""$$!!!""£¤‘ d !!""+ F &!""#+(!!"")##!"!"!"$$!!!"!"##&"!""$$&!!""? +¨©ªZ &"!""y &!!""+«¬&‘ !&"Z +% r&!""#+(!!"")#!!"!"£¤‘ d !!"%,"+ F &,!""Z &,!""#+(!!"%,")#!(!"%,"!)¥&!"%,"#!(!"%,"!)¦§&!"%,"#+(!!"%,")#&,!""YQ‘ !&" %!$"d‘ &!""#!!!""&!" !"!""y !!!""& Y! jp & &"!""#+(!"!"")#!!"!""&!!""#+(!!!"")#!!!!""K!!""##!"!""$$!!!""£¤‘ d !!""+ F &!""#+(!!"")#(#!"!""$$!!!"")!’#&"!""$$&!!""¡¢&‘ !$"Z® +% r&!""#+(!!"")#!!!""£¤‘ d !!"%,"+ F &,!""Z &,!""#+(!!"%,")#!!!"%,"¥&!"%,"#!!!"%,""!./01#$2./01345’67(8&!YQ‘ !$" %!’"d‘ &!""#!!"",-.!"""&!" !"!""y !!!""& Y! jp & &"!""#+(!"!"")#!"!"",-.!"""&!!""#+(!!!"")#!!!"",-.!"""K!!""##!"!""$$!!!""£¤‘ d !!""+ F &!""#+(!!"")#(#!"!""$$!!!""),-.!"""##!"!"",-.!"""$$!!!"",-.!"""&!""##&"!""$$&!!""¡¢&‘ !’"Z +% r&!""#+(!!"")#!!"",-.!"""£¤‘ d !!"%,"+ F &,!""Z &,!""#+(!!"%,")#!!"%,",-.!"""¥&!"%,"#!!"%,",-.(!"%,"")¦§&!"%,"’+(!!"%,")#&,!""¡¢&‘ !’" %!%"d‘ &!""#)!!""$*&!" !"!""y !!!""& r )&*% -& Y! j p &&"!""#+(!"!"")#)!"!""$*&!!""#+(!!!"")#)!!!""$*K!!""##!"!""$$!!!""£¤‘ d !!""+ F &!""#+(!!"")#)(#!"!""$$!!!"")$*’#&"!""$$&!!""¡¢&‘ !%"Z® +% r!!"#$%&’()*+,-’!&,!""#+(!!"%,")#)!!"%,"$*¥&!"%,"#)!!"%,"$*¦§&!"%,"#+(!!"%,")#&,!""¡¢&‘ !%" %!"%!#"#"!!" ‘ ’!""&!""#"-$"(-,##!!"%,"&&*-%¯r°+±-%!!"&!""#)!!""$*%!&"&!""#!!""$.!!"$""&&*.% -%!$"&!""#!!"!"%!’"&!""#!!,""&&*,%¯r°+±-%!%"&!""#!!%""%B "²iMZ¡³‘ +²iMZ®¡³‘ += 0 %!’!""&!""#"-$"(-,##!!"%,"&&*-%¯r°+±-%¡%´‘ µs¶w·+ :¸¹"rºµw·y»¼+ !!""&!!"%""&,& !!"%-"&¥yL½+ ¾ &YQ&´‘ Z¡³‘ %!!"&!""#)!!""$*%¡%´‘ µs¶w·+ :¸¹"rºµw·+!!""&¥yL½+ ¾ &Y Q&´‘ Z¡³‘ %!&"&!""#!!""$.!!"$""&&*.% -%¡%´‘ µ¿7w·!""+ : ÀÁ¹r¿7w·!""+ !!""&¥ÂÃÁ¹rL½w·!"$""w+ !!"$""&YQ´‘ Z®¡³‘ %!$"&!""#!!"!"%µ")!w&´‘ ds¶w·"w+ :Ä L½w·"!+ Y¹"&¡¢´‘ %®¡³‘ %!’"&!""#!!,""&&*,%¯r°+±-%XÅÆ,!$"&¿"*#w&‘ + : L½w·,"+ Y¹"&¡¢‘ %®¡³‘ %"!./01#$2./01345’67(8(!¡³‘ %!"&!#"$"!X ÈM‘ ’!""&!""#(-%",###,!!"%,"&&*##&#"&,&#-%"% -%!!"&!""#!#/+,"#&!"%""%#!&!"%!"$!!""%#/+,"#!!"%""&&*#&"#% -%BÉ&ÊËNO F /!""&= ‘ Z Ì"+Ì"+ÍÎZϤ+!’!""ÊËNO FZ‘ µ %ÊËNO4\!!""w+ :%dr´‘ &&ÊËNO F/!""#(-%",###,!!"%,"!!ÐZi M ÑÒr "##&ÓÔ%-+ ÕÓ4\& Z i M 012‘ %Ö r ##&#"&,&#-%"Ä% Õ+ -&YQ´‘ gZÌ"+%!!"bCÈ|pqÉU´‘ +ÊËNO F %pqi ’K !!""#!!""& ‘ + :&!""#/!""& /!""#!#/+,"#/!"%""%#!/!"%!"$!!""%#/+,"#!!"%"""##!/!#"#""#"!/!""#!#/+,"#%#/+,"###/+,"#"#!!/!!"#!#/+,"#(#/+,"#)%#!##!/+,!"#"#&!/!&"#!#/+,"#(#!/+,!"#)%#!(#/+,"#)##&/+,&"#×¢ØÙ&/!""##"/+,""#0!""!!pqÚ’ÛÜbQÝCÞßà!á r 3 â+pqÉ:‘ +ÊËNO F %d‘ + jp ãä3 â&U !"%!#1%"/+,"#$#!1%!"2!1"#!"%#1%"/+,"#"3!1"ã¥UW´‘ T å-4!1"#2!1"3!1"#"%#1%"/+,"#"%!#1%"/+,"#$#!1%!ÊËNO F /!""ZT å-4!1"+æ3 â& Þß+I !’&’"!ç’d #-.#$/012$ZI !’’’""&U /!""#5%"(4!1")##"/+,""#0!""¿§&È|pq!:+h³ZiO+%!!"#$%&’()*+,-)*!8&!""8#8/!"""!!""8#($7,###,/+,"#,!!"%,"#($7,###,/+,"#,8!!"%,"8#6($7,###,/+,"#,#6($7,##8#,8#6($7,##8#8,¿#+"w &‘ + :&!""#/!"""!!""#6"%8#8ÇZ ë+&YQ‘ ZÌ"+%ìz &X³#)"& ‘ Ì"%?íYC+ pq ZÌ" +"î& ë+ ïð ë+ :!4145"%!"’!#"&"!K /!""#-/!#"&/!""&/!!".#-&&!&".&É!""&"!""#/!"""/!""!!"&!!""#/!"""/!"""/!""!’Ék´, È|p q &i Z ñò óô+"îÉ&ÚZ 678974*+:õÎ/+.;½É&ºj !:EF+h³%pqi ’dÆ,!""& &"!""#(79#%7/!"%9"/!9"&bÉ:&"!#"#/!#"/!#"#<&"!""#/!#"/!""$/!""/!#"#"!&"!!"#/!#"/!!"$/!""/!""$/!!"/!#"#"#&"!&"#/!""/!!"$/!!"/!""#$&"!$"#/!!"/!!"#"!!dÆ,!!"&ö÷?&!!""#&"!"""/!""&øùóô+"îbÉ:&!!""#-!=&’$&%&&$$&!"&%&". ú+Ékûüýþÿ!"!:%pqÚ’öºÆ,!""+678974 4Z *>/#"/#%/"(:&#ä´ 4+h³%<!!"!!!"#!!$!!"!!öºÆ,!!"+678974 4ýþÿ!"!:%!"(!/!""$ "’%,!:&K !!""#-!!#"&!!""&!!!"&!!&".#-"&!&&&$.%!""É/!""+!E å-:/!9"%!!"É/!""y !!""+%E å-:/!!9"&=9::/!9"&:/!!9"+;<%。
数字信号处理_胡广书(第三版)_随书光盘关于光盘的使用说明数字信号处理_胡广书(第三版)_随书光盘.rar本光盘共包含六个子目录,其中三个是DSP_FORTRAN, DSP_C和DSP_MATLAB,另外三个是有关习题所需要的数据或文献。
DSP_FORTRAN和DSP_C各含有约40个信号处理的子程序,概括了书中所涉及到的绝大部分算法。
程序分别由FORTRAN语言和C语言编写(MA模型、ARMA模型及最小方差谱估计三个算法只给出了用C语言编写的程序, 没有给出相应的FORTRAN子程序),并在PC机上调试通过。
编译环境是FORTRAN77 V5. 10和TURBO C2. 0。
DSP_MATLAB含有近120多个用MA TLAB编写的信号处理程序,它们是本书各个章节的大部分例题,使用的是MA TLAB6.1。
FORTRAN子程序名称的长度全都是6位,扩展名为.for,C语言子程序的名称全部是7位,由相应的FORTRAN子程序在其名称前加字母m而形成,并将扩展名改为.c。
为了方便读者的使用,光盘中还给出了调用FORTRAN子程序的简单主程序。
读者只需将此主程序和主程序指定的子程序作编译、连接和运行,即可得出相应的结果。
FORTRAN主程序的名称为7位或8位,它是在原FORTRAN子程序前加字母h所构成的,扩展名仍是.for。
h后面的一个数(如果有的话)表示该程序是相应子程序的第几个主程序。
例如,子程序desiir.for是用来设计IIR滤波器的FORTRAN子程序,对应的C程序是mdesiir.c,调用desiir.for 的第一个主程序是h1desiir.for(设计低通IIR DF),依此类推。
用MATLAB编写的程序的名称由“exa”开头,接下来是所在的章、节及例题的序号,如exa010101,指的是第1章第1节(即1.1节)的第1个例题,即例1.1.1。
如果该程序是为了说明某一个m文件的应用,则在上述名称的后面跟一个下划线,再在后面加上所说明的MATLAB文件的名称,如exa011001_rand,即是例1.10.1,该例用来说明rand.m文件的应用。
203⎥⎦⎤⎢⎣⎡---=----)()()()(~01011010z H z z H z z H z H N N m Η (7.6.4b)利用(7.4.9b )的关系,有I ΗΗ210012~=⎥⎦⎤⎢⎣⎡=m m(7.6.5)这样,由(7.6.3)式,CQMFB 的分析滤波器组可以构成仿酉矩阵,其对应的系统也是仿酉系统。
由(7.6.4a )及(7.4.1)式有)1(2det ---=N m z Η(7.6.6)将这一结果代入(7.2.12)式,并令式中的k =0,则⎥⎦⎤⎢⎣⎡-----=--)()()()(0101)1(z H z H z H z H zN m G⎥⎦⎤⎢⎣⎡------=--------)()()()(2010)1(010)1()1(z H z H zz H z H z zN N N (7.6.7) 将(7.6.4a)及(7.6.7)代入(7.2.10)式,有X ΗG X T m m 21ˆ=X ⎥⎦⎤⎢⎣⎡---⎥⎦⎤⎢⎣⎡------=--------------)()()()()()()()(10)1(10)1(00010)1(010)1()1(z H z z H z z H z H z H z H zz H z H z zN N N N N X ⎥⎦⎤⎢⎣⎡-=--10012)1(2N z(7.6.8) 因此,实现了对X 的准确重建。
上面的结论说明,仿酉的调制矩阵m Η直接引出了对)(n x 的准确重建系统,也即CQMFB 。
由(7.6.7)式,可导出0G ,1G 和0H 的关系,即(7.4.2)式。
由上面的讨论可以看出,仿酉滤波器组总是包含了功率互补的关系。
需要指出的是,仿酉系统等效CQMFB ,可以实现准确重建。
但可实现准确重建的系统却并不一定是仿酉的。
现在利用上述讨论的结果来给出仿酉系统的多相表示形式。
记204)()()(20112000z E z z E z H -+= (7.6.9a ) )()()(21112101z E z z E z H -+=(7.6.9b ) )()()(20120010z R z R z z G +=- (7.6.9c ) )()()(21121011z R z R z z G +=-(7.6.9d )式中)(ij ij R E 的下标i 代表0H ,1H 的序号,j 代表多相结构的序号。
数字信号处理课后答案1.2 教材第一章习题解答1. 用单位脉冲序列()n 及其加权和表示题1图所示的序列。
解:()(4)2(2)(1)2()(1)2(2)4(3)0.5(4)2(6)x n nn n n n nnn n 2. 给定信号:25,41()6,040,nnx n n其它(1)画出()x n 序列的波形,标上各序列的值;(2)试用延迟单位脉冲序列及其加权和表示()x n 序列;(3)令1()2(2)x n x n ,试画出1()x n 波形;(4)令2()2(2)x n x n ,试画出2()x n 波形;(5)令3()2(2)x n x n ,试画出3()x n 波形。
解:(1)x(n)的波形如题2解图(一)所示。
(2)()3(4)(3)(2)3(1)6()6(1)6(2)6(3)6(4)x n nnnn n n n n n (3)1()x n 的波形是x(n)的波形右移2位,在乘以2,画出图形如题2解图(二)所示。
(4)2()x n 的波形是x(n)的波形左移2位,在乘以2,画出图形如题2解图(三)所示。
(5)画3()x n 时,先画x(-n)的波形,然后再右移2位,3()x n 波形如题2解图(四)所示。
3. 判断下面的序列是否是周期的,若是周期的,确定其周期。
(1)3()cos()78x n A n,A 是常数;(2)1()8()j n x n e 。
解:(1)3214,73w w ,这是有理数,因此是周期序列,周期是T=14;(2)12,168ww,这是无理数,因此是非周期序列。
5. 设系统分别用下面的差分方程描述,()x n 与()y n 分别表示系统输入和输出,判断系统是否是线性非时变的。
(1)()()2(1)3(2)y n x n x n x n;(3)0()()y n x n n ,0n 为整常数;(5)2()()y n x n ;(7)0()()n m y n x m 。
第5章信号的抽取与插值5.1前言至今,我们讨论的信号处理的各种理论、算法及实现这些算法的系统都是把抽样频率f视为恒定值,即在一个数字系统中只有一个抽样率。
但是,在实际工作中,我们经常会s遇到抽样率转换的问题。
一方面,要求一个数字系统能工作在“多抽样率(multirate)”状态,以适应不同抽样信号的需要;另一方面,对一个数字信号,要视对其处理的需要及其自身的特征,能在一个系统中以不同的抽样频率出现。
例如:1. 一个数字传输系统,即可传输一般的语音信号,也可传输播视频信号,这些信号的频率成份相差甚远,因此,相应的抽样频率也相差甚远。
因此,该系统应具有传输多种抽样率信号的能力,并自动地完成抽样率的转换;2. 如在音频世界,就存在着多种抽样频率。
得到立体声声音信号(Studio work)所用的抽样频率是48kHz,CD产品用的抽样率是44.1kHz,而数字音频广播用的是32kHz[15]。
3. 当需要将数字信号在两个具有独立时钟的数字系统之间传递时,则要求该数字信号的抽样率要能根据时钟的不同而转换;4.对信号(如语音,图象)作谱分析或编码时,可用具有不同频带的低通、带通及高通滤波器对该信号作“子带”分解,对分解后的信号再作抽样率转换及特征提取,以实现最大限度减少数据量,也即数据压缩的目的;5. 对一个信号抽样时,若抽样率过高,必然会造成数据的冗余,这时,希望能在该数字信号的基础上将抽样率减下来。
以上几个方面都是希望能对抽样率进行转换,或要求数字系统能工作在多抽样率状态。
近20年来,建立在抽样率转换理论及其系统实现基础上的“多抽样率数字信号处理”已成为现代信号处理的重要内容。
“多抽样率数字信号处理”的核心内容是信号抽样率的转换及滤波器组。
减少抽样率以去掉过多数据的过程称为信号的“抽取(decimatim)”,增加抽样率以增加数据的过程称为信号的“插值(interpolation)。
抽取、插值及其二者相结合的使用便可实现信号抽样率的转换。
1、%---filter求卷积,B(Z)/A(Z)=H(Z),已知B(Z)和A(Z),求y(n)=x(n)*h(n)----- clear;x=ones(100);t=1:100;b=[.001836,.007344,.011016,.007374,.001836];a=[1,-3.0544,3.8291,-2.2925,.55075];%y=filter(b,a,x);% 求所给系统的输出,本例实际上是求所给系统的阶跃响应;plot(t,x,'r.',t,y,'k-');grid on;ylabel('x(n) and y(n)')xlabel('n')1、%---filter求卷积,B(Z)/A(Z)=H(Z),已知B(Z)和A(Z),求y(n)=x(n)*h(n)----- clear;x=ones(100);t=1:100;b=[.001836,.007344,.011016,.007374,.001836];a=[1,-3.0544,3.8291,-2.2925,.55075];%y=filter(b,a,x);% 求所给系统的输出,本例实际上是求所给系统的阶跃响应;plot(t,x,'r.',t,y,'k-');grid on;ylabel('x(n) and y(n)')xlabel('n')第一章产生信号,求卷积和自相关函数1、%信号产生n=0:100;%工频f0=50;A=220;fs=400;x1=A*sin(2*pi*f0*n/fs);subplot(321);plot(n,x1);xlabel('n');ylabel('x1(n)') ;grid on;%率减正弦f0=2;A=2;alf=0.5;fs=16;x2=A*exp(-alf*n/fs).*sin(2*pi*f0*n/fs);subplot(323);plot(n,x2);xlabel('n');ylabel('x2(n)') ;grid on;%谐波信号f0=5;A1=1.0;A2=0.5;A3=0.2;fs=100;x3=A1*sin(2*pi*f0*n/fs)+A2*sin(2*pi*2*f0*n/fs)+A3*sin(2*pi*3*f0*n/fs); subplot(322);plot(n,x3);xlabel('n');ylabel('x3(n)') ;grid on;%哈明窗f0=10;fs=1000;x4=0.54-0.46*cos(2*pi*f0*n/fs);subplot(324);plot(n,x4);xlabel('n');ylabel('x4(n)') ;grid on;%采样n=-50:50;f0=10;fs=400;w=2*pi*f0*n/fs;x5=sinc(w);subplot(325);plot(n,x5);xlabel('n');ylabel('x5(n)') ;grid on;2、% 产生均匀分布的白噪信号,使均值为0,功率为p%-----------------------------------------------------------------clear;p=0.01;N=50000;u=rand(1,N);u=u-mean(u);a=sqrt(12*p);u1=u*a;power_u1=dot(u1,u1)/Nsubplot(211)plot(u1(1:200));grid on;ylabel('u(n)')xlabel('n')3、% 产生高斯分布的白噪信号,使功率为p,并观察数据分布的直方图%-----------------------------------------------------------------clear;p=0.1;N=500000;u=randn(1,N);a=sqrt(p);u=u*a;power_u=var(u);plot(u(1:200));grid on;ylabel('u(n)');xlabel('n')subplot(212)hist(u,50);grid on;ylabel('histogram of u(n)');4、% 产生一sinc 函数;%-----------------------------------------------------------------clear;n=200;stept=4*pi/n;t=-2*pi:stept:2*pi;y=sinc(t);plot(t,y,t,zeros(size(t)));ylabel('sinc(t)');xlabel('t=-2*pi~2*pi');grid on;5、% 产生一chirp 信号;% chirp(T0,F0,T1,F1):% T0: 信号的开始时间;F0:信号在T0时的瞬时频率,单位为Hz;% T1: 信号的结束时间;F1:信号在T1时的瞬时频率,单位为Hz;%-----------------------------------------------------------------clear;t=0:0.001:1;x=chirp(t,0,1,125);plot(t,x);ylabel('x(t)')xlabel('t')6、% 计算两个序列的线性卷积;%-----------------------------------------------------------------clear;N=5;%第一个序列的长度M=6;%第二个序列的长度L=N+M-1;h=[6,2,3,6,4,2];y=conv(x,h);nx=0:N-1;nh=0:M-1;ny=0:L-1;subplot(231);%绘制xstem(nx,x,'.k');xlabel('n');ylabel('x(n)');grid on;subplot(232);%绘制hstem(nh,h,'.k');xlabel('n');ylabel('h(n)');grid on;subplot(233);%绘制卷积stem(ny,y,'.k');xlabel('n');ylabel('y(n)');grid on;7、% 求两个序列的互相关函数,或一个序列的自相关函数;%-----------------------------------------------------------------clear;N=500;p1=1;p2=0.1;f=1/8;Mlag=50;%自相关的单边长度u=randn(1,N);n=[0:N-1];s=sin(2*pi*f*n);% 混有高斯白噪的正弦信号的自相关u1=u*sqrt(p1);%高斯白噪声x1=u1(1:N)+s;%混合信号rx1=xcorr(x1,Mlag,'biased');%自相关,无偏估计subplot(221);plot(x1(1:Mlag));title('信号x1');xlabel('n');ylabel('x1(n)');grid on;subplot(223);plot((-Mlag:Mlag),rx1);title('x1自相关');grid on;xlabel('m');ylabel('rx1(m)');% 高斯白噪功率由原来的p1减少为p2,再观察混合信号的自相关u2=u*sqrt(p2);%改变高斯白噪声x2=u2(1:N)+s;%新的混合信号rx2=xcorr(x2,Mlag,'biased');subplot(222);plot(x2(1:Mlag));title('信号x2');xlabel('n');ylabel('x2(n)');grid on;subplot(224);plot((-Mlag:Mlag),rx2);title('x2自相关');grid on;xlabel('m');ylabel('rx2(m)');8、%求序列的自相关函数clearN=500;Mlag=50;%单边长度nx=0:N-1;x=exp(-nx*0.1);rx=xcorr(x,Mlag,'biased');nrx=-Mlag:Mlag;%自相关序列的程度subplot(211);plot(nx,x);xlabel('n');ylabel('x(n)') ;grid on; subplot(212);plot(nrx,rx);xlabel('n');ylabel('rx(n)') ;grid on;9、%正弦加白噪声,自相关p=0.1;N=5000;Mlag=100;u=rand(1,N);u=u-mean(u);a=sqrt(12*p);u=u*a;power_u=dot(u,u)/Nnx=1:1000;x=1.414*sin(nx*pi/16.0);x1=x(1:1000)+5*u(1:1000);rx=xcorr(x1,Mlag,'biased');nrx=-Mlag:Mlag;subplot(211);plot(x1(1:200));xlabel('n');ylabel('x(n)') ;grid on; subplot(212);plot(nrx,rx);xlabel('n');ylabel('rx(n)') ;grid on;1、%产生信号,求卷积,FFT,求平均clear all;N=1024;%采样点数fs=100.0;%采样频率alf1=-1.0; f1=5;alf2=-1.5; f2=8;alf3=-0.7; f3=10;%产生x和w两个信号%产生xu=rand(1,N);u=u-mean(u);%均值为0的白噪声t=[0:1/fs:(N-1)/fs];x=1.0*exp(alf1*t).*sin(2*pi*f1*t)+1.0*exp(alf2*t).*sin(2*pi*f2*t)+1.0*exp(alf3*t).*sin(2*pi*f3 *t);x=x+u;x=x/max(x);%产生walf4=-1.0;w=1.0*exp(alf4*t);%x=x-mean(x);figure(1);subplot(211);plot(t,x,t,w,'r');title('x(t) w(t)')% 应用FFT 求频谱;f=0:fs/N:fs/N*(N-1);X=fft(x,N);X=abs(X);% X=20*log10(X);subplot(212);plot(f(1:N/2),X(1:N/2));title('x(t)频谱')% stem(f,X,'.');grid on;xlabel('Hz')%求卷积-----------------------------------------------y=x.*w;%时域相乘,频域卷积figure(2);subplot(211);plot(t,y);title('正弦加白噪声后与w时域相乘')Y=fft(y,N);Y=abs(Y);subplot(212);plot(f(1:N/2),Y(1:N/2));title('正弦加白噪声后与w时域相乘的FFT')xlabel('Hz')xlabel('Hz')%平均1000次----------------------------------------------figure(3);Y=zeros(1,N);u=rand(1,1000*N);u=u-mean(u);%零均值白噪声for i=0:999x=1.0*exp(alf1*t).*sin(2*pi*f1*t)+1.0*exp(alf2*t).*sin(2*pi*f2*t)+1.0*exp(alf3*t).*sin(2*pi*f3 *t);x=x+10*u(1+i*N:i*N+N); x=x/max(x);X=fft(x,N);X=abs(X);Y=Y+X;endY=Y/1000;subplot(211);plot(f(1:N/2),Y(1:N/2));title('正弦加白噪声的FFT 1000次平均')xlabel('Hz')Y=zeros(1,N);for i=0:999x=1.0*exp(alf1*t).*sin(2*pi*f1*t)+1.0*exp(alf2*t).*sin(2*pi*f2*t)+1.0*exp(alf3*t).*sin(2*pi*f3 *t);x=x+10*u(1+i*N:i*N+N); x=x/max(x);x=x.*w;X=fft(x,N);X=abs(X);Y=Y+X;endY=Y/1000;subplot(212);plot(f(1:N/2),Y(1:N/2));title('正弦加白噪声后与w时域相乘的FFT 1000次平均') xlabel('Hz')24681012-1-0.500.51x(t) w(t)05101520253035404550102030x(t)频谱Hz024681012-0.50.51正弦加白噪声后与w 时域相乘05101520253035404550510正弦加白噪声后与w 时域相乘的FFTHz51015202530354045501012141618正弦加白噪声的FFT 1000次平均Hz510152025303540455023456正弦加白噪声后与w 时域相乘的FFT 1000次平均Hz2、clear all;%三个正弦信号相加,分段函数,进行频谱分析 % 产生三个正弦相加的函数; N=512;f0=10;fs=100.0; t=[0:N-1];x=1.0*sin(2*pi*f0*t/fs)+1.0*sin(2*pi*2*f0*t/fs)+1.0*sin(2*pi*3*f0*t/fs); subplot(211);plot(t(1:N),x(1:N));title('x(t)')%加窗w=1-1.93*cos(2*pi*t/N)+1.29*cos(4*pi*t/N)-0.388*cos(6*pi*t/N)+0.0322*cos(8*pi*t/N); %w=1.0-cos(2*pi*t/N);x=w.*x;%加窗等于时域点乘 % 应用FFT 求频谱; f=0:fs/N:fs/N*(N-1);X=fft(x,N);%先点乘再进行傅里叶变换 X=abs(X)/N; subplot(212);plot(f(1:N/2),X(1:N/2));title('x(t)加窗之后的傅里叶变换') xlabel('Hz')100200300400500600-4-2024x(t)510152025303540455000.20.40.60.8x(t)加窗之后的傅里叶变换Hz%分段函数----------------------------------------------- M=170;L=N-2*M;x(1:M)=1.0*sin(2*pi*f0*t(1:M)/fs);x(M+1:2*M)=1.0*sin(2*pi*2*f0*t(1:M)/fs); x(2*M+1:N)=1.0*sin(2*pi*3*f0*t(1:L)/fs); figure;subplot(211);plot(t(1:N),x(1:N));title('x(t)为分段函数')%w=1-1.93*cos(2*pi*t/M)+1.29*cos(4*pi*t/M)-0.388*cos(6*pi*t/M)+0.0322*cos(8*pi*t/M); %x(1:M)=w(1:M).*x(1:M);X=fft(x,N); X=abs(X)/N; subplot(212);plot(f(1:N/2),X(1:N/2));title('x(t)频谱分析') xlabel('Hz')100200300400500600-1-0.500.51x(t)为分段函数510152025303540455000.050.10.150.2x(t)频谱分析Hz3、%采样长度不同对FFT 的影响---------------------------------------------------- clear all;% 观察数据长度N 的变化对DTFT 分辨率的影响 f1=2;f2=2.02;f3=2.07;fs=10; w=2*pi/fs;N=256;%N=256x=sin(w*f1*(0:N-1))+sin(w*f2*(0:N-1))+sin(w*f3*(0:N-1)); f=0:fs/N:fs/2-1/N; X=fft(x); X=abs(X); subplot(221);plot(f(45:60),X(45:60));title('N=256');grid on; subplot(223)plot(f(1:N/2),X(1:N/2));title('N=256');grid on; xlabel('Hz') %N=N*4;%N=1024x=sin(w*f1*(0:N-1))+sin(w*f2*(0:N-1))+sin(w*f3*(0:N-1)); f=0:fs/N:fs/2-1/N; X=fft(x); X=abs(X); subplot(222)plot(f(45*4:4*60),X(4*45:4*60));title('N=1024');grid on;xlabel('Hz')1.52 2.5050100150N=256024650100150N=256Hz1.52 2.5020*******N=1024Hz4、%补零的影响------------------------------------------------------------------- clear;% 计算长度为N 的原始信号的DTFT f1=2.67;f2=3.75;f3=6.75;fs=20;w=2*pi/fs; N=16;x=sin(w*f1*(0:N-1))+sin(w*f2*(0:N-1)+pi/2)+sin(w*f3*(0:N-1)); f=0:fs/N:fs/2-1/N; X=fft(x); X=abs(X);f=fs/N*(0:N/2-1); subplot(221)stem(f,X(1:N/2),'.');title('不补零');grid on; xlabel('Hz')% 在数据末补N 个零 x(N:2*N-1)=0;X=fft(x); X=abs(X); f=fs*(0:N-1)/(2*N); subplot(222)stem(f,X(1:N),'.');title('补N 个零');grid on; xlabel('Hz')% 在数据末补7*N 个零 x(N:8*N-1)=0;X=fft(x); X=abs(X); f=fs*(0:4*N-1)/(8*N); subplot(223)stem(f,X(1:4*N),'.');title('补7N 个零');grid on; xlabel('Hz')% 在数据末补29*N 个零 x(N:30*N-1)=0; X=fft(x); X=abs(X); f=fs*(0:15*N-1)/(30*N); subplot(224)plot(f,X(1:15*N));title('补29N 个零');grid on; xlabel('Hz')51002468不补零Hz 051002468补N 个零Hz 0510510补7N 个零Hz510510补29N 个零Hz5、%x(n)是两个正弦信号和一个白噪声相加,FFT 和IFFT ------------------ ----------- clear all;% 产生两个正弦加白噪声; N=256;f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs;x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N);% 应用FFT 求频谱; subplot(3,1,1);plot(x(1:N/4));title('x(n)') f=-0.5:1/N:0.5-1/N; X=fft(x); y=ifft(X); subplot(3,1,2);plot(f,fftshift(abs(X)));title('fft(x)') subplot(3,1,3);plot(real(x(1:N/4))); title('x(n)实部')010203040506070-10010x(n)-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50500fft(x)010203040506070-10010x(n)实部6、%fftfilt 和conv 比较 n=0:2; h=1./2.^n;for i=1:51x(i)=(i-1)/5; end for i=52:100x(i)=20-(i-1)/5; endy=fftfilt(h,x); % y=x*h z=conv(h,x);title('x(t)') hold; plot(x);subplot(312)title('f ftfilt 叠接相加法') hold; plot(y);subplot(313) plot(z);axis([0,100,0,20]); title ('conv 卷积')01020304050607080901000510x(t)010203040506070809010001020fftfilt 叠接相加法010********607080901001020conv 卷积7、clear;%补零后对频谱的影响% 计算长度为N 的原始信号的DFT f0=50;%信号频率 fs=200;%采样频率 N=16;%抽样点数w=2*pi/fs;%数字角频率 x=sin(w*f0*(0:N-1)); X=fft(x); X=abs(X);f=fs/N*(0:N/2-1);%N/2的数据stem(f,X(1:N/2),'.');title('未补零');grid on; xlabel('Hz')% 在数据末补N 个零 x(N:2*N-1)=0;X=fft(x); X=abs(X); f=fs*(0:N-1)/(2*N); subplot(212)stem(f,X(1:N),'.');title('补N 个零');grid on; xlabel('Hz')10203040506070809002468未补零Hz 010203040506070809010002468补N 个零Hz8、%抽样频率不同的影响% 计算长度为N 的原始信号的DFT % fs=100f0=50;fs=100;w=2*pi/fs; N=16;x=sin(w*f0*(0:N-1)); X=fft(x); X=abs(X);f=fs/N*(0:N/2-1); subplot(221)stem(f,X(1:N/2),'.');grid on; xlabel('Hz')如对您有帮助,欢迎下载支持,谢谢!% fs=150f0=50;fs=150;w=2*pi/fs; N=16;x=sin(w*f0*(0:N-1)); X=fft(x); X=abs(X);f=fs/N*(0:N/2-1); subplot(222)stem(f,X(1:N/2),'.');grid on; xlabel('Hz')% fs=200f0=50;fs=200;w=2*pi/fs; N=16;x=sin(w*f0*(0:N-1)); X=fft(x); X=abs(X);f=fs/N*(0:N/2-1); subplot(223)stem(f,X(1:N/2),'.');grid on; xlabel('Hz')02040600.511.5x 10-14Hz0204060802468Hz501002468Hz9、%周期延拓 clear all;M=128;N=1024; alf=-3.0;f0=10;fs=100.0;%u=rand(1,5000);u=u-mean(u);t=[0:1/fs:(M-1)/fs];x1=1.0*exp(alf*t).*sin(2*pi*f0*t);t=[0:1/fs:(N-1)/fs]; for i=0:7x(M*i+1:M*i+M)=x1(1:M); endx=x-mean(x); subplot(211); plot(t,x);% 应用FFT 求频谱; f=0:fs/N:fs/N*(N-1); X=fft(x,N); X=abs(X);%X=20*log10(X); subplot(212);plot(f(1:N/2),X(1:N/2)); % stem(f,X,'.');grid on; xlabel('Hz')024681012-1-0.500.510510152025303540455050100150Hz10、%正弦加白噪声并FFT 频谱分析% 产生4096点的两个正弦加白噪声; N=4096;f1=2.1;f2=2.2;fs=5; a1=5;a2=3; w=2*pi/fs;x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+10*randn(1,N);% 应用FFT 求频谱; f=fs/N*(0:N/2-1); X=fft(x,N); X=abs(X);subplot(211);stem(f,X(1:N/2),'.');title('x(n)');grid on; xlabel('Hz')f=-0.5:1/N:0.5-1/N; subplot(212);stem(f,fftshift(X),'.');title('FFT x(n)');grid on;00.511.522.5500010000x(n)Hz-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50500010000FFT x(n)11、clear;%补零对频谱分析的影响f1=10.8;f2=11.75;f3=12.55;fs=40;w=2*pi/fs;N=64;x=sin(w*f1*(0:N-1))+sin(w*f2*(0:N-1)+pi/2)+sin(w*f3*(0:N-1)); X=fft(x);X=abs(X);f=fs/N*(0:N/2-1);subplot(221)stem(f,X(1:N/2),'.');title('补0个零');grid on;xlabel('Hz')% 在数据末补3N个零x(N:4*N-1)=0;X=fft(x); X=abs(X);f=fs*(0:2*N-1)/(4*N);subplot(222)stem(f,X(1:2*N),'.');title('补3N个零');grid on;xlabel('Hz')% 在数据末补7*N个零x(N:8*N-1)=0;X=fft(x); X=abs(X);f=fs*(0:4*N-1)/(8*N);subplot(223)stem(f,X(1:4*N),'.');title('补7N个零');grid on;xlabel('Hz')% 在数据末补15*N个零x(N:16*N-1)=0;X=fft(x); X=abs(X);f=fs*(0:8*N-1)/(16*N);subplot(224)plot(f,X(1:8*N));title('补15N个零');grid on;xlabel('Hz')5101520010203040补0个零Hz 05101520010203040补3N 个零Hz 05101520010203040补7N 个零Hz5101520010203040补15N 个零Hz12、%窗函数的影响 clear;%窗函数对频谱的影响f1=10.8;f2=11.75;f3=12.55; fs=40; w=2*pi/fs;N=1024;t=[0:N-1]; f=fs/N*(0:N/2-1); %矩形窗x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); X=fft(x);X=2*abs(X)/N; subplot(221)stem(f,X(1:N/2),'.');grid on; xlabel('Hz') title('矩形窗') %升余弦窗x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); wt=0.5-0.5*cos(2*pi*t/N); x=wt.*x; X=fft(x);X=2*2*abs(X)/N; subplot(222)stem(f,X(1:N/2),'.');grid on;xlabel('Hz') title('升余弦窗') %平顶窗x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t);wt=1-1.93*cos(2*pi*t/N)+1.29*cos(4*pi*t/N)-0.388*cos(6*pi*t/N)+0.0322*cos(8*pi*t/N); x=wt.*x; X=fft(x);X=2*abs(X)/N; subplot(223)stem(f,X(1:N/2),'.');grid on; xlabel('Hz') title('平顶窗') %改进升余弦窗x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); wt=0.54-0.46*cos(2*pi*t/N); x=wt.*x; X=fft(x);X=1.852*2*abs(X)/N; subplot(224)stem(f,X(1:N/2),'.');grid on; xlabel('Hz')title('改进升余弦窗')051015200.51Hz 矩形窗051015200.51Hz升余弦窗0510152000.511.5Hz平顶窗051015200.51Hz改进升余弦窗13、%频谱分析,平均clear all;%产生正弦加白噪声信号,求频谱,平均% 产生两个正弦加白噪声;N=1024;f0=10;fs=100.0;u=rand(1,N);u=u-mean(u);t=[0:1/fs:(N-1)/fs];x=1.0*sin(2*pi*f0*t)+0.1*sin(2*pi*2*f0*t)+10.0*u;% 应用FFT 求频谱;f=0:fs/N:fs/N*(N-1);X=fft(x,N);X=abs(X);X=20*log10(X);subplot(211);plot(f(1:N/2),X(1:N/2));title('平均前')% stem(f,X,'.');grid on;xlabel('Hz')Y=zeros(1,N);u=rand(1,1000*N);u=u-mean(u);for i=0:999x=1.0*sin(2*pi*f0*t)+0.1*sin(2*pi*2*f0*t)+10*u(1+i*N:i*N+N);X=fft(x,N);X=abs(X);Y=Y+X;endY=Y/1000;Y=20*log10(Y);subplot(212);plot(f(1:N/2),Y(1:N/2));title('平均后')axis([0,50,0,60]);xlabel('Hz')05101520253035404550204060平均前Hz 05101520253035404550204060平均后Hz1高通 2带阻3低通 4带通 5切比雪夫 带通1、1高通,双线性变换法,巴特沃斯(两种方法)% to design 高通high-pass DF with s=2/Ts[(z-1)/(z+1)]%-给出:通带下限频率,阻带上限频率,通带衰减,阻带衰减,采样频率----------------------------------------------- %---------------------- clear allwp=.8*pi;%通带下限频率 ws=.44*pi;%阻带上限频率 rp=3;%通带衰减 rs=20;%阻带衰减 Fs=2000;%采样频率%设计模拟滤波器% Firstly to finish frequency prewarping;wap=2*Fs*tan(wp/2);%通带截止频率,模拟角频率,预畸变 was=2*Fs*tan(ws/2);%阻带截止频率,模拟角频率,预畸变 [n,wn]=buttord(wap,was,rp,rs,'s');%求取模拟低通滤波器阶数,n 是模拟低通滤波器的阶次,巴特沃斯滤波器阶数的选择,(最小阶数的选择)[z,p,k]=buttap(n);%设计模拟低通滤波器,极点,零点,增益 [b,a]=zp2tf(z,p,k);%零-极点增益模型转换为传递函数模型 [bt,at]=lp2hp(b,a,wap);%模拟低通滤波器转换成高通滤波器,G(p)转换成H (s ),wap 为低通的截止频率 %模拟转换成数字% Note: z=(2/ts)(z-1)/(z+1);ts=1,that is 2fs=1,fs=0.5; [bz,az]=bilinear(bt,at,Fs);%实现双线性变换,由模拟滤波器H (s )得到数字滤波器H (Z ),bz,az 分别是H (Z )的分子分母多项式的系数 [h,w]=freqz(bz,az,256,1);%求离散系统频响特性,H 包含了离散系统频响在 0~pi 范围内N 个频率等分点的值,w 则包含了范围内N 个频率等分点。
1、%---filter求卷积,B(Z)/A(Z)=H(Z),已知B(Z)和A(Z),求y(n)=x(n)*h(n)----- clear;x=ones(100);t=1:100;b=[.001836,.007344,.011016,.007374,.001836];a=[1,-3.0544,3.8291,-2.2925,.55075];%y=filter(b,a,x);% 求所给系统的输出,本例实际上是求所给系统的阶跃响应;plot(t,x,'r.',t,y,'k-');grid on;ylabel('x(n) and y(n)')xlabel('n')1、%---filter求卷积,B(Z)/A(Z)=H(Z),已知B(Z)和A(Z),求y(n)=x(n)*h(n)----- clear;x=ones(100);t=1:100;b=[.001836,.007344,.011016,.007374,.001836];a=[1,-3.0544,3.8291,-2.2925,.55075];%y=filter(b,a,x);% 求所给系统的输出,本例实际上是求所给系统的阶跃响应;plot(t,x,'r.',t,y,'k-');grid on;ylabel('x(n) and y(n)')xlabel('n')第一章产生信号,求卷积和自相关函数1、%信号产生n=0:100;%工频f0=50;A=220;fs=400;x1=A*sin(2*pi*f0*n/fs);subplot(321);plot(n,x1);xlabel('n');ylabel('x1(n)') ;grid on;%率减正弦f0=2;A=2;alf=0.5;fs=16;x2=A*exp(-alf*n/fs).*sin(2*pi*f0*n/fs);subplot(323);plot(n,x2);xlabel('n');ylabel('x2(n)') ;grid on;%谐波信号f0=5;A1=1.0;A2=0.5;A3=0.2;fs=100;x3=A1*sin(2*pi*f0*n/fs)+A2*sin(2*pi*2*f0*n/fs)+A3*sin(2*pi*3*f0*n/fs); subplot(322);plot(n,x3);xlabel('n');ylabel('x3(n)') ;grid on;%哈明窗f0=10;fs=1000;x4=0.54-0.46*cos(2*pi*f0*n/fs);subplot(324);plot(n,x4);xlabel('n');ylabel('x4(n)') ;grid on;%采样n=-50:50;f0=10;fs=400;w=2*pi*f0*n/fs;x5=sinc(w);subplot(325);plot(n,x5);xlabel('n');ylabel('x5(n)') ;grid on;2、% 产生均匀分布的白噪信号,使均值为0,功率为p%-----------------------------------------------------------------clear;p=0.01;N=50000;u=rand(1,N);u=u-mean(u);a=sqrt(12*p);u1=u*a;power_u1=dot(u1,u1)/Nsubplot(211)plot(u1(1:200));grid on;ylabel('u(n)')xlabel('n')3、% 产生高斯分布的白噪信号,使功率为p,并观察数据分布的直方图%-----------------------------------------------------------------clear;p=0.1;N=500000;u=randn(1,N);a=sqrt(p);u=u*a;power_u=var(u);plot(u(1:200));grid on;ylabel('u(n)');xlabel('n')subplot(212)hist(u,50);grid on;ylabel('histogram of u(n)');4、% 产生一sinc 函数;%-----------------------------------------------------------------clear;n=200;stept=4*pi/n;t=-2*pi:stept:2*pi;y=sinc(t);plot(t,y,t,zeros(size(t)));ylabel('sinc(t)');xlabel('t=-2*pi~2*pi');grid on;5、% 产生一chirp 信号;% chirp(T0,F0,T1,F1):% T0: 信号的开始时间;F0:信号在T0时的瞬时频率,单位为Hz;% T1: 信号的结束时间;F1:信号在T1时的瞬时频率,单位为Hz;%-----------------------------------------------------------------clear;t=0:0.001:1;x=chirp(t,0,1,125);plot(t,x);ylabel('x(t)')xlabel('t')6、% 计算两个序列的线性卷积;%-----------------------------------------------------------------clear;N=5;%第一个序列的长度M=6;%第二个序列的长度L=N+M-1;h=[6,2,3,6,4,2];y=conv(x,h);nx=0:N-1;nh=0:M-1;ny=0:L-1;subplot(231);%绘制xstem(nx,x,'.k');xlabel('n');ylabel('x(n)');grid on;subplot(232);%绘制hstem(nh,h,'.k');xlabel('n');ylabel('h(n)');grid on;subplot(233);%绘制卷积stem(ny,y,'.k');xlabel('n');ylabel('y(n)');grid on;7、% 求两个序列的互相关函数,或一个序列的自相关函数;%-----------------------------------------------------------------clear;N=500;p1=1;p2=0.1;f=1/8;Mlag=50;%自相关的单边长度u=randn(1,N);n=[0:N-1];s=sin(2*pi*f*n);% 混有高斯白噪的正弦信号的自相关u1=u*sqrt(p1);%高斯白噪声x1=u1(1:N)+s;%混合信号rx1=xcorr(x1,Mlag,'biased');%自相关,无偏估计subplot(221);plot(x1(1:Mlag));title('信号x1');xlabel('n');ylabel('x1(n)');grid on;subplot(223);plot((-Mlag:Mlag),rx1);title('x1自相关');grid on;xlabel('m');ylabel('rx1(m)');% 高斯白噪功率由原来的p1减少为p2,再观察混合信号的自相关u2=u*sqrt(p2);%改变高斯白噪声x2=u2(1:N)+s;%新的混合信号rx2=xcorr(x2,Mlag,'biased');subplot(222);plot(x2(1:Mlag));title('信号x2');xlabel('n');ylabel('x2(n)');grid on;subplot(224);plot((-Mlag:Mlag),rx2);title('x2自相关');grid on;xlabel('m');ylabel('rx2(m)');8、%求序列的自相关函数clearN=500;Mlag=50;%单边长度nx=0:N-1;x=exp(-nx*0.1);rx=xcorr(x,Mlag,'biased');nrx=-Mlag:Mlag;%自相关序列的程度subplot(211);plot(nx,x);xlabel('n');ylabel('x(n)') ;grid on; subplot(212);plot(nrx,rx);xlabel('n');ylabel('rx(n)') ;grid on;9、%正弦加白噪声,自相关p=0.1;N=5000;Mlag=100;u=rand(1,N);u=u-mean(u);a=sqrt(12*p);u=u*a;power_u=dot(u,u)/Nnx=1:1000;x=1.414*sin(nx*pi/16.0);x1=x(1:1000)+5*u(1:1000);rx=xcorr(x1,Mlag,'biased');nrx=-Mlag:Mlag;subplot(211);plot(x1(1:200));xlabel('n');ylabel('x(n)') ;grid on; subplot(212);plot(nrx,rx);xlabel('n');ylabel('rx(n)') ;grid on;1、%产生信号,求卷积,FFT,求平均clear all;N=1024;%采样点数fs=100.0;%采样频率alf1=-1.0; f1=5;alf2=-1.5; f2=8;alf3=-0.7; f3=10;%产生x和w两个信号%产生xu=rand(1,N);u=u-mean(u);%均值为0的白噪声t=[0:1/fs:(N-1)/fs];x=1.0*exp(alf1*t).*sin(2*pi*f1*t)+1.0*exp(alf2*t).*sin(2*pi*f2*t)+1.0*exp(alf3*t).*sin(2*pi*f3 *t);x=x+u;x=x/max(x);%产生walf4=-1.0;w=1.0*exp(alf4*t);%x=x-mean(x);figure(1);subplot(211);plot(t,x,t,w,'r');title('x(t) w(t)')% 应用FFT 求频谱;f=0:fs/N:fs/N*(N-1);X=fft(x,N);X=abs(X);% X=20*log10(X);subplot(212);plot(f(1:N/2),X(1:N/2));title('x(t)频谱')% stem(f,X,'.');grid on;xlabel('Hz')%求卷积-----------------------------------------------y=x.*w;%时域相乘,频域卷积figure(2);subplot(211);plot(t,y);title('正弦加白噪声后与w时域相乘')Y=fft(y,N);Y=abs(Y);subplot(212);plot(f(1:N/2),Y(1:N/2));title('正弦加白噪声后与w时域相乘的FFT')xlabel('Hz')xlabel('Hz')%平均1000次----------------------------------------------figure(3);Y=zeros(1,N);u=rand(1,1000*N);u=u-mean(u);%零均值白噪声for i=0:999x=1.0*exp(alf1*t).*sin(2*pi*f1*t)+1.0*exp(alf2*t).*sin(2*pi*f2*t)+1.0*exp(alf3*t).*sin(2*pi*f3 *t);x=x+10*u(1+i*N:i*N+N); x=x/max(x);X=fft(x,N);X=abs(X);Y=Y+X;endY=Y/1000;subplot(211);plot(f(1:N/2),Y(1:N/2));title('正弦加白噪声的FFT 1000次平均')xlabel('Hz')Y=zeros(1,N);for i=0:999x=1.0*exp(alf1*t).*sin(2*pi*f1*t)+1.0*exp(alf2*t).*sin(2*pi*f2*t)+1.0*exp(alf3*t).*sin(2*pi*f3 *t);x=x+10*u(1+i*N:i*N+N); x=x/max(x);x=x.*w;X=fft(x,N);X=abs(X);Y=Y+X;endY=Y/1000;subplot(212);plot(f(1:N/2),Y(1:N/2));title('正弦加白噪声后与w时域相乘的FFT 1000次平均') xlabel('Hz')24681012-1-0.500.51x(t) w(t)05101520253035404550102030x(t)频谱Hz024681012-0.50.51正弦加白噪声后与w 时域相乘05101520253035404550510正弦加白噪声后与w 时域相乘的FFTHz正弦加白噪声的FFT 1000次平均05101520253035404550Hz正弦加白噪声后与w时域相乘的FFT 1000次平均05101520253035404550Hz2、clear all;%三个正弦信号相加,分段函数,进行频谱分析% 产生三个正弦相加的函数;N=512;f0=10;fs=100.0;t=[0:N-1];x=1.0*sin(2*pi*f0*t/fs)+1.0*sin(2*pi*2*f0*t/fs)+1.0*sin(2*pi*3*f0*t/fs);subplot(211);plot(t(1:N),x(1:N));title('x(t)')%加窗w=1-1.93*cos(2*pi*t/N)+1.29*cos(4*pi*t/N)-0.388*cos(6*pi*t/N)+0.0322*cos(8*pi*t/N); %w=1.0-cos(2*pi*t/N);x=w.*x;%加窗等于时域点乘% 应用FFT 求频谱;f=0:fs/N:fs/N*(N-1);X=fft(x,N);%先点乘再进行傅里叶变换X=abs(X)/N;subplot(212);plot(f(1:N/2),X(1:N/2));title('x(t)加窗之后的傅里叶变换')xlabel('Hz')x(t)0100200300400500600x(t)加窗之后的傅里叶变换05101520253035404550Hz%分段函数-----------------------------------------------M=170;L=N-2*M;x(1:M)=1.0*sin(2*pi*f0*t(1:M)/fs);x(M+1:2*M)=1.0*sin(2*pi*2*f0*t(1:M)/fs);x(2*M+1:N)=1.0*sin(2*pi*3*f0*t(1:L)/fs);figure;subplot(211);plot(t(1:N),x(1:N));title('x(t)为分段函数')%w=1-1.93*cos(2*pi*t/M)+1.29*cos(4*pi*t/M)-0.388*cos(6*pi*t/M)+0.0322*cos(8*pi*t/M); %x(1:M)=w(1:M).*x(1:M);X=fft(x,N);X=abs(X)/N;subplot(212);plot(f(1:N/2),X(1:N/2));title('x(t)频谱分析')xlabel('Hz')x(t)为分段函数0100200300400500600x(t)频谱分析05101520253035404550Hz3、%采样长度不同对FFT的影响----------------------------------------------------clear all;% 观察数据长度N的变化对DTFT分辨率的影响f1=2;f2=2.02;f3=2.07;fs=10;w=2*pi/fs;N=256;%N=256x=sin(w*f1*(0:N-1))+sin(w*f2*(0:N-1))+sin(w*f3*(0:N-1));f=0:fs/N:fs/2-1/N;X=fft(x);X=abs(X);subplot(221);plot(f(45:60),X(45:60));title('N=256');grid on;subplot(223)plot(f(1:N/2),X(1:N/2));title('N=256');grid on;xlabel('Hz')%N=N*4;%N=1024x=sin(w*f1*(0:N-1))+sin(w*f2*(0:N-1))+sin(w*f3*(0:N-1));f=0:fs/N:fs/2-1/N;X=fft(x);X=abs(X);subplot(222)plot(f(45*4:4*60),X(4*45:4*60));title('N=1024');grid on;xlabel('Hz')1.52 2.5N=256N=256Hz1.52 2.5020*******N=1024Hz4、%补零的影响------------------------------------------------------------------- clear;% 计算长度为N 的原始信号的DTFT f1=2.67;f2=3.75;f3=6.75;fs=20;w=2*pi/fs; N=16;x=sin(w*f1*(0:N-1))+sin(w*f2*(0:N-1)+pi/2)+sin(w*f3*(0:N-1)); f=0:fs/N:fs/2-1/N; X=fft(x); X=abs(X);f=fs/N*(0:N/2-1); subplot(221)stem(f,X(1:N/2),'.');title('不补零');grid on; xlabel('Hz')% 在数据末补N 个零 x(N:2*N-1)=0;X=fft(x); X=abs(X); f=fs*(0:N-1)/(2*N); subplot(222)stem(f,X(1:N),'.');title('补N 个零');grid on; xlabel('Hz')% 在数据末补7*N 个零 x(N:8*N-1)=0;X=fft(x); X=abs(X); f=fs*(0:4*N-1)/(8*N); subplot(223)stem(f,X(1:4*N),'.');title('补7N 个零');grid on; xlabel('Hz')% 在数据末补29*N 个零 x(N:30*N-1)=0; X=fft(x); X=abs(X); f=fs*(0:15*N-1)/(30*N); subplot(224)plot(f,X(1:15*N));title('补29N 个零');grid on; xlabel('Hz')510不补零Hz 0510补N 个零Hz 0510补7N 个零Hz510510补29N 个零Hz5、%x(n)是两个正弦信号和一个白噪声相加,FFT 和IFFT ------------------ ----------- clear all;% 产生两个正弦加白噪声; N=256;f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs;x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N);% 应用FFT 求频谱;subplot(3,1,1);plot(x(1:N/4));title('x(n)')f=-0.5:1/N:0.5-1/N;X=fft(x);y=ifft(X);subplot(3,1,2);plot(f,fftshift(abs(X)));title('fft(x)')subplot(3,1,3);plot(real(x(1:N/4))); title('x(n)实部')x(n)010203040506070fft(x)-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5x(n)实部0102030405060706、%fftfilt和conv比较n=0:2;h=1./2.^n;for i=1:51x(i)=(i-1)/5;endfor i=52:100x(i)=20-(i-1)/5;endy=fftfilt(h,x); % y=x*hz=conv(h,x);title('x(t)')hold;plot(x);subplot(312)title('f ftfilt 叠接相加法')hold;plot(y);subplot(313)plot(z);axis([0,100,0,20]);title ('conv 卷积')x(t)0102030405060708090100fftfilt 叠接相加法0102030405060708090100conv 卷积01020304050607080901007、clear;%补零后对频谱的影响% 计算长度为N的原始信号的DFTf0=50;%信号频率fs=200;%采样频率N=16;%抽样点数w=2*pi/fs;%数字角频率x=sin(w*f0*(0:N-1));X=fft(x);X=abs(X);f=fs/N*(0:N/2-1);%N/2的数据stem(f,X(1:N/2),'.');title('未补零');grid on;xlabel('Hz')% 在数据末补N个零x(N:2*N-1)=0;X=fft(x); X=abs(X);f=fs*(0:N-1)/(2*N);subplot(212)stem(f,X(1:N),'.');title('补N个零');grid on;xlabel('Hz')未补零Hz补N个零Hz8、%抽样频率不同的影响% 计算长度为N的原始信号的DFT% fs=100f0=50;fs=100;w=2*pi/fs;N=16;x=sin(w*f0*(0:N-1));X=fft(x);X=abs(X);f=fs/N*(0:N/2-1);subplot(221)stem(f,X(1:N/2),'.');grid on;xlabel('Hz')% fs=150f0=50;fs=150;w=2*pi/fs; N=16;x=sin(w*f0*(0:N-1)); X=fft(x); X=abs(X);f=fs/N*(0:N/2-1); subplot(222)stem(f,X(1:N/2),'.');grid on; xlabel('Hz')% fs=200f0=50;fs=200;w=2*pi/fs; N=16;x=sin(w*f0*(0:N-1)); X=fft(x); X=abs(X);f=fs/N*(0:N/2-1); subplot(223)stem(f,X(1:N/2),'.');grid on; xlabel('Hz')0204060-14HzHzHz9、%周期延拓 clear all;M=128;N=1024;alf=-3.0;f0=10;fs=100.0;%u=rand(1,5000);u=u-mean(u);t=[0:1/fs:(M-1)/fs];x1=1.0*exp(alf*t).*sin(2*pi*f0*t);t=[0:1/fs:(N-1)/fs];for i=0:7x(M*i+1:M*i+M)=x1(1:M);endx=x-mean(x);subplot(211);plot(t,x);% 应用FFT 求频谱;f=0:fs/N:fs/N*(N-1);X=fft(x,N);X=abs(X);%X=20*log10(X);subplot(212);plot(f(1:N/2),X(1:N/2));% stem(f,X,'.');grid on;xlabel('Hz')024681012Hz10、%正弦加白噪声并FFT频谱分析% 产生4096点的两个正弦加白噪声;N=4096;f1=2.1;f2=2.2;fs=5;a1=5;a2=3;w=2*pi/fs;x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+10*randn(1,N);% 应用FFT 求频谱;f=fs/N*(0:N/2-1);X=fft(x,N);X=abs(X);subplot(211);stem(f,X(1:N/2),'.');title('x(n)');grid on;xlabel('Hz')f=-0.5:1/N:0.5-1/N;subplot(212);stem(f,fftshift(X),'.');title('FFT x(n)');grid on;x(n)HzFFT x(n)11、clear;%补零对频谱分析的影响f1=10.8;f2=11.75;f3=12.55;fs=40;w=2*pi/fs;N=64;x=sin(w*f1*(0:N-1))+sin(w*f2*(0:N-1)+pi/2)+sin(w*f3*(0:N-1)); X=fft(x);X=abs(X);f=fs/N*(0:N/2-1);subplot(221)stem(f,X(1:N/2),'.');title('补0个零');grid on;xlabel('Hz')% 在数据末补3N个零x(N:4*N-1)=0;X=fft(x); X=abs(X);f=fs*(0:2*N-1)/(4*N);subplot(222)stem(f,X(1:2*N),'.');title('补3N个零');grid on;xlabel('Hz')% 在数据末补7*N个零x(N:8*N-1)=0;X=fft(x); X=abs(X);f=fs*(0:4*N-1)/(8*N);subplot(223)stem(f,X(1:4*N),'.');title('补7N个零');grid on;xlabel('Hz')% 在数据末补15*N个零x(N:16*N-1)=0;X=fft(x); X=abs(X);f=fs*(0:8*N-1)/(16*N);subplot(224)plot(f,X(1:8*N));title('补15N个零');grid on;xlabel('Hz')补0个零Hz 05101520补3N 个零Hz 05101520补7N 个零Hz5101520010203040补15N 个零Hz12、%窗函数的影响 clear;%窗函数对频谱的影响f1=10.8;f2=11.75;f3=12.55; fs=40; w=2*pi/fs;N=1024;t=[0:N-1]; f=fs/N*(0:N/2-1); %矩形窗x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); X=fft(x);X=2*abs(X)/N; subplot(221)stem(f,X(1:N/2),'.');grid on; xlabel('Hz') title('矩形窗') %升余弦窗x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); wt=0.5-0.5*cos(2*pi*t/N); x=wt.*x; X=fft(x);X=2*2*abs(X)/N; subplot(222)stem(f,X(1:N/2),'.');grid on;xlabel('Hz') title('升余弦窗') %平顶窗x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t);wt=1-1.93*cos(2*pi*t/N)+1.29*cos(4*pi*t/N)-0.388*cos(6*pi*t/N)+0.0322*cos(8*pi*t/N); x=wt.*x; X=fft(x);X=2*abs(X)/N; subplot(223)stem(f,X(1:N/2),'.');grid on; xlabel('Hz') title('平顶窗') %改进升余弦窗x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); wt=0.54-0.46*cos(2*pi*t/N); x=wt.*x; X=fft(x);X=1.852*2*abs(X)/N; subplot(224)stem(f,X(1:N/2),'.');grid on; xlabel('Hz')title('改进升余弦窗')Hz 矩形窗Hz升余弦窗Hz平顶窗Hz改进升余弦窗13、%频谱分析,平均clear all;%产生正弦加白噪声信号,求频谱,平均% 产生两个正弦加白噪声;N=1024;f0=10;fs=100.0;u=rand(1,N);u=u-mean(u);t=[0:1/fs:(N-1)/fs];x=1.0*sin(2*pi*f0*t)+0.1*sin(2*pi*2*f0*t)+10.0*u;% 应用FFT 求频谱;f=0:fs/N:fs/N*(N-1);X=fft(x,N);X=abs(X);X=20*log10(X);subplot(211);plot(f(1:N/2),X(1:N/2));title('平均前')% stem(f,X,'.');grid on;xlabel('Hz')Y=zeros(1,N);u=rand(1,1000*N);u=u-mean(u);for i=0:999x=1.0*sin(2*pi*f0*t)+0.1*sin(2*pi*2*f0*t)+10*u(1+i*N:i*N+N);X=fft(x,N);X=abs(X);Y=Y+X;endY=Y/1000;Y=20*log10(Y);subplot(212);plot(f(1:N/2),Y(1:N/2));title('平均后')axis([0,50,0,60]);xlabel('Hz')5101520253035404550204060平均前Hz 05101520253035404550204060平均后Hz1高通 2带阻 3低通 4带通 5切比雪夫 带通1、1高通,双线性变换法,巴特沃斯(两种方法)% to design 高通high-pass DF with s=2/Ts[(z-1)/(z+1)]%-给出:通带下限频率,阻带上限频率,通带衰减,阻带衰减,采样频率----------------------------------------------- %---------------------- clear allwp=.8*pi;%通带下限频率 ws=.44*pi;%阻带上限频率 rp=3;%通带衰减 rs=20;%阻带衰减 Fs=2000;%采样频率%设计模拟滤波器% Firstly to finish frequency prewarping;wap=2*Fs*tan(wp/2);%通带截止频率,模拟角频率,预畸变 was=2*Fs*tan(ws/2);%阻带截止频率,模拟角频率,预畸变 [n,wn]=buttord(wap,was,rp,rs,'s');%求取模拟低通滤波器阶数,n 是模拟低通滤波器的阶次,巴特沃斯滤波器阶数的选择,(最小阶数的选择)[z,p,k]=buttap(n);%设计模拟低通滤波器,极点,零点,增益 [b,a]=zp2tf(z,p,k);%零-极点增益模型转换为传递函数模型 [bt,at]=lp2hp(b,a,wap);%模拟低通滤波器转换成高通滤波器,G(p)转换成H (s ),wap 为低通的截止频率 %模拟转换成数字% Note: z=(2/ts)(z-1)/(z+1);ts=1,that is 2fs=1,fs=0.5; [bz,az]=bilinear(bt,at,Fs);%实现双线性变换,由模拟滤波器H (s )得到数字滤波器H (Z ),bz,az 分别是H (Z )的分子分母多项式的系数 [h,w]=freqz(bz,az,256,1);%求离散系统频响特性,H 包含了离散系统频响在 0~pi 范围内N 个频率等分点的值,w 则包含了范围内N 个频率等分点。