MATLAB信號處理例題

1、共享知識分享快樂卑微如螻蟻、堅強似大象例1設(shè)方波的數(shù)學(xué)模型為:,基頻:4E11f (t)sin 0t Sin3 0tSin5 0t35用MATLA軟件完成該方波的合成設(shè)計 MATLAB源 程序t=-10:0.1:10; % 設(shè)定一個數(shù)組有201個點，方波周期為20e=5;w=pi/10; % 設(shè)定方波幅值為5, W代表w0m=-5*sign(t); % 給定幅值為5的方波函數(shù)y1=(-4*epi)*si n(w*t); % 計算 1 次諧波y3=(-4*epi)*(sin(w*t)+sin(3*w*t)3); %計算 3 次諧波y5=(-4*epi)*(sin(w*t)+sin(3*w*t)3

2、+sin(5*w*t)5); %計算 5 次諧波plot(t,y1,'y');hold; grid; %用黃色點線畫出1次諧波及網(wǎng)格線,并在同一張圖上畫其余曲線plot(t,y3,'g'); % plot(t,y5,'b'); % plot(t,m,'-k'); %用綠色點線畫出3次諧波 用藍色點線畫出5次諧波 用黑色實線畫方波title(' 方波合成');XIabeI('t')ylabel('f(t)'); %為圖形加上標(biāo)題n=50; %合成任意次方波，n決定方波的合成次數(shù)，在此給

3、定50yn=0; %設(shè)置初始值for i=1: nyn=yn+(-4*epi)*(1(2*i-1)*si n(2*i-1)*w*t);end; %計算n次諧波合成plot(t,yn,'r') %用紅色實線畫出n次諧波合成從圖中我們可以看到GibbS現(xiàn)象。在函數(shù)的間斷點附近，增加傅里葉級數(shù)的展開次數(shù)，雖然可以使其間斷點附近的微小振動的周期變小， 但振幅卻不能變小此現(xiàn)象在控制系統(tǒng)表現(xiàn)為：當(dāng)求控制系統(tǒng)對階躍函數(shù)的響應(yīng)時， 超調(diào)量總是存在 的。例2 (P110) MATLAB中函數(shù)FFT應(yīng)用舉例。%MATLAB中函數(shù)FFT應(yīng)用舉例t=0:0.001:0.6;X=Si n( 2*pi*5

4、0*t)+si n(2*pi*120*t);y=x+2*ra ndn( size(t)SUbPlot(2,1,1)plot(y(1:50)xlable(時間軸 t')ylable(信號值 f(t)')title(芷弦波 + 隨即噪聲','Fo ntSize ',10)y=fft(y,512);f=1000*(0:256)/512;subplot(2,1,2)plot(f,Y(1:257)set(gca,'Xtick ',0,50,100,150,200,250,300,350,400,450,500)set(gca,'XtickL

5、abel ',X50100150200250300350400450500') xlabel(頻率軸 omega')ylabel(頻譜幅值 F(omega)title(信號頻譜',XontSize)0)例3例3.8.3 有一二階系統(tǒng)，阻尼比 =0.47,固有頻率 n 500 HZ ,采樣間隔 t0.0004s采樣點數(shù)N = 256。試計算理論幅頻特性與由系統(tǒng)階躍響應(yīng)計算出的幅頻特性數(shù)據(jù)值，并畫出兩個計算結(jié)果的幅頻特性曲線。% example 3.8.3 MATLAB PRoGRAMN=256;dt=0.0004Wn=500;Seta=0.47;dw=2*pi(N

6、*dt);a=w n2;b=1,2*seta*w n,a;t=0:dt:(N-1)*dt;C=SteP(a,b,t);w=0:dw:(N-1)*dw;mag,phase=bode(a,b,w);ycw=fft(c);Re=real(ycw);Im=imag(ycw);for i=1:NRw(i,1)=1-Im(i,1)*(i-1)*dw*dt;Iw(i,1)=Re(i,1)*(i-1)*dw*dt;endffw=Rw+Iw*sqrt(-1);Aw=abs(ffw)SemilogX(W,20*log10(mag),'r-')axis(100,10000,-30,10)text(6

7、00,12,'幅頻特性')hold Onsemilogx(w,20*log10(Aw)axis(100,10000,-30,10)grid On幅頻特性例4例6.2.4用MATLAB中的函數(shù)XCoRR求出下列兩個周期信號的互相關(guān)函數(shù)，式中的 f=10Hz。x(t) Sin(2 ft), y(t) 0.5sin(2 ft 90 )%例6.2.4中計算兩個周期信號互相關(guān)函數(shù)的MATLAB程序N=500;Fs=500;n=0:N-1;t=n/Fs;Lag=200;X=S in (2*pi*10*t);y=0.5*si n(2*pi*10*t+90*pi180);c,lags=xcor

8、r(x,y丄ag,' Un biased');SUbPlot(2,1,1)plot(t,x,'r')hold Onplot(t,y,'b')XIabel('t')ylabel('x(t)y(t)');title('原周期信號')gridhold offsubplot(2,1,2);Plot(IagS/Fs,c,'k');XIabel('t')ylabel('Rxy(t)'); title('互相關(guān)函數(shù)');grid例5例6.2.5若有信

9、號為x(t) sin(2 fit) 2sin(2 f2t)(t)。設(shè)采樣頻256和式中，f150Hz, f2100Hz ,(t)為白噪聲(用 MATLAB中的函數(shù)產(chǎn)生)率FS 2000;試用周期法并應(yīng)用MATLAB編程計算，當(dāng)數(shù)據(jù)長度分別為N1N21024兩種情況下上述信號的功率譜。%例6.2.5中周期圖法計算信號功率譜的MATLAB程序CIfFs=2000;%情況1：數(shù)據(jù)長度N1=256N1=256;NIfft=256;n1=0:N1-1;t1= n1Fs;f1=50;f2=100;Xn 1=si n(2*pi*f1*t1)+si n( 2*pi*f2*t1)+ra ndn (1,N1);P

10、XX仁 10*log10(abs(fft(x n1,N1fft).2)N1);f1=(0:Ie ngth(Pxx1)-1)*Fsle ngth(Pxx1);subplot(2,1,1)plot(f1,Pxx1)ylabel('功率譜(dB)');title('數(shù)據(jù)長度N仁256')grid%情況2:數(shù)據(jù)長度N2=1024N2=1024;N2fft=1024;n2=0:N2-1;t2=n2Fs;f1=50;f2=100;Xn 2=si n(2*pi*f1*t2)+2*si n( 2*pi*f2*t2)+ra ndn (1,N2);Pxx2=10*log10(abs

11、(fft(x n2,N2fft).2)N2); f2=(0:Ie ngth(Pxx2)-1)*Fsle ngth(Pxx2); subplot(2,1,2)plot(f2,Pxx2)xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('數(shù)據(jù)長度 N2=1024')例6 (例3.8.1)分別用ConV和FFT算法計算序列：x(n)為在區(qū)間0 , 1上均勻分布的N 點隨機序列，表示為：x(n) = rand(1,N),h(n)為均值為零、方差為1的N點高斯分布隨機序列， 表示為：h(n) = rand(1,L),試求1

12、N 150時的卷積并比較其運算時間。9例3.8.1直接卷積和快速卷積的比較%Con v_time=ZeroS(1,150);fft_time=ZeroS(1,150);%for N=1:150tc=0;tf=0;%初始化L=2*N-1;%加長序列長度n u=rou nd(log10(L)log10(2)+0.45);L=2nu;%使點數(shù)成為2的幕次for 1=1:100h=ra ndn (1,N);%產(chǎn)生兩個隨機序列x=ran d(1,N);t0=clock;yc=c On v(h,x);%直接卷積計算t仁 etime(clock,t0);tc=tc+t1;%直接卷積運算的時間t0=clock

13、;y2=ifft(fft(h,L).*fft(x,L);%快速卷積計算t2=etime(clock,t0);tf=tf+t2;%快速卷積計算的時間end%con v_time(N)=tc/100;%直接卷積計算的平均時間fft_time(N)=tf/100;%快速卷積計算的平均時間end%n=1:150;SUbPlot(1,1,1);%圖形顯示上述兩種卷積的計算時間plot( n(25:150),co nv_time(25:150), n(25:150),fft_time(25:150)上述兩種卷積的計算時間的比較如圖3.8.5所示。圖385 兩種卷積計算時間比較例7 (例3.8.2 )運用F

14、FT求取矩形脈沖(t)1,0 t 10, t 1的譜，說明采樣頻率低引起的混疊現(xiàn)象。(1)先編寫有一定通用性的函數(shù)文件Cftbyfft.mcftbyfft.mfun ctio n AW,f=cftbyfft(wt,t,flag)%cftbyfft.m%本程序采用FFT計算連續(xù)時間FoUie變換。輸出幅頻數(shù)據(jù)對(f,AW)。%俞入量(wt,t)為已經(jīng)窗口化了的時間函數(shù)wt(t),它們分別是長度為N的向量。%寸于時限信號，應(yīng)使該取值時段與窗口長度相比足夠小。以提高頻率分辨率。%寸于非時限信號，窗口長度的選取應(yīng)使窗口外的函數(shù)值小%到可忽略，以提高近似精度。%俞入宗量flag控制輸出CFT的頻率范圍。

15、% flag 取非0時(缺省使用)，頻率范圍在0,fs);% flag 取0時，頻率范圍在-fs2,fs2)采樣點數(shù)，應(yīng)為2的幕次，以求快速。 窗口長度時間分辨率實施FFT變換算得0,fs 上的N點CFT值頻率分辨率if nargin=2;flag=1;e ndN=le ngth(t);%T=t(le ngth(t)-t(1);%dt=TN;%W0=fft(wt);%df=1T;%W=dt*W0;%n=0:1:(N-1)%e以上計算結(jié)果改寫到if flag=On=-N/2:(N/2-1);W=fftshift(W);%end f=n *df;%AW=abs(W);%if n argout=0p

16、lot(f,AW);grid,xlabel('-fs2,fs2范圍產(chǎn)生周期序列的頻譜頻率分量向量幅頻譜數(shù)據(jù)向量頻率 f)ylabel('w(f)')end(2)運行以下指令,繪制時域波形和幅頻譜M=5;%做2的幕次用。本例把 M設(shè)得較小，是為了觀察混疊。ten d=1;%波形取非零值的時間長度。T=10;%窗口化長度應(yīng)足夠大，以減小窗口化引起的泄露“旁瓣”效應(yīng)。N=2M;%采樣點數(shù)，取2的幕是為了使FFT運算較快。dt=TN;%以上T、N的取值應(yīng)使N/T = fs采樣頻率大于兩倍時間波形%帶寬，以克服采樣引起的頻譜混疊。%在本例中，據(jù)理論分析知W(f=7.5)=Sa (

17、7.5*pi)=1(7.5%*pi)<5%.因此，可近似認為本例時間信號帶寬為7.5H z.n=0:N-1;%采樣序列t=n *dt;%采樣點時間序列W=ZerOS(SiZe(t,2),1);Tow=fi nd(te nd-t)>O);%產(chǎn)生非零波形時段的相應(yīng)序列w(Tow,1)=o nes(le ngth(Tow),1);%在窗口時段內(nèi)定義的完整波形plot(t,w,'b','Li neWidth',2.5),title('Time WaVeform');XIabeI('t->')由上述程序畫出的被窗口化的時間

18、波形如圖3.8.12所示：圖3.8.12被窗口化的時間波形AW,f=cftbyfft(w,t,O);ff=f+eps;%為避免下面指令出現(xiàn)0/0而采取的措施AWW=abs(s in (pi*ff).(pi*ff);plot(f,AW,'b - ',ff,AWW,'r:');title('Aliasing CaUSed by Undersampling')XIabel('f->')ylabel('W(f)'),lege nd('by FFTTTheoreticaI')由上述程序畫出的“欠”采樣時

19、引起的混疊曲線如圖3.8.13所示：圖3813 “欠”采樣時引起的混疊例8 （例3.8.5 ）試分別求出含白噪聲干擾的正弦信號與白噪聲的自相關(guān)函數(shù)，并對它 們的結(jié)果進行比較。%example 3.8.5 MATLAB PRoGRAM N=1000;n=0:N-1;Fs=500;t=n/Fs;Lag=100;%x1=si n(2*pi*10*t)+0.6*ra ndn (1,le ngth(t);%c,lags=xcorr(x1,Lag,' Un biased')%subplot(2,2,1)%plot(t,x1)XIabel('t')ylabel('x1

20、(t)');title('含白噪聲的正弦信號x1');gridsubplot(2,2,2)%plot(lagsFs,c)XIabel('t')ylabel('Rxx1(t)');相關(guān)信號的最大延遲量 含白噪聲的正弦信號x1無偏自相關(guān)函數(shù)的計算 畫出x1曲線畫出x1自相關(guān)曲線title('x1 的自相關(guān)函數(shù)');gridx2=ra ndn (1,le ngth(t);%c,lags=xcorr(x2,Lag,' Un biased') %SUbPlot(2,2,3)%plot(t,x2)XIabel('

21、;t')ylabel('x2(t)');title(' 白噪聲 x2');gridsubplot(2,2,4)%plot(lagsFs,c)XIabel('t')ylabel('Rxx2(t)');title('x2 的自相關(guān)函數(shù)');發(fā)生白噪聲x2白噪聲x2的無偏自相關(guān)函數(shù)畫出x2曲線畫出x2自相關(guān)曲線grid由上述程序畫出的信號與自相關(guān)函數(shù)的曲線如圖3.8.19所示:含白噪聲的正弦信號Jd圖3.8.19信號與自相關(guān)函數(shù)由計算結(jié)果和和圖 3.8.19的曲線可以看出：同時具有周期性和白噪聲干擾的信號，其自相

22、關(guān)函數(shù)不但在T= 0時具有最大值，而且在T較大時自相關(guān)函數(shù)仍有明顯的周期性，它 的頻率和原周期信號的周期相同；而無周期的白噪聲，自相關(guān)函數(shù)在T= 0時也具有最大值，但當(dāng)T稍大時迅速衰減至零附近，利用自相關(guān)函數(shù)的這一特性可用來識別隨機信號中是否含 有周期成分和它的頻率。 例9（例4.4.3）試用MATLAB語言繪制巴特沃思低通模擬濾波器的平方幅頻相應(yīng)曲線，濾波器的階數(shù)分別為2, 5, 10, 20。解%例443的MATLAB 程序 n=0:0.01:2 for i=1:4SWitCh iCaSe 1N=2;CaSe 2N=5;CaSe 3N=10;CaSe 4N=20;%巴特沃思濾波器原型設(shè)計函

23、數(shù)，n:階數(shù);z,p,k:%濾波器零點，極點和增益%零極點增益轉(zhuǎn)換為傳遞函數(shù)%模擬濾波器頻率響應(yīng)函數(shù)%幅度平方函數(shù)end z,p,k=buttap(N);b,a=zp2tf( z, p,k); H,w=freqs(b,a, n); magH2=(abs(H).2; hold On plot(w,magH2);axis(0201);endxlabel( W/wc 'ylabel( H(jw) 2'' grid圖444 不同階次巴特沃思濾波器的幅度平方函數(shù) 例10（例446 ）試用MATLAB程序，確定一個模擬低通濾波器的階數(shù)N和截止頻率WC。設(shè)計指標(biāo)為：通帶邊界頻率200

24、 ,阻帶邊界頻率 300 ,通帶波紋為1dB,在 Z處，幅度衰減大于 18 dB。解%設(shè)計切比雪夫低通濾波器的MATLAB程序wp=200*pi;ws=300*pi;Rp=1;%通帶波紋Rs=16;%計算濾波器階數(shù)ebs=sqrt(10(Rp10)-1);A=10(Rs20);Wc=wp;Wr=WSwp;g=sqrt(A*A-1)ebs;%阻帶衰減%波紋系數(shù)%A為參變量%截止頻率%g為參變量N仁 log10(g+sqrt(g*g-1)log10(Wr+sqrt(Wr*Wr-1);%濾波器階數(shù)計算 N=CeiI(N1);%N 取整數(shù)運行上述程序后，可得濾波器的截止頻率WC和階數(shù)N為WC =628

25、.3185N =4 例11 (例444)繪制出階數(shù)分別為 2, 4, 6, 8的切比雪夫模擬低通濾波器的平方幅頻響應(yīng)曲線。解寫出其MATLAB程序如下：%繪制切比雪夫低通濾波器幅度平方函數(shù)的MATLAB程序n=0: 0.01: 2;for i=1:4SWitCh iCaSe 1N=2 ;CaSe 2N=4 ;CaSe 3N=6 ;CaSe 4N=8 ;endRp=1; z,p,k=cheblap(N,Rp);b,a=zp2tf(z,p,k); H,w=freqs(b,a, n); magH2=(abs(H).2;% OUtPUt%濾波器通帶波紋系數(shù)%設(shè)計切比雪夫模擬低通濾波器原型函數(shù),%乙p,

26、k分別為濾波器的零點、極點和增益%零點、極點和增益轉(zhuǎn)換為傳遞函數(shù)%模擬濾波器的頻率響應(yīng)%幅度平方函數(shù)%定義字符串變量posplot= 1'0 'u m2str(i); subplot(posplot) plot(w,magH2, k'); ylabel( H(jw)2 ' title( N= 'u m2str(N);gridend其幅度平方函數(shù)曲線如圖446所示。由圖446可以看出：切比雪夫濾波器在通帶內(nèi)具有等波紋起伏特性，而在阻帶內(nèi)則單 調(diào)下降。階次越高，特性越接近矩形。 例12(例447)試設(shè)計一巴特沃思模擬帶通濾波器，設(shè)計要求為：通帶頻率2kHz3

27、kHz,兩邊的過渡帶寬為0.5kHz ,若通帶波紋IdB ,阻帶衰減大于100dB.解%設(shè)計巴特沃思帶通濾波器MA TLAB程序%濾波器設(shè)計指標(biāo)wp=2000 3000*2*pi;ws=1500 3500*2*pi;Rp=1;RS=100;%計算階數(shù)和截止頻率N,Wn=buttord (wp,ws,Rp,Rs, ' s');FC=Wn/(2*pi);%計算濾波器傳遞函數(shù)多項式系數(shù)b,a=butter(N,Wn, ' S ' );%畫出濾波器幅頻特性w=li nspace(1,4000,1000)*2*pi;% 生成線性等間隔向量H=freqs(b,a,w);ma

28、gH=abs(H);PhaH=Un wrap(a ngle(H);Plot(W/(2*pi),20*log10(magH),' k ' );xlabel( '頻率(Hz)');ylabel( '幅度(dB)' );title( '巴特沃思模擬濾波器')grid On由上述程序畫出的巴特沃思模擬濾波器的幅頻特性如圖4411所示。圖4411巴特沃思帶通濾波器的幅頻特性 例13 (例448 )試設(shè)計一個切比雪夫模擬帶阻濾波器。要求的指標(biāo)為：阻帶上下邊界頻率為2KHz與3KHz ,兩側(cè)過渡帶為 0.5KHZ ,通帶波紋1dB,阻帶衰減大于

29、 60 dB。 解%設(shè)計切比雪夫帶阻濾波器MATLAB程序%濾波器設(shè)計指標(biāo)wp=2000 3000*2*pi;ws=1500 3500*2*pi;Rp=1;Rs=60;%計算階數(shù)和截止頻率N,Wn=cheblord (wp,ws,Rp,Rs, ' s');%計算濾波器傳遞函數(shù)多項式系數(shù)b,a=chebyl(N,Rp,Wn, ' stop ' , ' S ' );%畫出濾波器幅頻特性W=Ii nspace(1,4000,1000)*2*pi;% 生成線性等間隔向量H=freqs(b,a,w);% 相位展開magH=abs(H);PhaH=Un wr

30、ap(a ngle(H);Plot(W/(2*pi),20*log10(magH),' k ' );xlabel( '頻率(Hz)');ylabel( '幅度(dB)' );title( '切比雪夫模擬濾波器' )grid On由上述程序畫出的切比雪夫模擬濾波器的幅頻特性如圖4412所示。圖4.4.12切比雪夫帶通濾波器的幅頻特性 例14 （例5.3.4）一低通濾波器，其通帶和阻帶的技術(shù)指標(biāo)分別是0.2 IRP 0.25 dB0.3 IRZ 50dB試用頻率抽樣法設(shè)計一FIRN=20 。具有線性相位的濾波器，取jkjk由 H(k)

31、 Hd(e N ) Hd(e 20 ) , k 0, 1 ,19當(dāng)k=2時,正好是在通帶邊界頻率處有一個頻率抽樣點，即P 0.2下一個抽樣點為k=3 ,是阻帶上邊界頻率Z ，阻帶域通帶間無過渡帶，即Z 03則在通帶0P內(nèi)有三個抽樣點，阻帶上共有七個抽樣點。從而有Hk 1,1,1,0,由N=20 ,相位常數(shù)(N 1)/2(201）/29.5,其相位可表示為9k0.95 k ,0.95 (20 k),0 k 910 k 19h(n)，并有頻響內(nèi)插公式最根據(jù)式(5.3.37)得到H(k),利用離散傅里葉變換，求得程序如下:后可得FlR濾波器的頻響H(ej )。其MATLAB%例5.3.4中設(shè)計FIR

32、濾波器MATLAB程序N=20； alpha (N 1)/2;10: N 1 ; Wl (2*piN)*1HrS 1,1,1,ZerOSIl ,15),1,1;Hdr 1,1,0,0;Wdl 0,0.25,0.25,1;kl 0:floor(N 1)/2); k2 floor(N 1)/2) 1: N 1 ；angH = -alpha * ( 2 * Pi ) / N * kl ,alpha * (2 * Pi ) / N * ( N-k2 );H = Hrs. * exp (j * an gH);=real (ifft (H , N );db , mag , Pha ,grd , W = f

33、reqz -m ( h , 1 );subplot (1 , 1 , 1 )subplot ( 2,2 , 1 ) ; plot (wl( 1 : 11 ),o', wdl , Hdr );axis ( 0 , 1 , - 0.1 , 1.1 ); title ('頻率樣本： N = 20')xlabel (頻率(單位：Pi) ');ylabel ( Hr (k) JSet (gca , 'XtiCkMOde ' ,' manual',' XTiCk ',0,0.2,0.3,1 );Set (gca , YtiCk

34、MOde ' ,manual', Ytick',0,1); girdsubplot(2,2,2) ; Stem(I , h) ; axis(-1 , N ,-0.1 , 0.3)title (單位抽樣響應(yīng)')；ylabel( h(n)'); text (N+1 , -0.1 , n ')subplot (2,2,3 ) ; plot (ww/pi , Hr , wl (1:11)/pi , HrS (1:11)/pi , Hrs(1:11) ,6');axis ( 0 , 1 , -0.2,1.2 ) ; title (振幅響應(yīng)')xlabel (頻率(單位:Pi ) '); ylable ( Hr(W)')Set (gca , XTiCkMOde ', XtiCk ', 0,0.2,0.3 ,