成信工統(tǒng)計(jì)天氣預(yù)報(bào)實(shí)驗(yàn)指導(dǎo)03功率譜分析和濾波

實(shí)驗(yàn)三功率譜分析和濾波一、目的和要求：周期是大氣運(yùn)動(dòng)的基本理念之一，周期分析也是統(tǒng)計(jì)天氣分析中的基礎(chǔ)理論。功率譜分析的模型和濾波思想是最基本的周期分析統(tǒng)計(jì)模型，是必須理解和掌握的基本理論方法，是后續(xù)課程中許多天氣分析事實(shí)理解的基礎(chǔ)，也是畢業(yè)設(shè)計(jì)和論文中的基本分析方法。該方法作為天氣分析模型，在實(shí)際業(yè)務(wù)工作和科研工作中有極強(qiáng)的實(shí)用意義。通過該實(shí)驗(yàn)，能使學(xué)生深刻理解氣象時(shí)間序列周期分析的意義和周期分析的方法，為實(shí)際的業(yè)務(wù)和科研工作打下一定的基礎(chǔ)。二、實(shí)驗(yàn)的主要內(nèi)容：對(duì)原序列濾波前（或后）利用功率譜檢測序列的周期，并進(jìn)行顯著性檢驗(yàn)。本實(shí)驗(yàn)利用功率譜方法檢測原序列的周期，并進(jìn)行顯著性檢驗(yàn)；利用帶通濾波器提取所需要的周期分量，再利用功率譜檢驗(yàn)所提取周期的顯著性。三、步驟：3.1熟悉資料方法3.1.1資料本次實(shí)驗(yàn)使用國家氣候中心提供58年（1951年-2008年）夏季（6-8月）逐月降水資料。取得是長江中下游及江淮地區(qū)降水場第一主要因子的時(shí)間序列為例。3.1.2功率譜分析方法譜分析是時(shí)間序列在頻域上進(jìn)行分析的方法，亦稱頻譜分析或波譜分析。主要側(cè)重在時(shí)間序列波動(dòng)分析上。功率譜的概念是針對(duì)功率有限信號(hào)的，所表現(xiàn)的是單位頻帶內(nèi)信號(hào)功率隨頻率的變換情況。功率譜密度是指對(duì)于具有連續(xù)頻譜和有限平均功率的信號(hào)或噪聲，表示其頻譜分量的單位帶寬功率的頻率函數(shù)。氣象上常利用功率譜作周期分析，功率譜圖顯示不同頻率振動(dòng)的功率大小，同時(shí)也是方差貢獻(xiàn)的大小，因而可以從譜曲線中的譜值最大來確認(rèn)主要振動(dòng)及其對(duì)應(yīng)的周期。功率譜估計(jì)是指對(duì)一實(shí)測時(shí)間序列現(xiàn)實(shí)（t=1，2，…，n)的估計(jì)，可有兩種估計(jì)功率譜方法：離散功率譜估計(jì)和連續(xù)功率譜估計(jì)，本實(shí)驗(yàn)選擇連續(xù)功率譜估計(jì)。具體求連續(xù)功率譜估計(jì)方法步驟如下：(1)計(jì)算樣本落后自相關(guān)系數(shù)，m為最大步長，或最大落后時(shí)間長度，本實(shí)驗(yàn)m取18，通常m宜選擇在n/3～n/10之間。（為時(shí)間落后步長）(2)求粗譜估計(jì)波數(shù)與圓頻率的關(guān)系在取樣點(diǎn)數(shù)為m時(shí)有：波數(shù)k最大可取到m/2。習(xí)慣上令l=2k，波數(shù)l取值點(diǎn)可從0到m。(3)計(jì)算平滑功率譜也可（式中l(wèi)=0，1，…，m）(4)作譜圖，以波數(shù)l為橫軸，平滑功率譜密度估計(jì)值為縱坐標(biāo)作圖。波數(shù)l與周期關(guān)系為連續(xù)功率譜檢驗(yàn)包括紅色噪音過程和白色噪音過程，本實(shí)驗(yàn)選擇白色噪音過程。3.1.3濾波方法對(duì)氣象要素作譜分析過程中，一些規(guī)則周期占有很大的分量，但這種周期是眾所周知的，它的存在，壓低了其它周期的表現(xiàn)，一旦把它去掉之后，則可突出地表現(xiàn)其它周期的成分，即把感興趣的周期成分從原來的序列中識(shí)別和提取出來，這一過程就是濾波過程。實(shí)際上是原始序列經(jīng)一定的變換轉(zhuǎn)化為另一序列的過程。濾波的總體過程：輸入輸出x(t)濾波器y(t)（數(shù)字處理）其中x(t)為原始序列，y(t)為新的時(shí)間序列。要求：過濾系統(tǒng)具有時(shí)間不變性和穩(wěn)定性。濾波器可分為三類：用于從原序列中濾出低于某個(gè)頻率成分的濾波器，稱為低通濾波器；用于濾出高于某個(gè)頻率成分的濾波器稱為高通濾波器；濾出兩個(gè)頻率之間一個(gè)頻帶成分的濾波器稱為帶通濾波器。根據(jù)本實(shí)驗(yàn)的具體情況，選擇帶通濾波器，濾波范圍為2～8年。3.2編寫程序要求編寫主程序，其中包括資料讀入，范圍截取，子程序調(diào)用。Cccccccccccccccccccccccccccc（附程序，對(duì)關(guān)鍵部分標(biāo)志出）ccccccccccccccccccccccccc連續(xù)功率譜程序 parameter(n=50,m=n/3-1,ta=1.645,xx=12.59,pi=3.141592653589793) realx(n),r(0:m),ss(0:m),sr(0:m),sss(0:m),ssss(0:m),vv,ssr,sum,rc&tta,ave open(1,file='soi-1fil.txt') open(3,file='soi-1filsl1.txt')open(5,file='soi-1filtest1.txt') doi=1,n read(1,*)x(i) enddo close(1) **********************求方差***************** ex=0.0 doi=1,n ex=ex+x(i) enddoex=ex/n s=0.0 doi=1,n s=s+(x(i)-ex)**2 enddo s=sqrt(s/n)**********************求自相關(guān)系數(shù)***************** doj=0,m r(j)=0.0 doi=1,n-j r(j)=r(j)+(x(i)-ex)*(x(i+j)-ex)/(s*s) enddo r(j)=r(j)/(n-j)enddo*********************求sl***************************doi=0,m ssr=0.0dok=1,m-1 ssr=ssr+r(k)*(1+cos((pi*k)/m)*cos(i*pi*k/m))enddosr(i)=ssrif((i.ne.0).or.(i.ne.m))thenss(i)=(1.0/m)*(r(0)+sr(i))elsess(i)=1.0/2*m*(r(0)+sr(i)) endif enddodok=0,m write(3,'(f8.4)')ss(k) enddo close(2) **********檢驗(yàn)******************* sum=0.0 doi=0,m sum=sum+ss(i) enddo ave=sum/mtta=n-2rc=(-1+ta*sqrt(tta))/(n-1)vv=(2.0*n-m/2.0)/m write(*,*)ave,vv,rc,r(1) doi=0,m sss(i)=ave*(((1-r(1)**2))/((1+r(1)**2)-2*r(1)*cos(pi*i/m))) enddowrite(*,*)'ok' if(r(1)>rc)then doi=0,m ssss(i)=sss(i)*(xx/vv) enddopause elsedoi=0,mssss(i)=ave*(xx/vv) enddo endif doi=0,m write(5,*)ssss(i)enddoend濾波程序programhigh_filter PARAMETER(ND=50,m1=2,m2=8)parameter(nx=1,ny=1)DIMENSIONy(nd),y1(ND),x(nd),x1(ND),x2(ND)CC------------------------------------------------open(1,file='soi-7.txt')read(1,*)(x(i),i=1,nd)close(1)cc------------------------DO35IT=1,NDy1(IT)=X(IT)35CONTINUECALLFILTER(y1,ND,m1,x1) CALLFILTER(y1,ND,m2,x2)DO25IT=1,NDY(IT)=x2(IT)-x1(it)25CONTINUE30CONTINUEcc----------------------------open(2,file='soi-1fil.txt')write(2,'(f9.4)')(y(i),i=1,nd)END!********************低通濾波*********************************!************************************************!輸入一維序列Y(N)，其中N為樣本容量，M為要進(jìn)行光滑平均濾波的數(shù)目*!x1(N)為低通濾波后序列*!本程序中計(jì)算權(quán)重系數(shù)的子程序weight(wt,k)可以調(diào)整為其他分布*!比如:等權(quán)重系數(shù)，二項(xiàng)系數(shù)等*!*!****************************************************************** subroutinefilter(y,n,m,x1) real,dimension(n)::y,x1,x2 real,allocatable::wt(:),tp(:) x1=0;x2=0 if(mod(m,2)==1)then k=m/2 else k=m/2-1 endif allocate(wt(-k:k),tp(-k:k)) callweight(wt,k) doi=1,n if(i<=k)then sum=0 doj=-i+1,k sum=sum+wt(j) enddo doj=-i+1,k tp(j)=wt(j)/sum x1(i)=x1(i)+y(i+j)*tp(j) enddo elseif(i>n-k)then sum=0 doj=-k,n-i sum=sum+wt(j) enddo doj=-k,n-i tp(j)=wt(j)/sum x1(i)=x1(i)+y(i+j)*tp(j) enddo else doj=-k,k x1(i)=x1(i)+y(i+j)*wt(j) enddo endif enddoend!采用正態(tài)分布（高斯分布）計(jì)算滑動(dòng)平均系數(shù)!wt(-k:k)為權(quán)重系數(shù)，整個(gè)間隔范圍為m=2*k+1!******************************************************************** subroutineweight(wt,k) real,dimension(-k:k)::wt pi=3.1415926 sum=0 temp1=sqrt(2*pi) doi=-k,k wt(i)=exp(-(6.0*i/real(2*k+1))**2/2)/temp1 sum=sum+wt(i) enddo doi=-k,k wt(i)=wt(i)/sum enddo end3.3結(jié)果輸出(范例)范例所用資料濾波前后功率譜圖（圖1、圖2）僅供參考。圖1濾波前功率譜圖2濾波后功率譜（應(yīng)用EX