第一次接触滤波程序，不知道这个程序是不是这个意思？求帮忙看看

jolincai

程序如下，源程序应该是一个34年的逐日资料，it1=10,it2=30 意思是要从资料中滤出10-30天的波吗？我用上我的资料，是35年每个格点每年都有一个值，格点数为29*17，那么一共有29*17*35个数据，老师让我把它分离成2-9年的年际场和大于9年的年代际场。我就把程序里的时间改成35，it1改成2，,it2改成9，这样就可以得到2-9年的年际场了吗？可是改了参数怎么会出现forrtl severe 36的错误？，我的资料就是29*17*35的数据呀？如果我要得到大于9年的年代际场 时间上怎么设置，两个时间设置成it1=9,it2呢？35吗？

parameter(it=12419,iz=1,it1=10,it2=30,dt=1,ixy=81*81)

real x(ixy,it),y(ixy,it),w1,w2,h(it)

pi=3.1415926

w1=2*pi/real(it1)

w2=2*pi/real(it2)

      open(1,file='../data/rain.grd'

     & ,form='unformatted',access='direct',recl=ixy*4)

open(2,file='../data/rain-1030.grd',

     & form='unformatted',access='direct',recl=ixy*4)

irec=1

do k1=1,iz

print*,k1

do k=1,it

irec=(k-1)*iz+k1

read(1,rec=irec)(x(i,k),i=1,ixy)

enddo

do i=1,ixy

       call responseBF2(it,dt,w1,w2,a,b1,b2,h)

       call Bfilter2(it,x(i,1:it),y(i,1:it),a,b1,b2)

enddo

do k=1,it

irec=(k-1)*iz+k1

write(2,rec=irec)(y(i,k),i=1,ixy)

enddo

enddo

end

      subroutine Bfilter2(n,x,y,a,b1,b2)

C     Second Order Butterworth Band-Pass Filter

C         y(k)=S(j=0,N)a(j)*x(k-j)-S(j:1,L)b(j)*y(k-j)

C       where L=2 is called the order of the filter, and generally, N=L.

C

C **  Note that: Call the subroutine responseBF2 first before calling this subroutine.

C

C     Input parameters: n, x(n), a, b1, and b2

C        n: the number od data

C        x: the original series

C        a, b1 and b2 from the subroutine responseBF2

C     Output variable: y(n)

C        y: the final filtered series of x

C     Work array: y1(n)

C        y1: the initial filtered series of x, work array

c     By Dr. LI Jianping, March 8, 2000.

c-----*----------------------------------------------------6---------7--

      dimension x(n),y(n)

      dimension y1(n)                                !Work array

      y1(1)=0.

      y1(2)=0.

      do 10 i=3,n

        y1(i)=a*(x(i)-x(i-2))-(b1*y1(i-1)+b2*y1(i-2))

  10  continue

      y(n)=y1(n)

      y(n-1)=y1(n-1)

      do 20 i=n-2,1,-1

        y(i)=a*(y1(i)-y1(i+2))-(b1*y(i+1)+b2*y(i+2))

  20  continue

  30  continue

      return

      end

c-----*----------------------------------------------------6---------7--

      subroutine responseBF2(n,dt,w1,w2,a,b1,b2,h)

C     Frequency Response Function of Second-order Butterworth Band-pass Filter

C     

C **  Note that: Call the subroutine responseBF2 first before call the subroutine

C                Bfilter2() for the Second Order Butterworth Band-Pass Filter

C

C     Input parameters: n, dt, w1, and w2

C       dt: the sampling interval (e.g., dt=1. if you use daily data)

C       w1: lower (circular) cutoff frequency on the left of w0

C           w1=2*pi/t1 where t1 is corresponding period to w1

C       w2: upper (circular) frequency on the right of w0

C           w2=2*pi/t2 where t2 is corresponding period to w2

C       w0: central (circular) frequency

C           w0=sqrt(w1*w2), =2*pi/t0 where t0 is corresponding period to w0

C       w1<w0<w2, that is t1>t0>t2

C     Output variables: a, b1, b2, and h

C        h: =|w(z)|**2 the amplitude of the frequency response function w(z)

C           where z=exp(-i*omg*datt). Note that: For F90, z=exp((0,-omg)) should

C           be modified to F90 format, z=exp(cmplx(0,-omg)).

C           Here we compile it in F77 format.

C     Work array:

C        w: the frequency response function w(z), it is a complex array.

C     By Dr. LI Jianping, March 8, 2000.

      dimension h(n)

      complex w,z                       !Work array

      pi=3.1415926

      w0=sqrt(w1*w2)

c      dt=1.

      a1=sin(w1*dt)/(1.+cos(w1*dt))

      a2=sin(w2*dt)/(1.+cos(w2*dt))

      dw=2.*abs(a1-a2)

      omg2=4.*a1*a2

      c=4.+2.*dw+omg2

      a=2.*dw/c

      b1=2.*(omg2-4.)/c

      b2=(4.-2.*dw+omg2)/c

      do 10 i=1,n

        omg=2.*pi/float(i)

        omg=omg*dt

c        z=exp((0,-omg))!F77

        z=exp(cmplx(0,-omg)) !F90

        w=a*(1-z**2)/(1.+b1*z+b2*z**2)

        h(i)=abs(w)**2

  10  continue

      return

      end

jolincai

求助啊~~~不明真相