【文章內(nèi)容簡(jiǎn)介】 1,sdate)。 semilogy(cdate,logpred1,39。–39。,cdate,pop,39。+39。)。 grid on logres2 = log10(pop) –polyval(logp2,sdate)。 plot(cdate,logres2,39。+39。) 上面的圖不令人滿意，下面，我們用二階的對(duì)數(shù)分析： logp2 = polyfit(sdate,log10(pop),2)。 logpred2 = 10.^polyval(logp2,sdate)。 semilogy(cdate,logpred2,39。–39。,cdate,pop,39。+39。)。 grid on r = pop – 10.^(polyval(logp2,sdate))。 plot(cdate,r,39。+39。) 這種余量分析比多項(xiàng)式擬合的余量分析圖案要隨機(jī)的多（沒有很強(qiáng)的規(guī)律性），可以預(yù)見，隨著人數(shù)的增加，余糧所反映的不確定度也在增加，但總的說來，這種擬合方式要強(qiáng)好多！ 誤差邊界 誤差邊界常用來反映你所用的擬合方式是否適用于數(shù)據(jù)，為得到誤差邊界，只需在 polyfit()中傳遞第二個(gè)參數(shù)，并將其送入 polyval(). 下面是一個(gè)二階多項(xiàng)式擬合模型 ,年份已被標(biāo)準(zhǔn)化 ,下面的代碼用了 2σ ,對(duì)應(yīng)于 95%的可置信度 : [p2,S2] = polyfit(sdate,pop,2)。 [pop2,del2] = polyval(p2,sdate,S2)。 plot(cdate,pop,39。+39。,cdate,pop2,39。g–39。,cdate,pop2+2*del2,39。r:39。,... cdate,pop2–2*del2,39。r:39。), grid on 差分方程和濾波 MatLab 中的差分和濾波基本都是對(duì)向量而言的，向量則是存儲(chǔ)取樣信號(hào)或序列的． 函數(shù) y = filter(b, a, x)將用 a,b描述的濾波器處理向量 x,然后將其存儲(chǔ)在向量 y中 , filter()函數(shù)可看為是一差分方程 a1y(n)=b1*x(1)+b2*x(2)+...a2*y(2)... 如有以下差分方程 :y(n)=1/4*x(n)+1/4*x(n1)+1/4*x(n2)+1/4*x(n3),則 a = 1。 b = [1/4 1/4 1/4 1/4]。 我們載入數(shù)據(jù) ,取其第一列 ,并計(jì)算有 : load x = count(:,1)。 y = filter(b,a,x)。 t = 1:length(x)。 plot(t,x,39。–.39。,t,y,39。–39。), grid on legend(39。Original Data39。,39。Smoothed Data39。,2) 實(shí)現(xiàn)所表示的就是濾波后的數(shù)據(jù) ,它代表了 4 小時(shí)的平均車流量 MatLab 的信號(hào)處理工具箱中提供了很多用來濾波的函數(shù) ,可用來處理實(shí)際問 題 ! 快速傅立葉變換 (FFT) 傅立葉變換能把信號(hào)按正弦展開成不同的頻率值 ,對(duì)于取樣信號(hào) ,用的是離散傅立葉變換 . FFT 是離散傅立葉變換的一種高速算法 ,在信號(hào)和圖像處理中有極大的用處 ! fft 離散傅立葉變換 fft2 二維離散傅立葉變換 fftn n 維離散傅立葉變換 ifft 離散傅立葉反變換 ifft2 二維離散傅立葉反變換 ifftn n 維離散傅立葉反變換 abs 幅度 angle 相角 unwrap 相位按弧度展開 ,大于π的變換為 2π的補(bǔ)角 fftshift 把零隊(duì)列移至功率譜中央 cplxpair 把數(shù)據(jù)排成復(fù)數(shù)對(duì) nextpow2 下兩個(gè)更高的功率 向量 x 的 FFT 可以這樣求 : x = [4 3 7 –9 1 0 0 0]’ y = fft(x) y = – – – – + – – – + + x雖然是實(shí)數(shù) ,但 y是復(fù)數(shù) ,其中 ,第一個(gè)是因?yàn)樗浅?shù)相加的結(jié)果 ,第五個(gè)則對(duì)應(yīng)于奈奎斯特頻率 ,后三個(gè)數(shù)是由于負(fù)頻率的影響 ,它們是前面三個(gè)數(shù)的共軛值 ! 下面 ,讓我們來驗(yàn)證一下太陽(yáng)黑子活動(dòng)周期是 11 年 !Wolfer數(shù)記錄了 300年太陽(yáng)黑子的數(shù)量及大小 : load year = sunspot(:,1)。 wolfer = sunspot(:,2)。 plot(year,wolfer) title(’Sunspot Data’) 現(xiàn)在來看看其 FFT: Y = fft(wolfer)。 Y 的幅度是功率譜 ,畫出功率譜和頻率的對(duì)應(yīng)關(guān)系就得出了周期圖 ,去掉第一點(diǎn) ,因?yàn)樗皇撬袛?shù)據(jù)的和 ,畫圖有 : N = length(Y)。 Y(1) = []。 power = abs(Y(1:N/2)).^2。 nyquist = 1/2。 freq = (1:N/2)/(N/2)*nyquist。 plot(freq,power), grid on xlabel(’cycles/year’) title(’Periodogram’) 上面的圖看起來不大方便 ,下面我們畫出頻譜 周期圖 period = 1./freq。 plot(period,power), axis([0 40 0 2e7]), grid on ylabel(’Power’) xlabel(’Period(Years/Cycle)’) 為了得出精確一點(diǎn)的解 ,如下 : [mp index] = max(power)。 period(index) ans = 變換后的幅度和相位 abs()和 angle()是用來計(jì)算幅度和相位的 先創(chuàng)建一信號(hào) ,再進(jìn)行分析 ,unwarp()把相位大于π的變換為 2π的補(bǔ)角 : t = 0:1/99:1。 x = sin(2*pi*15*t) + sin(2*pi*40*t)。 y = fft(x)。 m = abs(y)。 p = unwrap(angle(y))。 f = (0:length(y)–1)39。*99/length(y)。 subplot(2,1,1), plot(f,m), ylabel(39。Abs. Magnitude39。), grid on subplot(2,1,2), plot(f,p*180/pi) ylabel(39。Phase [Degrees]39。), grid on xlabel(39。Frequency [Hertz]39。) 可以發(fā)現(xiàn)幅度曲線關(guān)于奈奎斯特頻率對(duì)稱 ,只有 050Hz的信息是有用的 ! FFT 的長(zhǎng)度與速度 可以為 FFT 加上第二個(gè)參數(shù) ,告訴 MatLab這是 n點(diǎn) y = fft(x,n),若 x長(zhǎng)度大于 n,軟件自動(dòng)補(bǔ) 0,否則截取 x. 若 : 1. n 為 2 的冪 ,軟件將執(zhí)行基 2 快速傅立葉算 法 ,這時(shí)的運(yùn)算速度是最快的 2. n 為合數(shù) ,軟件將 n 分解為素?cái)?shù)來算 ,計(jì)算量與 n 的值有關(guān) .n 為 1013將比 1000點(diǎn)的速度慢的多 ! 3. n 為素?cái)?shù) ,軟件執(zhí)行 DFT 的公式 ,此時(shí)最慢 矩陣和線性代數(shù) MatLab 中的矩陣 MatLab 中有好多函數(shù)可以產(chǎn)生不同的矩陣 ,下面就讓我們產(chǎn)生兩個(gè) 3*3的矩陣 ,這一章中 ,我們的學(xué)習(xí)就靠她們了 !!! A = pascal(3) A = 1 1 1 1 2 3 1 3 6 B = magic(3) B = 8 1 6 3 5 7 4 9 2 還有一個(gè) 3*2 的隨機(jī)矩陣 : C = fix(10*rand(3,2)) C = 9 4 2 8 6 7 看看列矩陣，行矩陣，以及常數(shù)的表達(dá)： u = [3。 1。 4] v = [2 0 —1] s = 7 產(chǎn)生的矩陣是： u = 3 1 4 v = 2 0 —1 s = 7 加減法 X = A + B X = 9 2 7 4 7 10 5 12 8 Y = X –A Y = 8 1 6 3 5 7 4 9 2 若二矩陣維數(shù)不統(tǒng)一，則會(huì)出錯(cuò)！ X = A + C Error using == + Matrix dimensions must agree. 向量的乘積與轉(zhuǎn)置 x = v*u x = 2 X = u*v X = 6 0 —3 2 0 —1 8 0 —4 X = B39。 X = 8 3 4 1 5 9 6 7 2 如ｘ與ｙ均是列向量，則 x*y 無解 ,但下二表達(dá)式卻可以 : x39。*y y39。*x 稱 內(nèi)積 或 點(diǎn)積 . 下面的語(yǔ)句產(chǎn)生 單位矩陣 eye(m,n) 若用 eye(n)則產(chǎn)生 n*n 的方陣 解線性方程 情況一 : x = A\u x = 10 —12 5 又如： X = A\B X = 19 –3 —1 —17 4 13 6 0 —6 情況二； y 是不同時(shí)刻 t 時(shí)的觀測(cè)值 : t = [0 .3 .8 ]39。 y = [.82 .72 .63 .60 .55 .50]39。 若函數(shù)形式是 :y(t)=