成信工统计天气预报实验指导03功率谱分析和滤波

实验三功率谱分析和滤波一、目的和要求：周期是大气运动的基本理念之一，周期分析也是统计天气分析中的基础理论。功率谱分析的模型和滤波思想是最基本的周期分析统计模型，是必须理解和掌握的基本理论方法，是后续课程中许多天气分析事实理解的基础，也是毕业设计和论文中的基本分析方法。该方法作为天气分析模型，在实际业务工作和科研工作中有极强的实用意义。通过该实验，能使学生深刻理解气象时间序列周期分析的意义和周期分析的方法，为实际的业务和科研工作打下一定的基础。二、实验的主要内容：对原序列滤波前（或后）利用功率谱检测序列的周期，并进行显著性检验。本实验利用功率谱方法检测原序列的周期，并进行显著性检验；利用带通滤波器提取所需要的周期分量，再利用功率谱检验所提取周期的显著性。三、步骤：3.1熟悉资料方法3.1.1资料本次实验使用国家气候中心提供58年（1951年-2008年）夏季（6-8月）逐月降水资料。取得是长江中下游及江淮地区降水场第一主要因子的时间序列为例。3.1.2功率谱分析方法谱分析是时间序列在频域上进行分析的方法，亦称频谱分析或波谱分析。主要侧重在时间序列波动分析上。功率谱的概念是针对功率有限信号的，所表现的是单位频带内信号功率随频率的变换情况。功率谱密度是指对于具有连续频谱和有限平均功率的信号或噪声，表示其频谱分量的单位带宽功率的频率函数。气象上常利用功率谱作周期分析，功率谱图显示不同频率振动的功率大小，同时也是方差贡献的大小，因而可以从谱曲线中的谱值最大来确认主要振动及其对应的周期。功率谱估计是指对一实测时间序列现实（t=1，2，…，n)的估计，可有两种估计功率谱方法：离散功率谱估计和连续功率谱估计，本实验选择连续功率谱估计。具体求连续功率谱估计方法步骤如下：(1)计算样本落后自相关系数，m为最大步长，或最大落后时间长度，本实验m取18，通常m宜选择在n/3～n/10之间。（为时间落后步长）(2)求粗谱估计波数与圆频率的关系在取样点数为m时有：波数k最大可取到m/2。习惯上令l=2k，波数l取值点可从0到m。(3)计算平滑功率谱也可（式中l=0，1，…，m）(4)作谱图，以波数l为横轴，平滑功率谱密度估计值为纵坐标作图。波数l与周期关系为连续功率谱检验包括红色噪音过程和白色噪音过程，本实验选择白色噪音过程。3.1.3滤波方法对气象要素作谱分析过程中，一些规则周期占有很大的分量，但这种周期是众所周知的，它的存在，压低了其它周期的表现，一旦把它去掉之后，则可突出地表现其它周期的成分，即把感兴趣的周期成分从原来的序列中识别和提取出来，这一过程就是滤波过程。实际上是原始序列经一定的变换转化为另一序列的过程。滤波的总体过程：输入输出x(t)滤波器y(t)（数字处理）其中x(t)为原始序列，y(t)为新的时间序列。要求：过滤系统具有时间不变性和稳定性。滤波器可分为三类：用于从原序列中滤出低于某个频率成分的滤波器，称为低通滤波器；用于滤出高于某个频率成分的滤波器称为高通滤波器；滤出两个频率之间一个频带成分的滤波器称为带通滤波器。根据本实验的具体情况，选择带通滤波器，滤波范围为2～8年。3.2编写程序要求编写主程序，其中包括资料读入，范围截取，子程序调用。Cccccccccccccccccccccccccccc（附程序，对关键部分标志出）ccccccccccccccccccccccccc连续功率谱程序 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)**********************求自相关系数***************** 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) **********检验******************* 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为要进行光滑平均滤波的数目*!x1(N)为低通滤波后序列*!本程序中计算权重系数的子程序weight(wt,k)可以调整为其他分布*!比如:等权重系数，二项系数等*!*!****************************************************************** 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!采用正态分布（高斯分布）计算滑动平均系数!wt(-k:k)为权重系数，整个间隔范围为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结果输出(范例)范例所用资料滤波前后功率谱图（图1、图2）仅供参考。图1滤波前功率谱图2滤波后功率谱（应用EX