<p>　　下面我們以美國人口普查的數據來研究一下有關曲線擬合的問題（MatLab是別人的,教學文檔是別人的,例子也是別人的，我只是一個翻譯而已．．．）</p><p class="code">load census</p><p>這樣我們得到了兩個變量，<span class="explain">cdate</span>是1790至1990年的時間列向量(10年一次),<span class="explain">pop</span>是相應人口數列向量.</p><p>上一小節所講的多項式擬合可以用函數<span class="explain">polyfit()</span>來完成,數字指明了階數</p><p class="code">p = polyfit(cdate,pop,4) <br>  　　Warning: Matrix is close to singular or badly scaled. <br>  　　Results may be inaccurate. RCOND = 5.429790e–20 <br>  p = <br>  　　1.0e+05 * <br>  　　0.0000 –0.0000 0.0000 –0.0126 6.0020</p><p>產生警告的原因是計算中的cdata值太大，在計算中的Vandermonde行列式使變換產生了問題，解決的方法之一是使<span class="explain">數據標準化</span>．</p><h3>預處理：標準化數據</h3><p>數據的標準化是對數據進行縮放，以使以后的計算能更加精確，一種方法是使之成為0均值：</p><p class="code">sdate = (cdate – mean(cdate))./std(cdate)</p><p>現在再進行曲線擬合就沒事了！</p><p class="code">p = polyfit(sdate,pop,4) <br>  p = <br>  0.7047 0.9210 23.4706 73.8598 62.2285<br>  pop4 = polyval(p,sdate); <br>  plot(cdate,pop4,'–',cdate,pop,'+'), grid on </p><p><img src="image/data5.jpg" width="496" height="367"></p><p>在上面的數據標準化中，也可以有別的算法，如令１７９０年的人口數為０．</p><h3>余量分析</h3><p class="code">p1 = polyfit(sdate,pop,1);<br>  pop1 = polyval(p1,sdate);<br>  plot(cdate,pop1,'–',cdate,pop,'+')<br>  <img src="image/data6_1.jpg" width="447" height="351"> </p><p class="code">res1 = pop – pop1;<br>  figure, plot(cdate,res1,'+')<br>  <img src="image/data6_2.jpg" width="449" height="357"> </p><p class="code">p = polyfit(sdate,pop,2);<br>  pop2 = polyval(p,sdate);<br>  plot(cdate,pop2,'–',cdate,pop,'+')<br>  <img src="image/data6_3.jpg" width="454" height="353"> </p><p class="code">res2 = pop – pop2;<br>  figure, plot(cdate,res2,’+’)<br>  <img src="image/data6_4.jpg" width="448" height="354"> </p><p class="code">p = polyfit(sdate,pop,4);<br>  pop4 = polyval(p,sdate);<br>  plot(cdate,pop4,'–',cdate,pop,'+')<br>  <img src="image/data6_5.jpg" width="561" height="386"> </p><p class="code">res4 = pop – pop4;<br>  figure, plot(cdate,res4,'+')<br>  <img src="image/data6_6.jpg" width="465" height="382"> </p><p>　　可以看出，多項式擬合即使提高了階次也無法達到令人滿意的結果</p><h3>指數擬合</h3><p>　　從人口增長圖可以發現人數的增長基本是呈指數增加的，因此我們可以用年份的對數來進行擬合，這兒，年數是標準化后的！</p><p class="code">logp1 = polyfit(sdate,log10(pop),1); <br>  logpred1 = 10.^polyval(logp1,sdate); <br>  semilogy(cdate,logpred1,'–',cdate,pop,'+'); <br>  grid on</p><p><img src="image/data6_7.jpg" width="590" height="460"></p><p class="code">logres2 = log10(pop) –polyval(logp2,sdate);<br>  plot(cdate,logres2,'+')</p><p><img src="image/data6_8.jpg" width="472" height="372"> </p><p>上面的圖不令人滿意，下面，我們用二階的對數分析：</p><p class="code">logp2 = polyfit(sdate,log10(pop),2); <br>  logpred2 = 10.^polyval(logp2,sdate); <br>  semilogy(cdate,logpred2,'–',cdate,pop,'+'); <br>  grid on</p><p><img src="image/data6_9.jpg" width="593" height="459"></p><p class="code">r = pop – 10.^(polyval(logp2,sdate));<br>  plot(cdate,r,'+')</p><p><img src="image/data6_10.jpg" width="471" height="368"> </p><p>這種余量分析比多項式擬合的余量分析圖案要隨機的多（沒有很強的規律性），可以預見，隨著人數的增加，余糧所反映的不確定度也在增加，但總的說來，這種擬合方式要強好多！</p><h3>誤差邊界</h3><p>誤差邊界常用來反映你所用的擬合方式是否適用于數據，為得到誤差邊界，只需在<span class="explain">polyfit()中傳遞第二個參數，并將其送入polyval()</span>.<br>  下面是一個二階多項式擬合模型,年份已被標準化,下面的代碼用了2σ,對應于95%的可置信度:</p><p class="code">[p2,S2] = polyfit(sdate,pop,2); <br>  [pop2,del2] = polyval(p2,sdate,S2); <br>  plot(cdate,pop,'+',cdate,pop2,'g–',cdate,pop2+2*del2,'r:',... cdate,pop2–2*del2,'r:'),   <br>  grid on </p><p><img src="image/data7.jpg" width="519" height="406"></p><p>&nbsp;</p><h2>差分方程和濾波</h2><p> 　　MatLab中的差分和濾波基本都是對向量而言的，向量則是存儲取樣信號或序列的．<br>  　　函數<span class="explain">y = filter(b, a, x)</span>將用a,b描述的濾波器處理向量x,然后將其存儲在向量y中,   <br>  　　filter()函數可看為是一差分方程a1y(n)=b1*x(1)+b2*x(2)+...-a2*y(2)-...<br>  　　如有以下差分方程:y(n)=1/4*x(n)+1/4*x(n-1)+1/4*x(n-2)+1/4*x(n-3),則</p><p class="code">a = 1; <br>  b = [1/4 1/4 1/4 1/4]; </p><p>我們載入數據,取其第一列,并計算有:</p><p class="code">load count.dat<br>  x = count(:,1);<br>  y = filter(b,a,x);<br>  t = 1:length(x); <br>  plot(t,x,'–.',t,y,'–'), <br>  grid on <br>  legend('Original Data','Smoothed Data',2) </p><p><img src="image/data8.jpg" width="572" height="435"></p><p>實現所表示的就是濾波后的數據,它代表了4小時的平均車流量<br>  MatLab的信號處理工具箱中提供了很多用來濾波的函數,可用來處理實際問題!</p><h2>快速傅立葉變換(FFT)</h2><p>傅立葉變換能把信號按正弦展開成不同的頻率值,對于取樣信號,用的是離散傅立葉變換.<br>  FFT是離散傅立葉變換的一種高速算法,在信號和圖像處理中有極大的用處!</p><table width="100%" border="0" cellspacing="0" cellpadding="0" height="349">  <tr>    <td width="23%">fft</td>    <td width="77%">離散傅立葉變換</td>  </tr>  <tr>    <td width="23%">fft2</td>    <td width="77%">二維離散傅立葉變換</td>  </tr>  <tr>    <td width="23%">fftn</td>    <td width="77%">n維離散傅立葉變換</td>  </tr>  <tr>    <td width="23%">ifft</td>    <td width="77%">離散傅立葉反變換</td>  </tr>  <tr>    <td width="23%">ifft2</td>    <td width="77%">二維離散傅立葉反變換</td>  </tr>  <tr>    <td width="23%">ifftn</td>    <td width="77%">n維離散傅立葉反變換</td>  </tr>  <tr>    <td width="23%">abs</td>    <td width="77%">幅度</td>  </tr>  <tr>    <td width="23%">angle</td>    <td width="77%">相角</td>  </tr>  <tr>    <td width="23%">unwrap</td>    <td width="77%">相位按弧度展開,大于π的變換為2π的補角</td>  </tr>  <tr>    <td width="23%">fftshift</td>    <td width="77%">把零隊列移至功率譜中央</td>  </tr>  <tr>    <td width="23%">cplxpair</td>    <td width="77%">把數據排成復數對</td>  </tr>  <tr>    <td width="23%">nextpow2</td>    <td width="77%">下兩個更高的功率</td>  </tr></table><p>向量x的FFT可以這樣求:</p><p class="code">x = [4 3 7 –9 1 0 0 0]’<br>  y = fft(x)<br>  y = <br>  　　6.0000 <br>  　　11.4853　–2.7574i <br>  　–2.0000 　–12.0000i <br>  　–5.4853　　+11.2426i <br>  　　18.0000 <br>  　–5.4853 　–11.2426i <br>  　–2.0000　　+12.0000i <br>  　　11.4853 　+2.7574i </p><p>x雖然是實數,但y是復數,其中,第一個是因為它是常數相加的結果,第五個則對應于奈奎斯特頻率,后三個數是由于負頻率的影響,它們是前面三個數的共軛值!</p><p>下面,讓我們來驗證一下太陽黑子活動周期是11年!Wolfer數記錄了300年太陽黑子的數量及大小:</p><p class="code">load sunspot.dat <br>  year = sunspot(:,1); <br>  wolfer = sunspot(:,2); <br>  plot(year,wolfer) <br>  title(’Sunspot Data’)</p><p><img src="image/data9_1.jpg" width="618" height="485"></p><p>現在來看看其FFT:</p><p class="code">Y = fft(wolfer);</p><p>Y的幅度是功率譜,畫出功率譜和頻率的對應關系就得出了周期圖,去掉第一點,因為他只是所有數據的和,畫圖有:</p><p class="code">N = length(Y); <br>  Y(1) = []; <br>  power = abs(Y(1:N/2)).^2; <br>  nyquist = 1/2; <br>  freq = (1:N/2)/(N/2)*nyquist; <br>  plot(freq,power), <br>  grid on <br>  xlabel(’cycles/year’) <br>  title(’Periodogram’)</p><p><img src="image/data9_2.jpg" width="596" height="497"></p><p>上面的圖看起來不大方便,下面我們畫出頻譜-周期圖</p><p class="code">period = 1./freq; <br>  plot(period,power), <br>  axis([0 40 0 2e7]), <br>  grid on <br>  ylabel(’Power’) <br>  xlabel(’Period(Years/Cycle)’)</p><p><img src="image/data9_3.jpg" width="606" height="494"></p><p>為了得出精確一點的解,如下:</p><p class="code">[mp index] = max(power); <br>  period(index) <br>  ans = <br>  　　 11.0769</p><h3>變換后的幅度和相位</h3><p><span class="explain">abs()和angle()</span>是用來計算幅度和相位的</p><p>先創建一信號,再進行分析,unwarp()把相位大于π的變換為2π的補角:</p><p class="code">t = 0:1/99:1; <br>  x = sin(2*pi*15*t) + sin(2*pi*40*t);<br>  y = fft(x);<br>  m = abs(y); <br>  p = unwrap(angle(y)); <br>  f = (0:length(y)–1)'*99/length(y); <br>  subplot(2,1,1), plot(f,m), <br>  ylabel('Abs. Magnitude'), grid on <br>  subplot(2,1,2), plot(f,p*180/pi) <br>  ylabel('Phase [Degrees]'), grid on <br>  xlabel('Frequency [Hertz]') </p><p><img src="image/data9_4.jpg" width="644" height="483"></p><p>可以發現幅度曲線關于奈奎斯特頻率對稱,只有0-50Hz的信息是有用的!</p><h3>FFT的長度與速度</h3><p>可以為FFT加上第二個參數,告訴MatLab這是n點FFT.如y = fft(x,n),若x長度大于n,軟件自動補0,否則截取x.</p><p>若:</p><p>1. n為2的冪,軟件將執行基2快速傅立葉算法,這時的運算速度是最快的<br>  2. n為合數,軟件將n分解為素數來算,計算量與n的值有關.n為1013將比1000點的速度慢的多!<br>  3. n為素數,軟件執行DFT的公式,此時最慢!</p><p>&nbsp;</p><p>&nbsp;</p><p>&nbsp;</p></body></html>