低通滤波程序和带通滤波程序都是针对格点资料的吗？站点资料的可以用吗？

jolincai

我有几份滤波程序 ，不过好像都是处理再分析资料的，我把OPEN语句改成打开站点的TXT文本就不行，不知道是读取方式不对还是程序不能用在站点资料上？有站点资料可以用的滤波程序吗？还是我的程序两者都可以用 只是要修改一下？下面我贴一下我有的程序：（程序应该没错）
1.带通滤波
  parameter(it=35,iz=1,it1=2,it2=9,dt=1,ixy=29*17)
  real x(ixy,it),y(ixy,it),w1,w2,h(it)
  integer k,k1
  pi=3.1415926
  w1=2*pi/real(it1)
  w2=2*pi/real(it2)
      open(1,file='35eratmin.grd'
     & ,form='unformatted',access='direct',recl=ixy)
open(2,file='35eratmin2-9.grd',
     & form='unformatted',access='direct',recl=ixy)
irec=1
do k1=1,iz
print*,k1
  do k=1,it
  read(1,rec=k)(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
  write(2,rec=k)(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



2.低通滤波
        parameter(ixy=12*11,it=12419,iz=60,iw=31)
        real x(ixy,it),y(ixy,it)
        open(1,file='../data/qv-s/quv.grd',form='unformatted',
     &  access='direct',recl=ixy*4)
        open(2,file='../data/qv-s/quv-31.grd',form='unformatted',
     &  access='direct',recl=ixy*4)      
        do kz=1,iz
        print*,kz,'/60'
        do k=1,it
        irec=(k-1)*iz+kz
        read(1,rec=irec)(x(i,k),i=1,ixy)
        enddo   
        do i=1,ixy
               call guassfilter_2(it,iw,x(i,1:it),y(i,1:it))
        enddo
        do k=1,it
        irec=(k-1)*iz+kz
        write(2,rec=irec)(y(i,k),i=1,ixy)
        enddo
        enddo
        close(1)
        close(2)
        end
        subroutine guassfilter_2(n,m,x,y)
c     M-term Guassian-Type Filter
c     Input variables: n, x(n), m
c        m: the term number used to running mean
c           it must be an odd number.
c     Output variables: y(n),z(n-m-1)
c        y: the filtered series of x.
c     Work parameters and array: c, cgm and ck(-(m-1)/2:(m-1)/2)
c        c: a tunable parameter, generally, c>2.0.
c      cgm: variance of Guassian distribution.
c     By Dr. LI Jianping, April 6, 2001.
      dimension x(n),y(n)
      dimension xw((-(m-1)/2+1):(n+(m-1)/2)),c k(-(m-1)/2:(m-1)/2)
!work array
      undef=-9.99e33
      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

lqouc

1,可以读txt。open里的form改一下，read也要改。
2.程序对任何数据通用。

不会的地方随便找本教材翻翻

jolincai

lqouc 发表于 2016-9-5 10:52
1,可以读txt。open里的form改一下，read也要改。
2.程序对任何数据通用。

请问下 文本TXT一般是顺序读取  而二进制是直接读取，那REC RECL都是对直接读取起作用的，那么直接不要RECL=ixy了吗？read语句直接写成read(1,*) .......就可以了吗？

lqouc

jolincai 发表于 2016-9-5 15:04
请问下 文本TXT一般是顺序读取  而二进制是直接读取，那REC RECL都是对直接读取起作用的，那么直接不要RE ...

都不要了。就可以了。

jolincai

lqouc 发表于 2016-9-5 15:41
都不要了。就可以了。

谢谢~~好了。可是不懂滤波出来的结果是什么意思，负值滤波后不一定还是负值，而且滤波后值很小了，请问这个数值大小和正负怎么去分析呢？

lqouc

负值变正和数值变小都是正常的，毕竟也算是平滑的一种（尤其低通）。
结果的分析有很多方法，几句话也说不清楚，建议找几篇研究目的相似的文章看看。

zwk

请问楼主低通滤波里的parameter(ixy=12*11,it=12419,iz=60,iw=31)，分别代表什么，ixy是空间x*y,it是时间（天数），iz是高度层次，然后iw是滤掉31天周期以上的信号，保留月份尺度信号的意思吗？

zwk

lqouc 发表于 2016-9-5 19:55
负值变正和数值变小都是正常的，毕竟也算是平滑的一种（尤其低通）。
结果的分析有很多方法，几句话也说不 ...

請問大神该楼主低通濾波程序裏的parameter(ixy=12*11,it=12419,iz=60,iw=31)，分別代表什麼，ixy是空間x*y,it是時間（天數），iz是高度層次，然後iw是濾掉31天周期以上的信號，保留月份尺度信號的意思嗎？

lqouc

zwk 发表于 2016-12-23 16:05
請問大神该楼主低通濾波程序裏的parameter(ixy=12*11,it=12419,iz=60,iw=31)，分別代表什麼，ixy是空間x* ...

是的，保留月尺度以及更长时间尺度的信号。

独行侠

我想请教一下，您的带通滤波程序中ixy应该是空间格点数吧？那为什么您定义的数组二维的而不是三维的呢？不应该是对某一格点上的it长度的时间序列进行滤波吗？