空间点的时间低频滤波问题，求达人们帮助啊

飘_Icy

要对全国夏季降水距平百分率各站点60-13年进行高斯滤波，也就是说每个站点都做一次时间上的滤波，我用下面的程序做出来了1000站的滤波结果，为了验证我又将其中的某站点作了单独的滤波，
可是两者滤波结果不一样，按理说应该是一样才对的啊！
难道要我每个站单独的作一维滤波么，那岂不是太浩大了点

PROGRAM stnFFTT
parameter(s1=1055,n1=54,m1=9)
dimension r(s1,n1),r1(s1,n1),x1(n1),y1(n1)
INTEGER i,sta(s1)
!----------------------------------------------------
        open(20,file='china-60-13-summer-pre-anomaly-percentage-stn.dat')
    do i=1,s1
            do j=1,n1
             READ(20,*)sta(i),r(i,j)
            enddo
        enddo
        close(20)
!--------------------------------------------------------
open(30,file='china-60-13-summer-pre-anomaly-percentage-stn1055-lvbo.dat')        
do  i=1,s1
    do j=1,n1
          x1(j)=r(i,j)
        call guassfilter_2(n1,m1,x1,y1)
           r1(i,j)=y1(j)
                 write(30,"(i8,f10.2)") sta(i),r1(i,j)
        enddo
enddo
end
!-------------------------------------------------
      subroutine guassfilter_2(n,m,x,y)
!     M-term Guassian-Type Filter
!    Input variables: n, x(n), m
!        m: the term number used to running mean
!           it must be an odd number.
!     Output variables: y(n),z(n-m-1)
!       y: the filtered series of x.
!     Work parameters and array: c, cgm and ck(-(m-1)/2:(m-1)/2)
!        c: a tunable parameter, generally, c>2.0.
!      cgm: variance of Guassian distribution.
!     By Dr. LI Jianping, April 6, 2001.

      dimension x(n),y(n)
      dimension xw((-(m-1)/2+1):(n+(m-1)/2)),ck(-(m-1)/2:(m-1)/2) !work array
      undef=3276.6
      pi=3.1415926
      c=2.15
      nl=(m-1)/2
      cgm=float(nl)/c

      x1=x(1)
      xn=x(n)

      if(x(1).eq.undef)x1=x(2)
      if(x(n).eq.undef)xn=x(n-1)
      do i=-nl+1,1
        xw(i)=x1
      enddo
      do i=2,n-1
        xw(i)=x(i)
      enddo
      do i=n,n+nl
        xw(i)=xn
      enddo
      c1=1./(cgm*sqrt(2.*pi))
      ck(0)=c1
      do 10 i=1,nl
        ck(i)=c1*exp(-i*i/(2.*cgm*cgm))
        ck(-i)=ck(i)
  10  continue
      do 20 i=1,n
        y(i)=0.
        do j=-nl,nl
          y(i)=y(i)+ck(j)*xw(i+j)
        enddo
  20  continue
      if(x(1).eq.undef)y(1)=undef
      if(x(n).eq.undef)y(n)=undef
      return
      end

已50136站为例，绿色是原始的多年pa，红色线条是单独做的滤波，而黑色曲线是我利用上面的程序做的结果，
虽然趋势大致一样，但是为什么这样啊，想不明白了~请求大神们帮助啊

飘_Icy

自己调了调终于发现错误了，原来数据赋值出现错误了，改进的新程序


open(30,file='china-60-13-summer-pre-anomaly-percentage-stn1055-lvbo.dat')       
do 70 i=1,s
     do j=1,n1
          x1(j)=r(i,j)
        enddo
!        print*,x1(1:n1)
        do j=1,n1
        call guassfilter_2(n1,m1,x1,y1)
           r1(i,j)=y1(j)
                 write(30,"(i8,f10.2)") sta(i),r1(i,j)
        enddo

  70  continue

end

nightelfxu

你好，我也正在做类似的工作。请问你是先进行尺度分离然后分别eof么？

xlmooo

好东西 留个爪印

我就是-猴子

谢谢楼主 分享

full◆mind◇

{:5_213:}{:5_213:}

full◆mind◇

{:5_213:}{:5_213:}

full◆mind◇

{:5_213:}{:5_213:}

飘逸雪

留爪，学习