Lanczos带通滤波程序问题

麦田_smile

c                          LANCZOS FILTER
c        _________________________________________________________
c        ndat    | No of time points
c        t1 &t2  | 15 & 90 for a 15-90 day band pass filter
c        nwt        | No. of weights (should be > 1.3/((1/t1)-(1/t2))
c        _________________________________________________________
         parameter(ndat=744,t1=12,t2=120,nwt=35,undef=-999.)
       parameter(lon=81,lat=9,m=lon*lat)
         real c(ndat),p(m,ndat),u(m),p2(m)
         real t(ndat),a(m,ndat),b(m,ndat)


         open(10,file='e:/lan/u850jp.grd',form='binary')
         open(101,file='e:/lan/u850jp_filter.grd',form='binary')

       read(10) ((p(i,it),i=1,m),it=1,ndat)

         do i = 1,m
         do n = 1, ndat
               t(n) = p(i,n)
               enddo ! n
         call lanczos(t,ndat,nwt,t2,t1,c)
         do n = 1,ndat  
           b(i,n) = c(n)
           enddo
         enddo  ! i

         do n = 1, ndat
              write(101) (b(i,n),i=1,m)
             enddo  ! n
       stop
       end
         subroutine lanczos(u,it,N,t2,t1,p)


c         Program for LANCZOS filtering (25/5/2001,Prince K xavier)
c        for given values of cut-in (fc1) and cut-out (fc2) frequency
c        and the number of weights (N)
c        Ref: Duchon (1979), J.Appl.Meteor.,18
c         *************************************************************

c        it        : number of time points (365|366)
c        idt        : time step
c        N        : number of weights
c        fc1        : cut-in frequency
c        fc2        : cut-out frequency
c        sigma        : sigma factors
c        u()        : input time series data
c        w()        : computed weights
c        y()        : filtered time series
c        R()        : Response function
c        *************************************************************
        real u(it),v(it+N*2), w(2*N+1), y(it+N*2), R(it),p(it)

        l = it+2*N

c        Defining frequencies       
c        --------------------
c        write(*,*)'subroutine!'
        fc1 = 1.0/t2
        fc2 = 1.0/t1

c        Applying the padding
c        --------------------
        do i = 1, l
        if(i.gt.it+N)then
c        v(i) = u(it+N)/(i-l+N+1)
        v(i) = 0.0
        endif
        if(i.le.N)then
c        v(i) = u(N+1)/(N-i+2)
        v(i) = 0.0
        else
        v(i)= u(i-N)
        endif
        enddo


c        Defining constants
c        ------------------
        pi = 22.0/7.0
        idt = 1
       
c        Checking the criteria for N
c        ---------------------------
        temp = 1.3/(fc2-fc1)
        if(N.lt.temp)then
        write(*,*)'Program can not be run'
        else

c        Calculating the weights
c        -----------------------
        do k = -N,N
        if(k.eq.0)then
            w(N+k+1) = 2*(fc2-fc1)
        else
            sigma = sin(pi*k/N)/(pi*k/N)
            t01 = sin(2*pi*fc2*k)/(pi*k)
            t02 = sin(2*pi*fc1*k)/(pi*k)
            w(N+k+1) = (t01-t02)*sigma
        endif
c        write(*,*)k,w(N+k+1)
        enddo

c        Calculating the Filtered time series
c        ------------------------------------       
c        NEW

        do i = 1, l
        y(i) = 0.0
          if((i.le.N).or.(i.gt.(it+N)))then
           y(i) = 0.0
          else
        do k = -N, N
        y(i) =  y(i)+w(N+k+1)*v(i+k*idt)
        enddo
          endif        
c        write(*,*)i,y(i)
        enddo

        do i = 1+N, it+N
        p(i-N) = y(i)
c        write(*,*)y(i)
        enddo

c        The Response Function
c        ---------------------       
        do i = 1, it
        R(i) = 0.0
        f=1.0/i
        do k = -N, N
        R(i) =  R(i)+w(N+k+1)*cos(2*pi*f*k)
        enddo        
c        write(*,*)f,R(i)
        enddo
        endif
        return
        end

麦田_smile

一直数组越界 唉 烦死了

王淑莉

我试了一下你的这个程序，是你Applying the padding部分应该用elseif，而不是endif之后再if

/xin村儿/

能否给个什么是Lanczos带通滤波的文字说明？可以滤掉哪一部分？

爆浆炸鸡

这个可以滤出大气瞬变波吗

zhuzhejl

求详细介绍啊啊

小虫小虫

  不大懂滤波