📄 one.f90
字号:
!===================================================================
!****** ****
!*** Funciton : 1-D Filter For TSP **
!****** ****
!====================================================================
!
! Writer:ShenHongyan
! Date :2005.4.21
! :2005.7.20
! :2005.9.30
!
!*******************************
!********************************
program Filter_1D
character*20 :: InFile,OutFile1,OutFile2,OutFile4,OutFile5
integer*4 :: m1,m2,m,n,n0
integer :: Spec_Trace,Win_type
integer :: Func
integer*2 :: head(120)
real df
real :: Fmin,Fmax,Cut_F1,Cut_f2
InFile='Fx4.sgy'
OutFile1='TSP_Dx1.sgy'
OutFile2='Ampl1.dat'
OutFile4='Freq.dat'
OutFile5='x.txt'
m1=1
m2=1000
Spec_Trace=2
Win_type=2
Func=3
Fmin=30.0
Fmax=200.0
Cut_F1=10.0
Cut_f2=100.0
call One_Filter_FFT(InFile,OutFile2,OutFile4,OutFile5,m1,m2,&
Spec_Trace,Win_type,df,head,m,n0,n)
call One_Filter_IFFT(OutFile5,OutFile1,head,m,n0,n,m1,m2,df,Func,&
Fmin,Fmax,Cut_F1,Cut_f2)
end program
subroutine One_Filter_FFT(InFile,OutFile2,OutFile4,OutFile5,m1,m2,&
Spec_Trace,Win_type,df,head,m,n0,n)
implicit none
character*20 :: InFile,OutFile2,OutFile4,OutFile5
integer*4 :: L_head,Len
integer*4 :: m1,m2,m,n,n0,nw,n1,n2,Nmin,nx
integer :: Spec_Trace,Win_type
integer :: i,j,k,k1
integer*2 :: head(120)
real :: t_step,df
real :: ax,bx,cx
real :: time1,Time2
real,allocatable :: N_Point1(:),N_Point3(:),Ampl(:),&
Phase(:),N_Point_z1(:)
complex,allocatable :: N_Point_x(:),N_Point_z2(:)
!c================================================
! Read parameters of head
L_head=240
open(12,file=InFile,access='direct',recl=L_head,status='unknown')
read(12,rec=1) head
n0=head(58)
t_step=head(59)
close(12)
Time1=t_step
Time2=n0*t_step
t_step=t_step/1e6 ! sampling interval in s
df=1.0/(n0*t_step)
Nmin=1.0/df
if(Nmin>n0) then
nx=Nmin
else
nx=n0
endif
!=================================================
! the parameters of spectrum time window
n1=Time1/(1000*t_step)
n2=Time2/(1000*t_step)
if(n1<=0.0.or.n1>n0) then
if(n1==0.0) then
n1=1
else
! write(*,*) 'Time1 is error!!!'
stop
endif
endif
if(n2>n0.or.n2<=0.0) then
if(n2>n0) then
n2=n0
else
! write(*,*) 'Time2 is error!!!'
stop
endif
endif
if(n2<=n1) then
! write(*,*) 'Time2 is error!!!'
stop
endif
!=================================================
! Extending n
if(nx.le.64) then
n=64
elseif(nx.le.128) then
n=128
elseif(nx.le.256) then
n=256
elseif(nx.le.512) then
n=512
elseif(nx.le.1024) then
n=1024
elseif(nx.le.2048) then
n=2048
elseif(nx.le.4096 ) then
n=4096
elseif(nx.le.8192 ) then
n=8192
elseif(nx.le.16384 ) then
n=16384
end if
nw=n/2
!=================================================
allocate(N_Point1(n0),N_Point3(n),Ampl(n),&
Phase(n),N_Point_x(n),N_Point_z1(n),N_Point_z2(n))
Len=n0*4+L_head
open(15,file=InFile,access='direct',recl=Len,status='unknown')
open(22,file=OutFile2,status="replace")
open(30,file=OutFile4,status="replace")
open(33,file=OutFile5,status="replace")
! Read Sum trace(m)
j=0
do while(.not.eof(15))
j=j+1
read(15,rec=j) head
enddo
m=j
if(m2>m) m2=m
! write(*,*)"Sum Trace:",m2
! write(*,*)
! Time sampling
df=1.0/(n*t_step)
!=================================================
! Fast Fourier Transformation modual
k=0
do j=m1,m2
k=k+1
! Read seismic data
read(15,rec=k) head,(N_Point1(i),i=1,n0)
! Read spectrum anslysis data
if(j==Spec_Trace)then
if(Win_type==2) then
k1=0
do i=n1,n2
k1=k1+1
N_Point_z1(k1)=N_Point1(i)
enddo
endif
endif
! One-D filter data type change
do i=1,n0
N_Point_x(i)=CMPLX(N_Point1(i),0.0)
enddo
do i=n0+1,n
N_Point_x(i)=CMPLX(0.0,0.0)
enddo
! Fast Freior Tramsform
call fft(n,N_Point_x,-1)
do i=1,n
ax=REAL(N_Point_x(i))
bx=AIMAG(N_Point_x(i))
cx=sqrt(ax*ax+bx*bx)
N_Point3(i)=cx
if(cx==0.0) cx=1e-6
Ampl(i)=20.0*log10(cx)
if(ax==0.0) ax=1e-6
Phase(i)=atan(bx/ax)
enddo
!=================================================
! Output one trace frequence(txt)
if(Win_type==1) then
if(j==Spec_Trace)then
do i=1,nw
write(22,*) i*df,Ampl(i)
write(30,*) i*df,N_Point3(i)
enddo
endif
endif
do i=1,n
write(33,*) j,i,N_Point_x(i)
end do
enddo
!
!====================
! Output one trace frequence(txt)
do i=1,n
N_Point_z2(i)=CMPLX(N_Point_z1(i),0.0)
enddo
! Fast Freior Tramsform
call fft(n,N_Point_z2,-1)
do i=1,nw
ax=REAL(N_Point_z2(i))
bx=AIMAG(N_Point_z2(i))
cx=sqrt(ax*ax+bx*bx)
N_Point3(i)=cx
if(cx==0.0) cx=1e-6
Ampl(i)=20.0*log10(cx)
write(22,*) i*df,Ampl(i)
write(30,*) i*df,N_Point3(i)
enddo
write(*,*)
deallocate(N_Point1,N_Point3,N_Point_x,&
Ampl,Phase,N_Point_z1,N_Point_z2)
close(15)
close(22)
close(30)
close(33)
end
!===============================
!=======================
!***********
SUBROUTINE FFT(LX,CX,IG)
COMPLEX CX(lx),CW,CTEMP
L=LX/2
SC=SQRT(1./LX)
A=3.1415926535*IG
J=1
DO 30 I=1,LX
IF(I.LE.J) THEN
CTEMP=CX(J)*SC
CX(J)=CX(I)*SC
CX(I)=CTEMP
ENDIF
M=L
20 if(j.gt.m) then
J=J-M
M=M/2
IF(M.GE.1) GO TO 20
ENDIF
30 j=j+m
L=1
40 istep=2*l
AB=A/L
SA=0.0
CA=1.0
SB=SIN(AB)
CB=COS(AB)
cw=cmplx(ca,sa)
DO 50 M=1,L
DO 45 I=M,LX,ISTEP
CTEMP=CW*CX(I+L)
CX(I+L)=CX(I)-CTEMP
45 CX(I)=CX(I)+CTEMP
CALL SWING(CA,SA,CB,SB)
CW=CMPLX(CA,SA)
50 continue
L=ISTEP
IF(L.LT.LX) GO TO 40
RETURN
endsubroutine
!===============
!c==================================================
!***********
SUBROUTINE SWING(CA,SA,CB,SB)
D=CB-1.0
Z1=CA+D*CA-SA*SB
Z2=SA+D*SA+CA*SB
T=1.5-.5*(Z1*Z1+Z2*Z2)
CA=T*Z1
SA=T*Z2
RETURN
endsubroutine
!!========================
!!******
!==============================================================
!***** *****
!** Function: Filter Factor for Frequency Space **
!***** *****
!==============================================================
!
! Writer:ShenHonyan
! Date :2005.4.22
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
subroutine Sub_Factor(Func,n,Cut_F1,cut_F2,f1,f3,df,h)
implicit none
parameter Pi=3.1415926
integer :: i,n,Func
real :: f1,f2,f3,f4
real :: df,Cut_F1,Cut_F2,mx1,mx2,mx3,mx4
real :: h(n)
!======================================
Select case(Func)
! Low Pass
case(1)
f2=f1+Cut_F1
mx1=f1/df
mx2=f2/df
do i=1,n
if(i<=mx1.or.i>=n-mx1) then
h(i)=1.0
elseif(i>mx1.and.i<=mx2) then
h(i)=(sin(Pi*(i-mx2)/(2.0*(mx1-mx2))))**2
elseif(i>n-mx2.and.i<n-mx1) then
h(i)=(sin(Pi*((n-i)-mx2)/(2.0*(mx1-mx2))))**2
else
h(i)=0.0
endif
enddo
! High Pass
case(2)
f4=f3-Cut_F2
mx3=f3/df
mx4=f4/df
do i=1,n
if(i>=mx3.and.i<=n-mx3) then
h(i)=1.0
elseif(i>mx4.and.i<=mx3) then
h(i)=(sin(Pi*(i-mx4)/(2.0*(mx3-mx4))))**2
elseif(i>n-mx3.and.i<n-mx4) then
h(i)=(sin(Pi*((n-i)-mx4)/(2.0*(mx3-mx4))))**2
else
h(i)=0.0
endif
enddo
! Band Pass
case(3)
f2=f1-Cut_F1
f4=f3+Cut_F2
mx1=f1/df
mx2=f2/df
mx3=f3/df
mx4=f4/df
if(mx1==mx3) mx3=mx1+0.1
do i=1,n
if((i>=mx1.and.i<=mx3).or.(i>=n-mx3.and.i<=n-mx1)) then
h(i)=1.0
elseif(i>=mx2.and.i<=mx1) then
h(i)=(sin(Pi*(i-mx2)/(2.0*(mx1-mx2))))**2
elseif(i>=n-mx1.and.i<=n-mx2) then
h(i)=(sin(Pi*((n-i)-mx2)/(2.0*(mx1-mx2))))**2
elseif(i>=mx3.and.i<=mx4) then
h(i)=(sin(Pi*(i-mx4)/(2.0*(mx3-mx4))))**2
elseif(i>=n-mx4.and.i<=n-mx3) then
h(i)=(sin(Pi*((n-i)-mx4)/(2.0*(mx3-mx4))))**2
else
h(i)=0.0
endif
enddo
! Band Stop
case(4)
f2=f1+Cut_F1
f4=f3-Cut_F2
mx1=f1/df
mx2=f2/df
mx3=f3/df
mx4=f4/df
if(mx1==mx3) mx3=mx1+0.1
do i=1,n
if((i<=mx1.or.(i>=mx3.and.i<=n-mx3).or.i>=n-mx1)) then
h(i)=1.0
elseif(i>=mx1.and.i<=mx2) then
h(i)=(sin(Pi*(i-mx2)/(2.0*(mx1-mx2))))**2
elseif(i>=n-mx2.and.i<=n-mx1) then
h(i)=(sin(Pi*((n-i)-mx2)/(2.0*(mx1-mx2))))**2
elseif(i>=mx4.and.i<=mx3) then
h(i)=(sin(Pi*(i-mx4)/(2.0*(mx3-mx4))))**2
elseif(i>=n-mx3.and.i<=n-mx4) then
h(i)=(sin(Pi*((n-i)-mx4)/(2.0*(mx3-mx4))))**2
else
h(i)=0.0
endif
enddo
! Defaut
case default
write(*,*)"Filter Function select error!!!"
write(*,*)"Filter Function select error!!!"
write(*,*)"Filter Function select error!!!"
write(*,*)
stop
end select
endsubroutine
!===============
!c==================================================
!***********
subroutine One_Filter_IFFT(OutFile5,OutFile1,head,m,n0,n,m1,m2,df,Func,&
Fmin,Fmax,Cut_F1,Cut_f2)
implicit none
character*20 :: OutFile5,OutFile1
integer*4 :: L_head,Len
integer*4 :: m1,m2,m,n,n0
integer :: Func
integer :: i,j,k,k1,k2
integer*2 :: head(120)
real :: df
real :: ax,bx
real :: Fmin,Fmax,Cut_F1,Cut_f2
real,allocatable :: N_Point(:),N_Point2(:),h(:)
complex,allocatable :: N_Point_x(:),N_Point_y(:)
!
!=====================================
! Judge Fmax > Fmin ???
if((Func==3.or.Func==4).and.Fmin>=Fmax) then
! write(*,*) "Error! "
! write(*,*) " Fmin>=Fmax "
! write(*,*) "Error! "
! write(*,*) " Fmin>=Fmax "
! write(*,*)
stop
endif
! Judge Fmin >=0 and Fmax >=0 ???
if(Fmin<=0.or.Fmax<=0) then
! write(*,*) "Error! "
! write(*,*) " Fmin<=0 or Fmax<=0 "
! write(*,*) "Error! "
! write(*,*) " Fmin<=0 or Fmax<=0 "
! write(*,*)
stop
endif
!c================================================
allocate(N_Point(n0),N_Point2(n),N_Point_x(n),N_Point_y(n),h(n))
L_head=240
Len=n0*4+L_head
open(15,file=OutFile5,status='unknown')
open(20,file=OutFile1,access='direct',recl=Len,status="replace")
!=================================================
! Add Filter Factor
call Sub_Factor(Func,n,Cut_F1,Cut_F2,Fmin,Fmax,df,h)
k=0
do j=m1,m2,1
k=k+1
do i=1,n
read(15,*) k1,k2,N_Point_x(i)
end do
do i=1,n
N_Point_y(i)=h(i)*N_Point_x(i)
N_Point2(i)=h(i)*abs(ax)
enddo
!=================================================
! Inverse Fast Freior Tramsform
call fft(n,N_Point_y,1)
do i=1,n0
ax=REAL(N_Point_y(i))
bx=AIMAG(N_Point_y(i))
N_Point(i)=ax/n
enddo
!=================================================
! Output seismic data(SEG-Y)
write(20,rec=k) head,(N_Point(i),i=1,n0)
! if(mod(k,50)==1) write(*,*) k," Trace"
enddo
deallocate(N_Point,N_Point2,N_Point_x,N_Point_y,h)
close(15)
close(20)
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -