lab01_Matlab平臺(tái)的使用_記得帶耳機(jī)_2014

1、Matlab平臺(tái)的使用1.總論Matlab是信號(hào)處理中功能強(qiáng)大且廣泛使用的工具。本實(shí)驗(yàn)的主要目標(biāo)是使你們熟悉使用Matlab。我們假定你們有過(guò)其他某種語(yǔ)言的編程經(jīng)驗(yàn)，應(yīng)該已經(jīng)了解變量，循環(huán)和函數(shù)等基本概念。這里介紹一下Matlab的基本編程結(jié)構(gòu)。本實(shí)驗(yàn)的三個(gè)目標(biāo)：1.學(xué)習(xí)基本的Matlab命令和語(yǔ)法，包括“幫助”系統(tǒng)2.學(xué)習(xí)編寫(xiě)你自己的Matlab的腳本文件，并像命令一樣運(yùn)行它們3.學(xué)習(xí)一點(diǎn)Matlab的高級(jí)編程技術(shù)，例如向量化2. 初識(shí)Matlab1. 啟動(dòng)Matlab，運(yùn)行Matlab的幫助。在命令窗口輸入“doc”，或在菜單中選“help”。2. 在命令窗口中，試試：helpwin he

2、lpwin plot helpwin colon %<- a VERY IMPORTANT notation helpwin ops helpwin zeros helpwin ones lookfor filter %<- keyword search 3.運(yùn)行Matlab的demo：在命令窗口輸入“demo”4.把Matlab當(dāng)計(jì)算器用，試試：pi*pi - 10 sin(pi/4) ans 2 %<- "ans" holds the last result 5.在Matlab中給變量賦值，試試：xx = sin( pi/5 ); cos( pi/5

3、) %<- assigned to what? yy = sqrt( 1 - xx*xx ) ans 注意觀察語(yǔ)句末尾有分號(hào)和無(wú)分號(hào)有什么差別。6.Matlab中的復(fù)數(shù)運(yùn)算的基本操作，試試：zz = 3 + 4i, ww = -3 + 4j real(zz), imag(zz) abs(zz,ww) %<- Vector constructor conj(zz+ww) angle(zz) exp( j*pi ) exp(j* pi/4, 0, -pi/4 ) 3. 再識(shí)Matlab3.1Matlab數(shù)組索引(a).冒號(hào)的含義試試如下命令的結(jié)果：jkl = 0 : 6 jkl = 2

4、 : 4 : 17 jkl = 99 : -1 : 88 ttt = 2 : (1/9) : 4 tpi = pi * 0:0.1:2 ; %注意這里有分號(hào)(b).從向量提取數(shù)字和或插入數(shù)字，思考如下xx的定義：xx = zeros(1,3), linspace(0,1,5), ones(1,4) xx(4:6) size(xx) length(xx) xx(2:2:length(xx) xx(2:2:end) 解釋以上代碼最后4行的結(jié)果(c).觀察如下賦值的結(jié)果：yy = xx; yy(4:6) = pi*(1:3) 現(xiàn)在寫(xiě)一條語(yǔ)句,xx用(b)中定義的方法，把xx的偶數(shù)索引的元素(即xx(

5、2),xx(4)等)的值替換為常數(shù)pp。使用向量替換，不要用循環(huán)。3.2Matlab腳本文件(a)看看關(guān)于向量的實(shí)驗(yàn)。把向量看作一個(gè)數(shù)字集合，試試：xk = cos( pi*(0:11)/4 ) %<-comment: compute cosines 解釋一下余弦函數(shù)的值是怎樣存儲(chǔ)在向量xk中的。xk(1)是什么意思？有定義xk(0)嗎？注：語(yǔ)句后沒(méi)有分號(hào)，所以語(yǔ)句會(huì)把xk的計(jì)算結(jié)果依次顯示在命令窗口，在百分號(hào)后面的文字是注釋?zhuān)梢院雎浴＃╞）Matlab中可以寫(xiě)循環(huán)，但是這不是最有效率的解決辦法。盡量避免循環(huán)，并使用冒號(hào)會(huì)提高效率。下面的代碼是用一個(gè)循環(huán)來(lái)計(jì)算余弦函數(shù)的值的： yy =

6、 ; %<- initialize the yy vector to be empty for k=-5:5 yy(k+6) = cos( k*pi/3 ) end yy 解釋為什么需要寫(xiě)yy(k+6)。如果用yy(k)會(huì)發(fā)生什么？不使用循環(huán)，而用3.1的方法，利用冒號(hào)，重新實(shí)現(xiàn)這個(gè)功能。（c）Matlab中實(shí)數(shù)和復(fù)數(shù)的繪制都很容易?；纠L圖命令繪制向量yy和向量xx的關(guān)系，試試：xx = -3 -1 0 1 3; yy = xx.*xx - 3*xx; plot( xx, yy ) zz = xx + yy*sqrt(-1) plot( zz ) %<- complex valu

7、es: plot imag vs. real 在命令窗口輸入“helpwin arith”，了解當(dāng)xx是一個(gè)向量時(shí)，操作xx.*xx是如何工作的。當(dāng)xx是一個(gè)矩陣時(shí)，操作xx.*xx又是如何工作的。當(dāng)你沒(méi)把握某個(gè)命令時(shí)，可隨時(shí)使用幫助系統(tǒng)。（d）使用Matlab的編譯器創(chuàng)建腳本文件，命名為mylab1.m，文件內(nèi)容如下：clear all; close all;tt = -1 : 0.01 : 1; xx = cos(5*pi*tt); zz = 1.4*exp(j*pi/2)*exp(j*5*pi*tt); plot(tt, xx, 'b-', tt, real(zz), &

8、#39;r-') %<- plot a sinusoid grid on title('TEST PLOT of a SINUSOID') xlabel('TIME (sec)') 解釋為什么real(zz)的圖形是一個(gè)正弦曲線。它的相位和幅度是多少？由所繪圖形的時(shí)移量計(jì)算相位。（e）在命令窗口輸入“mylab1”即可運(yùn)行上述m文件，即mylab1 %<-will run the commands in the file type mylab1 %<-will type out the contents of % mylab1.m to

9、 the screen 3.3 Matlab與聲音相關(guān)的函數(shù)命令這一節(jié)練習(xí)是關(guān)于聲音信號(hào)的，所以需要帶耳機(jī)來(lái)機(jī)房聽(tīng)效果。（a）運(yùn)行Matlab的聲音demo：在命令窗口輸入“xpsound”。你可以聽(tīng)到多種聲音，并可以查看對(duì)應(yīng)的時(shí)間序列，功率譜密度和聲譜圖（聲譜圖詳見(jiàn)7.3）。（b）現(xiàn)在在Matlab中生成一個(gè)音調(diào)（例如正弦），用soundsc()命令 soundsc(xx,fs)函數(shù)需用兩個(gè)參數(shù)，,xx是播放文件的數(shù)據(jù)的向量，fs是播放樣本的采樣率。另外，soundsc(xx,fs)會(huì)自動(dòng)把xx的值尺度變化成1到1之間，然后調(diào)用sound(xx,fs)播放信號(hào)。聽(tīng)一聽(tīng)在3.2（d）的第2行和

10、第3行創(chuàng)建了一個(gè)向量xx為2.5Hz的正弦序列。而你們要生成的正弦信號(hào)的頻率為2000Hz，時(shí)長(zhǎng)為0.9秒，采樣率fs為11025個(gè)樣本秒。采樣率說(shuō)明了采樣點(diǎn)之間的時(shí)間間隔，所以時(shí)間向量應(yīng)定義為：tt = 0:(1/fs):dur; 其中fs為采樣率，dur為時(shí)長(zhǎng)（以秒為單位）。可以在幫助系統(tǒng)中查找sound()和soundsc()了解它們的用法和含義。那么，你們的tt向量的長(zhǎng)度（即樣本數(shù)）是多少？4.用Matlab處理正弦信號(hào)在實(shí)驗(yàn)報(bào)告中要包括本節(jié)的圖形和實(shí)驗(yàn)結(jié)果描述。寫(xiě)一個(gè)Matlab的m文件來(lái)完成步驟（a）到步驟（d）。你們的實(shí)驗(yàn)報(bào)告也要包括這個(gè)m文件的代碼。（a）生成一個(gè)時(shí)間向量tt，

11、需要覆蓋下面的（b）中頻率為4000Hz的正弦信號(hào)的兩個(gè)周期。tt的定義方式使用3.2（d）的方式。如果我們使用T表示正弦的周期，定義向量tt的開(kāi)始時(shí)間為T(mén)，結(jié)束時(shí)間為T(mén)，那兩個(gè)周期就會(huì)包括tt=0。最后，確保你們正弦波的每個(gè)周期有至少25個(gè)樣本。也就是說(shuō)，當(dāng)你們使用冒號(hào)操作符定義時(shí)間變量時(shí)，要使增量足夠小，才能保證每個(gè)周期可以產(chǎn)生至少25個(gè)樣本。（b）生成以下兩個(gè)4000Hz的正弦序列：x1(t) = A1 cos(2(4000)(t  tm1)  x2(t) = A2 cos(2(4000)

12、(t  tm2) 其中A1是你的年齡，A2=1.2A1，tm1 = (37.2/M)T ，tm2 = (41.3/D)T ，這里M和D分別是你生日的月和日，T是周期繪制以上兩個(gè)信號(hào)在-TtT的圖形。為了輸出下面（d）的要求，使用subplot(3,1,1)和subplot(3,1,2)來(lái)使圖在一個(gè)圖形窗口內(nèi)，詳見(jiàn)helpwin subplot一般用法是：subplot(2,1,2);plot(xx);（c）創(chuàng)建第3個(gè)正弦信號(hào)為 x3(t) = x1(t) +

13、0;x2(t). 在Matlab中，這表示把兩個(gè)正弦向量中的值對(duì)應(yīng)相加。繪制 x3(t) 在-TtT的圖形，使用subplot(3,1,3)繪圖。（d）對(duì)每個(gè)圖都添加一個(gè)圖名，圖名都要包含你的姓名，用title函數(shù)詳見(jiàn) helpwin title, helpwin subplot4.1 復(fù)數(shù)幅度用一條語(yǔ)句生成以上正弦信號(hào)x1(t) 的值，可以使用復(fù)數(shù)幅度表示方法：x1(t) = ReXejwt  其中X和w為某個(gè)常數(shù)。5.復(fù)指數(shù)簡(jiǎn)介本節(jié)目標(biāo)是讓大家熟悉復(fù)數(shù)，并了解怎么使用它們表示正弦信號(hào)(例如 x(t) 

14、= Acos(t + ) )為復(fù)指數(shù)形式（例如z(t) = Aejfejwt），方法如下所示：x(t) = Acos(t + ) = Re Aejfejwt  5.2 Matlab中的復(fù)數(shù)Matlab可用于計(jì)算復(fù)值式子，并可以向量顯示結(jié)果，以下是一些Matlab的復(fù)數(shù)操作符：conj復(fù)共軛abs幅度 angle以弧度為單位的角度或相位 real 實(shí)部imag虛部 i,j 已預(yù)定義為x = 3 + 4i &#

15、160;exp(j*theta)  復(fù)指數(shù)函數(shù) ej 每個(gè)函數(shù)都以向量或矩陣作為輸入?yún)⒘?。注意：Matlab中沒(méi)有 mag()和phase() 函數(shù)。5.3 向量化Matlab最重要的部分就是矩陣向量語(yǔ)法。大多數(shù)情況下，循環(huán)都可以用向量操作代替，因?yàn)樾稳鏴xp() 和cos() 等函數(shù)都被定義為向量輸入，例如cos(vv) = cos(vv(1), cos(vv(2), cos(vv(3), . cos(vv(N) 其中vv是N個(gè)元素的行向量。向量化可用于簡(jiǎn)化代碼。如果有如下代碼來(lái)繪制某個(gè)信號(hào)：M = 200; for k=

16、1:M x(k) = k; y(k) = cos( 0.001*pi*x(k)*x(k) ); end plot( x, y, 'ro-' ) 你可以把循環(huán)替換掉，而只用如下3行代碼即可：M = 200; y = cos( 0.001*pi*(1:M).*(1:M) ); plot( 1:M, y, 'ro-' ) 用這種向量化的方法編寫(xiě)23行代碼完成以下Matlab代碼，不使用循環(huán)。（注：當(dāng)xx是向量時(shí)，xx*xx和xx.*xx是不同的）%- make a plot of a weird signal N = 200; for k=1:N xk(k) = k/

17、50; rk(k) = sqrt( xk(k)*xk(k) + 2.25 ); sig(k) = exp(j*2*pi*rk(k); end plot( xk, real(sig), 'mo-' ) 5.4函數(shù)函數(shù)是一種特別的M-文件，可以接受輸入（矩陣或向量），也可以得到輸出。關(guān)鍵詞function必須是文件的第一個(gè)詞，并且定義了函數(shù)，M-文件的第一行定義函數(shù)如何傳遞輸入和輸出參量。文件擴(kuò)展名必須使用小寫(xiě)的“m”，例如my_func.m。如下函數(shù)有一些小錯(cuò)誤，請(qǐng)找找錯(cuò)在哪里：（至少有3處錯(cuò)誤）matlab mfile xx,tt = badcos(ff,dur) %BADCO

18、S Function to generate a cosine wave % usage: % xx = badcos(ff,dur) % ff = desired frequency in Hz % dur = duration of the waveform in seconds % tt = 0:1/(100*ff):dur; %- gives 100 samples per period badcos = cos(2*pi*freeq*tt); 錯(cuò)誤有：第一個(gè)詞必須是“function”。而且“freeq”在使用前沒(méi)有定義。最后，第一行指明要“xx”作輸出，因此“xx”需要在函數(shù)體中出

19、現(xiàn)在至少一個(gè)賦值行的左邊。正確的函數(shù)寫(xiě)法應(yīng)該是：function xx,tt = goodcos(ff,dur) tt = 0:1/(100*ff):dur; %- gives 100 samples per period xx = cos(2*pi*ff*tt); 可在命令窗口輸入函數(shù)的文件名（注意不是函數(shù)名）來(lái)調(diào)用它（也可以在另一個(gè)m文件調(diào)用它），例如，如果以上goodcos函數(shù)的文件名為testcos.m，且ff=1，dur=10，則調(diào)用時(shí)應(yīng)輸入testcos(1,10)。你可以嘗試一下，如果你輸入goodcos(1,10),matlab會(huì)有什么反應(yīng)。6.復(fù)指數(shù)6.1 生成正弦信號(hào)的M-

20、文件寫(xiě)一個(gè)可以生成單一正弦信號(hào)x(t) = Acos(t+)的函數(shù)，使用4個(gè)輸入?yún)⒘浚悍華，頻率w，相位f和時(shí)長(zhǎng)dur。函數(shù)應(yīng)當(dāng)返回兩個(gè)輸出參量：正弦信號(hào)的值x和對(duì)應(yīng)的時(shí)間t。確保函數(shù)生成的正弦信號(hào)在每個(gè)周期有20個(gè)值，函數(shù)名為one_cos()。提示：可借鑒上面的goodcos()函數(shù)。繪制你們的one_cos()函數(shù)，參數(shù)選為：A=95，w=200p弧度秒，f=p/5弧度，時(shí)長(zhǎng)為0.025秒。推導(dǎo)所繪圖形的周期和相位是否正確。如果周期以毫秒為單位是多少？7線性調(diào)頻脈沖chirp線性調(diào)頻脈沖chirp信號(hào)（或線性掃頻信號(hào)）是一個(gè)正弦信號(hào)，其頻率從一個(gè)初值線性變化到終值。

21、Chirp信號(hào)的表達(dá)式可寫(xiě)為：x(t)=cos(t)， 其中 (t) = 2  t2+ 2 f0 t +  對(duì)(t) 求導(dǎo)后得到隨時(shí)間線性變化的瞬時(shí)頻率：fi (t) = 2t + f0 fi (t) 的斜率等于 2，截距為f0。如果信號(hào)于t=0開(kāi)始，f0也是初始頻率。這種時(shí)變角度造成的頻率變化稱(chēng)為頻率調(diào)制。這類(lèi)信號(hào)是頻率調(diào)制FM信號(hào)的一個(gè)例子。更一般的，我們通?？紤]它們?yōu)楦笠活?lèi)的角度調(diào)制信號(hào)。

22、最后，因?yàn)轭l率的線性變化可以生成一個(gè)類(lèi)似汽笛聲或鳥(niǎo)叫聲（chirp），所以線性FM信號(hào)也可以稱(chēng)為chirp信號(hào)。7.1chirp信號(hào)的Matlab合成方法以下的Matlab代碼可以合成一個(gè)chirp信號(hào)fsamp = 11025; dt = 1/fsamp; dur = 1.8; tt = 0 : dt : dur; psi = 2*pi*(100 + 200*tt + 500*tt.*tt); xx = real( 7.7*exp(j*psi) ); soundsc( xx, fsamp ); （a）確定合成信號(hào)的總時(shí)長(zhǎng)（秒），確定tt向量的長(zhǎng)度（樣本數(shù)）（b）在Matlab中，只能合成離散

23、時(shí)間信號(hào)，所以對(duì)于chirp信號(hào)，我們處理為：x(tn)=Acos( 2  tn2+ 2 f0 tn+ )其中，tn=nTs表示離散時(shí)間值。在以上的Matlab代碼中，tn的值是多少？A, m, f0和f的值是多少？（c）確定以上Matlab代碼合成的頻率的范圍（Hz為單位）。手繪瞬時(shí)頻率隨時(shí)間的變化情況。聽(tīng)到的最小頻率和最大頻率是多少？（d）聽(tīng)聽(tīng)信號(hào)的頻率是上升還是下降（使用soundsc()）。注意soundsc()需要知道信號(hào)創(chuàng)建時(shí)的采樣率。詳見(jiàn)幫助。7.2chirp的函數(shù)使用以下代碼來(lái)寫(xiě)一個(gè)可以合成chirp信號(hào)的函

24、數(shù)，根據(jù)注釋的提示編寫(xiě)。function xx,tt = mychirp( f1, f2, dur, fsamp ) %MYCHIRP generate a linear-FM chirp signal % % usage: xx = mychirp( f1, f2, dur, fsamp ) % % f1 = starting frequency % f2 = ending frequency % dur = total time duration % fsamp = sampling frequency (OPTIONAL: default is 11025) % % xx = (vector of) samples of the chirp signal % tt = vector of time instants for t=0 to t=dur % if( nargin < 4 ) %- Allow optional input argument fsamp = 11025; end 請(qǐng)生成一個(gè)chirp，其頻率起始于2500Hz，終止于500Hz