data_as.htm

对matlab的初学者也许有用

<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>