mfilt.f90

基于matlab的探地雷达信号处理开源程序

! Depending on imported flag 'MODE', applies:
! A 2-D moving-window median filter if MODE = 'MEDIFILT'
! or
! A 2-D moving-window mean filter if MODE = 'MEANFILT'
!
! ==> Data are imported after preprocessing by MATLAB or OCTAVE (see MFILT.M)
! ==> Results are exported for post-processing by MATLAB or OCTAVE (see MFILT.M)
!
! Author  : Andreas Tzanis
!           Department of Geophysics, University of Athens
! Created : 12 September 2003
!      

$debug
character mode*8
character buffer*255, infile*255, outfile*255
integer*2 status
integer   ns, ntr, nslow, nshigh,  ntrlow, ntrhigh, err_alloc
! Implements dynamic array allocation : Use only a F90 and higher compiler
real, allocatable :: dpad(:,:), dpadf(:,:), subvec(:)

! input and output file names are passed as arguments
call getarg(int2(1),buffer,status)
infile = adjustl(buffer(1:status))
call getarg(int2(2),buffer,status)
outfile = adjustl(buffer(1:status))
! Input file (INFILE) has been prepared by MATLAB or OCTAVE (or some other program)
open(1,file=infile,form='binary')
read(1) mode
read(1) ns, nslow, nshigh, ntr, ntrlow, ntrhigh
ih = ntrlow+ntr+ntrhigh
it = nslow+ns+nshigh
is = (ntrlow+ntrhigh)*(nslow+nshigh)
! Allocate data arrays 
if(mode.eq.'medifilt'.or.mode.eq.'MEDIFILT') then
    allocate( dpad(it,ih), dpadf(it,ih), subvec(is), STAT = err_alloc  )
elseif(mode.eq.'meanfilt'.or.mode.eq.'MEANFILT') then
    allocate( dpad(it,ih), dpadf(it,ih), STAT = err_alloc  )
else
    print*, ' Wrong filter MODE - No filter action'
	stop
endif
if (err_alloc .ne. 0) then 
    print *, "Memory allocation error - No further action"
     stop
endif
! Get data 
do i=1,it
read(1) (dpad(i,j), j=1,ih)
enddo
close(1,STATUS='delete')

percent = 0.
ddec    = 0.
!!!!!!!!!!!!!!!!    MEDIAN FILTER
if(mode.eq.'medifilt'.or.mode.eq.'MEDIFILT') then
do iy = nslow,ns+nslow

!!!!! Report progress
      percent= 100.0*(ns-(ns-iy-nslow+1))/ns
      if (percent > ddec) then
          write(*,'(''  Done '',i3,'' % '')')	int(percent)
          ddec = ddec + 10.0
      endif
!!!!!
        do ix = ntrlow,ntr+ntrlow
               ivec = 0
               do k = iy-nslow+1, iy+nshigh
                      do j = ix-ntrlow+1, ix+ntrhigh
                            ivec = ivec + 1
                            subvec(ivec) = dpad(k,j)
                      enddo
               enddo
               call mdian1(subvec,ivec,dpadf(iy,ix))
        enddo
enddo
!!!!!!!!!!!!!!!!    MEAN FILTER
elseif(mode.eq.'meanfilt'.or.mode.eq.'MEANFILT') then
do iy = nslow,ns+nslow

!!!!! Report progress
      percent= 100.0*(ns-(ns-iy-nslow+1))/ns
      if (percent > ddec) then
          write(*,'(''  Done '',i3,'' % '')')	int(percent)
          ddec = ddec + 10.0
      endif
!!!!!

        do ix = ntrlow,ntr+ntrlow
               ivec         = 0
               dpadf(iy,ix) = 0.0
               do k = iy-nslow+1, iy+nshigh
                      do j = ix-ntrlow+1, ix+ntrhigh
                            ivec = ivec + 1
                            dpadf(iy,ix) = dpadf(iy,ix) + dpad(k,j)
                      enddo
               enddo
               dpadf(iy,ix) = dpadf(iy,ix)/float(ivec)
        enddo
enddo
endif


! Done - export filtered data for post-processing by MATLAB or OCTAVE  
open(1,file=outfile,form='binary')
do i=1,ih
write(1) (dpadf(j,i), j=1,it)
enddo
close(1)

end

! *** F77 routines from Numerical Recipes to compute median
! *** Modified to rid them of GO TO's and numbering
 SUBROUTINE MDIAN1(X,N,XMED)
 DIMENSION X(N)
 CALL SORT(N,X)
 N2=N/2
 IF(2*N2.EQ.N)THEN
     XMED=0.5*(X(N2)+X(N2+1))
   ELSE
     XMED=X(N2+1)
  ENDIF
  RETURN
  END
      
 SUBROUTINE SORT(N,RA)
  DIMENSION RA(N)
   L=N/2+1
   IR=N
    do while(.TRUE.)
     IF(L.GT.1)THEN
       L=L-1
       RRA=RA(L)
     ELSE
       RRA=RA(IR)
       RA(IR)=RA(1)
       IR=IR-1
       IF(IR.EQ.1)THEN
         RA(1)=RRA
         RETURN
       ENDIF
     ENDIF
     I=L
     J=L+L
     do while (j.le.ir)
       IF(J.LT.IR)THEN
         IF(RA(J).LT.RA(J+1))J=J+1
       ENDIF
       IF(RRA.LT.RA(J))THEN
         RA(I)=RA(J)
         I=J
         J=J+J
       ELSE
         J=IR+1
       ENDIF
     enddo
     RA(I)=RRA
   enddo
   END