地震数据处理方法实训

褶积滤波

这里使用的子进程为实训手册中已有的子进程。
计算使用FORTRAN。
绘图使用MATLAB。

使用理想低通滤波器和带通滤波器进行滤波。对应滤波保留的频率范围为：
低通滤波器：0-70Hz
带通滤波器：10-80Hz

均使用时域下的滤波因子得到滤波器，所有滤波过程在时域之下进行。

program ch1_conv_fliter
real, dimension(100) :: X,r_low_end=0., coun=[(i,i=1,100)],r_band_end=0.
real,dimension (101) :: h_low=0., h_band=0., t1
real,dimension (101+100-1) :: r_low=0., r_band=0.
real :: wc1=70.0,wc2=10.0,wc3=80.0 ! 滤波器的滤波边界
complex, dimension(100,2) :: C   ! 定义一个复数数组C，用于表示X。
character(10) :: filename='INPUT1.DAT'
integer :: i, num=100
real, parameter :: pi=3.1415926535, dt=2
! 读取input1的数据并输出到屏幕
open(1, file=filename, form='formatted', status='old')
read (1,*)(X(i), i=1,100)
print *,'################################################'
print *,'This is the data in the ',filename ,'.'
print *,'################################################'
do i=1,100
PRINT *,X(i)
end do
print *,'################################################'

! 理想低通滤波器
low_pass_fliter: do i=1,101
t1(i) = dt*(i-100/2-1)/1000.
! t1(i) = dt*(i-1)/1000.
	if ((t1(i).lt.0.0001) .and. (t1(i).gt.-0.0001)) then
	h_low(i)=2*pi*wc1/pi
	else
	h_low(i)=sin(2*pi*wc1*t1(i))/(pi*t1(i))
	endif
end do low_pass_fliter

call CON(X,h_low,r_low,100,101,100+101-1)
! 得到最后的输出，r_low_end为最终的滤波结果
do i=1,100
r_low_end(i)=r_low(i+50)
end do

print *,'################################################'
print *,'This is the last result, low pass 70.'
print *,'################################################'
do i=1,100
PRINT *,r_low_end(i)
end do
print *,'################################################'
! 写入70HZ低通滤波器结果
! open(2,file='con_low_filter_70.dat', status='replace', action='write')
! write_70: do i=1,100
! 	write(2,100) coun(i),r_low_end(i)
	
! end do write_70


! 理想带通滤波器
band_pass_fliter: do i=1,101
	if ((t1(i).lt.0.0001) .and. (t1(i).gt.-0.0001)) then
	h_band(i)=2*pi*(wc3-wc2)/pi
	else
	h_band(i)=sin(2*pi*wc3*t1(i))/(pi*t1(i))-sin(2*pi*wc2*t1(i))/(pi*t1(i))
	endif
end do band_pass_fliter

call CON(X,h_band,r_band,100,101,100+101-1)
! 得到最后的输出，r_band_end为最终的滤波结果
do i=1,100
r_band_end(i)=r_band(i+50)
end do

print *,'################################################'
print *,'This is the last result，band pass 80-10hz'
print *,'################################################'
do i=1,100
PRINT *,r_band_end(i)
end do
print *,'################################################'
! 写入滤波器结果
open(3,file='ch1_result.dat', status='replace', action='write')
! write(3,*) "number,origin,low_pass_70,band_pass_80_10"
write_8: do i=1,100
	write(3,100) coun(i),(i-1.)*dt ,X(i),r_low_end(i),r_band_end(i) ! number,t,origin,low_pass_70,band_pass_80_10
100 format(F4.0,X,F6.0,x,2(F14.7,X), F14.7)
end do write_8
! 未经过截取的数据
open(4,file='ch1_result_no.dat', status='replace', action='write')
! write(3,*) "number,origin,low_pass_70,band_pass_80_10"
write_9: do i=1,200
	write(4,200) (i-1.)*dt ,r_low(i),r_band(i) ! number,t,origin,low_pass_70,band_pass_80_10
200 format(2(F14.7,X), F14.7)
end do write_9

end program ch1_conv_fliter


SUBROUTINE DFT(N,C,S)
	COMPLEX C(N),Y(512),W
	IF(S.EQ.1.0)Z=1.0
	IF(S.EQ.-1.0)Z=1.0/FLOAT(N)
	DO 20 I=1,N
	Y(I)=(0.0,0.0)
	DO 20 J=1,N
	A=COS(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))
	B=SIN(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))*(-S)
	W=CMPLX(A,B)
 20	Y(I)=Y(I)+C(J)*W
	DO 30 I=1,N
 30	C(I)=Y(I)*Z
	RETURN
	END

	SUBROUTINE CON(A,B,C,I,J,K)
	DIMENSION A(I),B(J),C(K)
	DO 1 K1=1,K
  1	C(K1)=0.0
	DO 2 I1=1,I
	DO 2 I2=1,J
	II=I1+I2-1
  2	C(II)=C(II)+A(I1)*B(I2)*0.004
	RETURN
	END

gfortran ch1_conv_fliter  # 编译
./a.out # 执行

clear
load ch1_result.dat
load ch1_result_no.dat
figure(1)
hold on
subplot(3,1,1)
plot(ch1_result(:,2), ch1_result(:,3))
title('Origin')

subplot(3,1,2)
plot(ch1_result(:,2), ch1_result(:,4))
title('low_pass-70')

subplot(3,1,3)
plot(ch1_result(:,2), ch1_result(:,5))
title('band_pass_10-80')

figure(2)
hold on
subplot(3,1,2)
plot(ch1_result_no(:,1), ch1_result_no(:,2))
title('low_pass-70')

subplot(3,1,3)
plot(ch1_result_no(:,1), ch1_result_no(:,3))
title('band_pass_10-80')

快变滤波

program ch2_fliter
real, dimension(100) :: X
real, dimension(:), allocatable:: X_fft_output,H_low,H_band,X_real, X_image
COMPLEX, dimension(:), allocatable:: C, C_low,C_band ! 定义一个复数数组C，用于表示X。C_low表示使用低通滤波器得到的最终频率域数据。
real, dimension(:,:), allocatable:: X_low, X_band
real :: wc1=70.0,wc2=10.0,wc3=80.0,dt=2.  ! 滤波器的滤波边界! 采样间隔2ms
character(10) :: filename='INPUT1.DAT'
integer :: i, num=100, nfft ! 频率域采样间隔
real, parameter :: pi=3.1415926535 
real :: df, k, f
! 读取input1的数据并输出到屏幕
open(1, file=filename, form='formatted', status='old')
read (1,*)(X(i), i=1,100)
print *,'################################################'
print *,'This is the data in the ',filename ,'.'
print *,'################################################'
do i=1,100
PRINT *,X(i)
end do
print *,'################################################'

! 创建文件,输出原始数据
open(unit=10, status='replace',action='write', file='ch2_result.dat')
! write(10,*)'Origin data of seis.'
! write(10,*)'t(ms),amplitude'
do i=1,num
t=(i-1)*dt
! write(10,300) t,X(i)
enddo
300 format(F4.0,',',F15.6)

dt=dt/1000. ! 转换单位为s
! 计算FFT点数
k=log(real(num))/log(2.)
if(num .gt. 2**k) k=int(k+1)
nfft=2**k
print *, 'FFT点数为：',nfft
! 定义可变数组大小
allocate(X_fft_output(nfft),H_low(nfft),H_band(nfft),X_real(nfft), X_image(nfft))
allocate(C(nfft), C_low(nfft),C_band(nfft),X_low(nfft,2), X_band(nfft,2))


! 频率域采样间隔
df=1./(nfft*dt)
print *, '频率域采样间隔为：', df ,'HZ'

print *,'---------------------------------------------'
print *, '数据重组,为了使用fft将输入数据补零'
print *,'---------------------------------------------'
! 进行数据重组，将X转换为fft需要的数组C
real_to_complex: do i=1,nfft
if (i.le.num) then
	C(i)=CMPLX(X(i), 0.0)
else
	C(i)=CMPLX(0.0, 0.0)
endif
end do real_to_complex

print *,'---------------------------------------------'
print *,'补零后的复数数组：'
print *, C
print *,'---------------------------------------------'

! 进行FFT
print *,'---------------------------------------------'
print *, '执行FFT'
print *,'---------------------------------------------'
call FFT(nfft,C,1.0)
! 得到实部和虚部
complex_to_real: do i=1,nfft
 	X_real(i)=real(C(i))
	X_image(i)=aimag(C(i))
end do complex_to_real
print *,'FFT后的频谱数据。'
print *,'-real-------------------------------------------'
print *, X_real
print *,'-image-------------------------------------------'
print *, X_image
print *,'---------------------------------------------'

print *,'---------------------------------------------'
print *, '设计滤波器'
print *,'---------------------------------------------'

! 设置频率域滤波器
! 低通滤波器
low_pass_filter: do i=1,(nfft/2+1)
	f=(i-1)*df
	if((f.le.wc1).and.(f.ge.(-wc1))) then 
		H_low(i)=1.
	else
		H_low(i)=0.
	end if
end do low_pass_filter
low_pass_filter2: do i=(nfft/2)+2,nfft
	f=(i-1)*df
	H_low(i)=H_low(nfft-i+2)
end do low_pass_filter2

! 带通滤波器
band_pass_filter: do i=1,(nfft/2+1)
	f=(i-1)*df
	if((f.le.wc3.and.f.ge.wc2) .or. (f.ge.(-wc3).and.(f.le.(-wc2)))) then 
		H_band(i)=1.
	else
		H_band(i)=0.
	end if
end do band_pass_filter
band_pass_filter2: do i=(nfft/2)+2,nfft
	f=(i-1)*df
	H_band(i)=H_band(nfft-i+2)
end do band_pass_filter2


100 format(I4, ',' , 3(f15.6,','),f15.6)
! 进行频率域滤波
print *,'---------------------------------------------'
print *, '执行频率域滤波'
print *,'---------------------------------------------'
! 低通滤波,得到频率域结果
! write(10,*) "low pass frequency"
! write(10,*) "i,frequence,real,image,amplitude"
do i=1,nfft
	f=(i-1)*df
	X_low(i,1)=X_real(i)*H_low(i)
	X_low(i,2)=X_image(i)*H_low(i)
	f=(i-1)*df
	amp=sqrt(X_low(i,1)*X_low(i,1)+X_low(i,2)*X_low(i,2))*dt
	! write(10,100) i, f, X_low(i,1),X_low(i,2),amp
	C_low(i)=complex(X_low(i,1),X_low(i,2)) ! 将X转化为复数数组
enddo 

! 带通滤波,得到频率域结果
! write(10,*) "band pass frequency"
! write(10,*) "i,frequence,real,image,amplitude"
do i=1,nfft
	f=(i-1)*df
	X_band(i,1)=X_real(i)*H_band(i)
	X_band(i,2)=X_image(i)*H_band(i)
	f=(i-1)*df
	amp=sqrt(X_band(i,1)*X_band(i,1)+X_band(i,2)*X_band(i,2))*dt
	! write(10,100) i, f, X_band(i,1),X_band(i,2),amp
	C_band(i)=complex(X_band(i,1),X_band(i,2)) ! 将X转化为复数数组
enddo 
print *,'---------------------------------------------'
print *, '频率域滤波执行完成'
print *,'---------------------------------------------'


! write(10,*) "low pass time"
! write(10,*) "i,time(ms),origin, amplitude_low,amplitude_band"

! 进行反FFT
print *,'---------------------------------------------'
print *, '执行反FFT变换，得到时间域'
print *,'---------------------------------------------'
call FFT(nfft,C_low,-1.0)
call FFT(nfft,C_band,-1.0)

do i=1,num
t=(i-1)*dt*1000.
write(10,200) i, t, X(i), real(C_low(i)), real(C_band(i))
enddo
200 format(I4,x,F15.6,x,2(F15.6,x),F15.6)
print *,'---------------------------------------------'
print *, '写入完成'
print *,'---------------------------------------------'


deallocate(X_fft_output)
end program ch2_fliter


SUBROUTINE DFT(N,C,S)
	COMPLEX C(N),Y(512),W
	IF(S.EQ.1.0)Z=1.0
	IF(S.EQ.-1.0)Z=1.0/FLOAT(N)
	DO 20 I=1,N
	Y(I)=(0.0,0.0)
	DO 20 J=1,N
	A=COS(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))
	B=SIN(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))*(-S)
	W=CMPLX(A,B)
 20	Y(I)=Y(I)+C(J)*W
	DO 30 I=1,N
 30	C(I)=Y(I)*Z
	RETURN
	END

	SUBROUTINE CON(A,B,C,I,J,K)
	DIMENSION A(I),B(J),C(K)
	DO 1 K1=1,K
  1	C(K1)=0.0
	DO 2 I1=1,I
	DO 2 I2=1,J
	II=I1+I2-1
  2	C(II)=C(II)+A(I1)*B(I2)*0.004
	RETURN
	END

	SUBROUTINE FFT(LX,CX,SIGNI)
	COMPLEX CX(LX),CARG,CEXP,CW,CTEMP
	J=1
	SC=1.0/LX
	IF(SIGNI.EQ.1.0)SC=1.0
	SIG=-SIGNI
	DO 30 I=1,LX
	IF(I.GT.J)GO TO 10
	CTEMP=CX(J)*SC
	CX(J)=CX(I)*SC
	CX(I)=CTEMP
 10     M=LX/2
 20	IF(J.LE.M)GO TO 30
	J=J-M
	M=M/2
	IF(M.GE.1)GO TO 20
 30	J=J+M
	L=1
 40	ISTEP=2*L
	DO 50 M=1,L
	CARG=(0.0,1.0)*(3.14159265*SIG*(M-1))/L
	CW=CEXP(CARG)
	DO 50 I=M,LX,ISTEP
	CTEMP=CW*CX(I+L)
	CX(I+L)=CX(I)-CTEMP
 50	CX(I)=CX(I)+CTEMP
	L=ISTEP
	IF(L.LT.LX)GO TO 40
	RETURN
	END

clear
load ch2_result.dat

figure(2)
hold on
subplot(3,1,1)
plot(ch2_result(:,2), ch2_result(:,3))
title('Origin')

subplot(3,1,2)
plot(ch2_result(:,2), ch2_result(:,4))
title('lowpass 70')

subplot(3,1,3)
plot(ch2_result(:,2), ch2_result(:,5))
title('bandpass 10-80')

褶积滤波与递归滤波的比较

program ch3_compare
implicit none
real :: sum_ax, sum_bx
real, dimension(50) :: X, con_fliter_result_end, con_fliter_zero_result_end, X_digui, X_digui_zero
real, parameter :: dt=2
! 褶积滤波的滤波因子
real, dimension(21) :: hn,hn_r
integer :: i, j, p, numht_zero, num_con_fliter_zero_result, num_con_fliter_result
real, dimension(:), allocatable:: ht_zero, con_fliter_result, con_fliter_zero_result

real, dimension(5) :: a
real, dimension(4) :: b
real(kind=4) :: t1, t1e, t2, t2e, t3, t3e, t4, t4e  ! 定义运行时间

hn=[1.0, 5.254, 11.6, 13.7, 8.47, 0.775, -3.325, -2.764, &
-0.364, 1.099, 0.97, 0.15, -0.37, -0.344, -0.06, -0.126, 0.122, 0.025, -0.042, -0.043,0.0]

a=[1,4,6,4,1]
b=[-1.254,0.987,-0.341,0.0524]

! 返序的hn，用于制作零相位滤波因子
do i=1,size(hn)
    hn_r(size(hn)-i+1)=hn(i)
enddo
open(unit=2, file='ch3_test_result.dat', status='replace', action='write')
open(unit=3, file='ch3_result.dat', status='replace', action='write')
open(1, file='INPUT3.DAT', form='FORMATTED', status='OLD')
read(1,*) (X(i), i=1,50)

numht_zero=2*size(hn)-1
allocate(ht_zero(numht_zero))

! 计算零相位滤波因子ht_zero，其为h(t)*h(-t)
call CON(hn, hn_r, ht_zero, size(hn), size(hn_r), numht_zero)
! 开启零相位褶积滤波结果数组
num_con_fliter_zero_result=size(X)+numht_zero-1
allocate(con_fliter_zero_result(num_con_fliter_zero_result))
! 输出原始数据
write(2,*) 'Original data.'
do i=1, size(X)
    write (2,100) i, dt*(i-1)/1000., X(i)
enddo
call cpu_time(t1)


! 计算零相位褶积滤波
do p=1,1000000
call CON(X, ht_zero, con_fliter_zero_result, size(X), size(ht_zero), num_con_fliter_zero_result)
! write(2,*) 'Zero phase convolution filter result1.'
do i=(numht_zero-1)/2+1, num_con_fliter_zero_result-(numht_zero-1)/2
    ! write (2,100) i-(numht_zero-1)/2-1+1, dt*(i-(numht_zero-1)/2-1)/1000., con_fliter_zero_result(i) ! 舍去两端(size(ht)-1)/2个点
	con_fliter_zero_result_end(i-((numht_zero-1)/2+1)+1)=con_fliter_zero_result(i)
enddo
end do
call cpu_time(t1e)


! 计算褶积滤波
num_con_fliter_result=size(X)+size(hn)-1
allocate(con_fliter_result(num_con_fliter_result))
call cpu_time(t2)
do p=1,1000000
call CON(X, hn, con_fliter_result, size(X), size(hn), num_con_fliter_result)
! write(2,*) 'convolution filter result.'
do i=(size(hn)-1)/2+1, num_con_fliter_result-(size(hn)-1)/2
    ! write (2,100) i-(size(hn)-1)/2-1+1, dt*(i-(size(hn)-1)/2-1)/1000., con_fliter_result(i) ! 舍去两端(size(ht)-1)/2个点
	con_fliter_result_end(i+1-((size(hn)-1)/2+1))=con_fliter_zero_result(i)
enddo
enddo
call cpu_time(t2e)


! 创建递归滤波
call cpu_time(t3)
do p=1,1000000
do i=1,size(X)
sum_ax=0.
sum_bx=0.
	if(i.lt.(size(a)-1)) then
		if(i.eq.1) then
			sum_ax=a(i)*X(i)
			sum_bx=0.
		else
			do j=1,i
				sum_ax=sum_ax+X(j)*a(i+1-j)
			end do

			do j=1,i-1
				sum_bx=sum_bx+X_digui(j)*b(i-j)
			end do
		endif
	else 
		do j=1,size(a)
			sum_ax=sum_ax+X(i-j+1)*a(j)
		end do
		do j=1,size(b)
			sum_bx=sum_bx+X_digui(i-j)*b(j)
		end do
	endif
 	X_digui(i)=sum_ax-sum_bx
enddo
enddo
call cpu_time(t3e)

! write (2, *) 'Digui filter data'
! do i=1,size(X_digui)
!     write (2,100) i, (i-1)*dt/1000., X_digui(i)
! enddo


call cpu_time(t4)
! 创建零相位递归滤波
do p=1,1000000
do i=size(X),1,-1
	sum_ax=0.
	sum_bx=0.
	if(i.gt.(size(X)-size(a)+2)) then
		if(i.eq.size(X)) then
			sum_ax=a(1)*X_digui(i)
			sum_bx=0.
		else
			do j=1,size(X_digui)-i+1
				sum_ax=sum_ax+X_digui(i+j-1)*a(j)
			end do

			do j=1,size(X_digui)-i
				sum_bx=sum_bx+X_digui(i+j)*b(j)
			end do
		endif
	else 
		do j=1,size(a)
			sum_ax=sum_ax+X_digui(i+j-1)*a(j)
		end do
		do j=1,size(b)
			sum_bx=sum_bx+X_digui(i+j)*b(j)
		end do
	endif
 	X_digui_zero(i)=sum_ax-sum_bx
enddo
enddo
call cpu_time(t4e)

! write (2, *) 'Digui zero filter data'
! do i=1,size(X_digui)
!     write (2,100) i, (i-1)*dt/1000., X_digui_zero(i)
! enddo


do i=1,size(X_digui)
    write(3,200) i, (i-1)*dt,X(i), con_fliter_result_end(i), con_fliter_zero_result_end(i), X_digui(i),  X_digui_zero(i)
enddo
100 format(I4, 2(',', F15.6))
200 format(I4,x,F5.1,5(x, F15.4))
write(2,300) t1e-t1, t2e-t2, t3e-t3, t4e-t4
300 format(4(F29.27,x)) 
end program ch3_compare

	SUBROUTINE CON(A,B,C,I,J,K)
	DIMENSION A(I),B(J),C(K)
	DO 1 K1=1,K
  1	C(K1)=0.0
	DO 2 I1=1,I
	DO 2 I2=1,J
	II=I1+I2-1
  2	C(II)=C(II)+A(I1)*B(I2)*0.004
	RETURN
	END

clear
load ch3_result.dat

figure(2)
hold on
subplot(5,1,1)
plot(ch3_result(:,2), ch3_result(:,3))
title('Origin')

subplot(5,1,2)
plot(ch3_result(:,2), ch3_result(:,4))
title('褶积滤波')

subplot(5,1,3)
plot(ch3_result(:,2), ch3_result(:,5))
title('零相位褶积滤波')

subplot(5,1,4)
plot(ch3_result(:,2), ch3_result(:,6))
title('递归滤波')

subplot(5,1,5)
plot(ch3_result(:,2), ch3_result(:,7))
title('零相位递归滤波')

设计高通滤波因子进行频域分析

program ch4_highpass
real, parameter :: dt=4., pi=3.1415926535, wc=30.
real, dimension(101) :: ht, t, f, hwr, hwi
complex, dimension(101) :: hw   ! 定义一个复数数组hw，用于表示频率域的滤波因子。
open(6,file='ch4_result.dat', status='replace',action='write')

! 设计时域滤波因子ht
df=1000/dt/size(f)
do i=1,size(ht)
    t(i)=dt*(i-1.-(size(ht)-1.)/2.)/1000.
	! t(i)=dt*(i-1.)/1000.
    f(i)=df*(i-1)
	! if (t(i).ne.0.) then
	! 	f(i)=1./t(i)
	! else
	! 	f(i)=1./dt*1000.
	! endif
    if (t(i).ne.0) then
        ht(i)=0.-1.*sin(2*pi*t(i)*wc)/pi/t(i)! /(pi*t(i))
        hw(i)=CMPLX(ht(i), 0.0)
    elseif(t(i).eq.0) then
        ht(i)=1000.-2.*wc
        hw(i)=CMPLX(ht(i), 0.0)
    else
        continue
    endif
enddo

call DFT(size(f), hw, 1.0)

complex_to_real: do i=1,size(f)
 	hwr(i)=real(hw(i))
	hwi(i)=aimag(hw(i))
end do complex_to_real

! 输出：序号，时间，时域滤波因子，频率，频率域滤波器
do i=1,size(f)
    write (6,100) i, t(i), ht(i), f(i), sqrt(hwr(i)**2+hwi(i)**2)
enddo
100 format(I4,6(x,F20.4))


end program ch4_highpass


	SUBROUTINE DFT(N,C,S)
	COMPLEX C(N),Y(512),W
	IF(S.EQ.1.0)Z=1.0
	IF(S.EQ.-1.0)Z=1.0/FLOAT(N)
	DO 20 I=1,N
	Y(I)=(0.0,0.0)
	DO 20 J=1,N
	A=COS(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))
	B=SIN(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))*(-S)
	W=CMPLX(A,B)
 20	Y(I)=Y(I)+C(J)*W
	DO 30 I=1,N
 30	C(I)=Y(I)*Z
	RETURN
	END

clear
load ch4_result.dat
load ch4_test.dat

figure(2)
hold on
subplot(2,1,1)
plot(ch4_result(:,2), ch4_result(:,3))
title('Origin')

subplot(2,1,2)
plot(ch4_result(:,4), ch4_result(:,5))
title('frequency')

不同信号的频谱分析

program ch5_freanalysis
integer, parameter :: N=128, dt=4, pi=3.1415926535
integer :: k, nfft
complex, dimension(N,2) :: x_per_fre, x_pulse_fre, x_seis_fre
real, dimension(N) :: x_per, x_pulse, x_seis
real :: df,f,x_per_amp, x_pulse_amp, x_seis_amp

open(1, file='WAVE.DAT', form='formatted', status='old')
read(1,*) (x_seis(i), i=1,N)

do i=1,N
    x_per(i)=sin(2.*pi*i/64.)
    if(i.eq.1) then
        x_pulse(i)=1
    else 
        x_pulse(i)=0   
    endif
    x_per_fre(i,:)=CMPLX(x_per(i),0.0)
    x_pulse_fre(i,:)=CMPLX(x_pulse(i),0.0)
    x_seis_fre(i,:)=CMPLX(x_seis(i),0.0)
enddo


call FFT(N, x_per_fre, 1.0)
call FFT(N, x_pulse_fre, 1.0)
call FFT(N, x_seis_fre, 1.0)
open(3, file='ch5_result.dat',status='replace',action='write')

df = 1000./dt/N
do i=1,N
    f = (i-1)*df
    x_per_amp=sqrt(real(x_per_fre(i,1))**2+aimag(x_per_fre(i,2))**2)
    x_pulse_amp=sqrt(real(x_pulse_fre(i,1))**2+aimag(x_pulse_fre(i,2))**2)
    x_seis_amp=sqrt(real(x_seis_fre(i,1))**2+aimag(x_seis_fre(i,2))**2)
    write(3,100) i, f, x_per_amp, x_pulse_amp, x_seis_amp
enddo
100 format(I4, 4(x,F15.4))
end program ch5_freanalysis





	SUBROUTINE FFT(LX,CX,SIGNI)
	COMPLEX CX(LX),CARG,CEXP,CW,CTEMP
	J=1
	SC=1.0/LX
	IF(SIGNI.EQ.1.0)SC=1.0
	SIG=-SIGNI
	DO 30 I=1,LX
	IF(I.GT.J)GO TO 10
	CTEMP=CX(J)*SC
	CX(J)=CX(I)*SC
	CX(I)=CTEMP
 10     M=LX/2
 20	IF(J.LE.M)GO TO 30
	J=J-M
	M=M/2
	IF(M.GE.1)GO TO 20
 30	J=J+M
	L=1
 40	ISTEP=2*L
	DO 50 M=1,L
	CARG=(0.0,1.0)*(3.14159265*SIG*(M-1))/L
	CW=CEXP(CARG)
	DO 50 I=M,LX,ISTEP
	CTEMP=CW*CX(I+L)
	CX(I+L)=CX(I)-CTEMP
 50	CX(I)=CX(I)+CTEMP
	L=ISTEP
	IF(L.LT.LX)GO TO 40
	RETURN
	END

load ch5_result.dat

figure(5)
hold on
subplot(3,1,1)
plot(ch5_result(:,2), ch5_result(:,3))
title('周期')
xlabel('frequency Hz')
ylabel('amplitude')

subplot(3,1,2)
plot(ch5_result(:,2), ch5_result(:,4))
title('脉冲')
xlabel('frequency Hz')
ylabel('amplitude')

subplot(3,1,3)
plot(ch5_result(:,2), ch5_result(:,5))
title('地震波')
xlabel('frequency Hz')
ylabel('amplitude')

FFT补零的影响分析

program ch6_zeros

integer, parameter :: N=60, dt=4, pi=3.1415926535, N1=64, N2=128
complex, dimension(N,2) :: x_per_fre, x_seis_fre
complex, dimension(N1,2) :: x_per_fre1, x_seis_fre1
complex, dimension(N2,2) :: x_per_fre2, x_seis_fre2
real, dimension(N) :: x_per, x_seis
real :: df,f,x_per_amp, x_seis_amp

open(1, file='WAVE.DAT', form='formatted', status='old')
read(1,*) (x_seis(i), i=1,N)
! 设置存储文件
open(3, file='ch6_result_nozero.dat',status='replace',action='write')
open(4, file='ch6_result_64.dat',status='replace',action='write')
open(5, file='ch6_result_128.dat',status='replace',action='write')

! No zero
do i=1,N
    x_per(i)=sin(2.*pi*i/30.)
    x_per_fre(i,:)=CMPLX(x_per(i),0.0)
    x_seis_fre(i,:)=CMPLX(x_seis(i),0.0)
enddo
call DFT(N, x_per_fre, 1.0)
call DFT(N, x_seis_fre, 1.0)
df = 1000./dt/N
do i=1,N
    f = (i-1)*df
    x_per_amp=sqrt(real(x_per_fre(i,1))**2+aimag(x_per_fre(i,2))**2)
    x_seis_amp=sqrt(real(x_seis_fre(i,1))**2+aimag(x_seis_fre(i,2))**2)
    write(3,100) i, f, x_per_amp, x_seis_amp
end do

! TO 64 with zero
do i=1,N1
    if (i.le.N) then
        x_per_fre1(i,:)=CMPLX(x_per(i),0.0)
        x_seis_fre1(i,:)=CMPLX(x_seis(i),0.0)
    else
        x_per_fre1(i,:)=CMPLX(0.0,0.0)
        x_seis_fre1(i,:)=CMPLX(0.0,0.0)
    endif
end do
call DFT(N, x_per_fre1, 1.0)
call DFT(N, x_seis_fre1, 1.0)
df = 1000./dt/N1

do i=1,N1
    f = (i-1)*df
    x_per_amp=sqrt(real(x_per_fre1(i,1))**2+aimag(x_per_fre1(i,2))**2)
    x_seis_amp=sqrt(real(x_seis_fre1(i,1))**2+aimag(x_seis_fre1(i,2))**2)
    write(4,100) i, f, x_per_amp, x_seis_amp
enddo

! TO 128 with zero
do i=1,N2
    if (i.le.N) then
        x_per_fre2(i,:)=CMPLX(x_per(i),0.0)
        x_seis_fre2(i,:)=CMPLX(x_seis(i),0.0)
    else
        x_per_fre2(i,:)=CMPLX(0.0,0.0)
        x_seis_fre2(i,:)=CMPLX(0.0,0.0)
    endif
end do
call DFT(N, x_per_fre2, 1.0)
call DFT(N, x_seis_fre2, 1.0)
df = 1000./dt/N2
do i=1,N2
    f = (i-1)*df
    x_per_amp=sqrt(real(x_per_fre2(i,1))**2+aimag(x_per_fre2(i,2))**2)
    x_seis_amp=sqrt(real(x_seis_fre2(i,1))**2+aimag(x_seis_fre2(i,2))**2)
    write(5,100) i, f, x_per_amp, x_seis_amp
end do

100 format(I4, 3(x,F15.4))

end program ch6_zeros


	SUBROUTINE FFT(LX,CX,SIGNI)
	COMPLEX CX(LX),CARG,CEXP,CW,CTEMP
	J=1
	SC=1.0/LX
	IF(SIGNI.EQ.1.0)SC=1.0
	SIG=-SIGNI
	DO 30 I=1,LX
	IF(I.GT.J)GO TO 10
	CTEMP=CX(J)*SC
	CX(J)=CX(I)*SC
	CX(I)=CTEMP
 10     M=LX/2
 20	IF(J.LE.M)GO TO 30
	J=J-M
	M=M/2
	IF(M.GE.1)GO TO 20
 30	J=J+M
	L=1
 40	ISTEP=2*L
	DO 50 M=1,L
	CARG=(0.0,1.0)*(3.14159265*SIG*(M-1))/L
	CW=CEXP(CARG)
	DO 50 I=M,LX,ISTEP
	CTEMP=CW*CX(I+L)
	CX(I+L)=CX(I)-CTEMP
 50	CX(I)=CX(I)+CTEMP
	L=ISTEP
	IF(L.LT.LX)GO TO 40
	RETURN
	END

SUBROUTINE DFT(N,C,S)
	COMPLEX C(N),Y(512),W
	IF(S.EQ.1.0)Z=1.0
	IF(S.EQ.-1.0)Z=1.0/FLOAT(N)
	DO 20 I=1,N
	Y(I)=(0.0,0.0)
	DO 20 J=1,N
	A=COS(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))
	B=SIN(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))*(-S)
	W=CMPLX(A,B)
 20	Y(I)=Y(I)+C(J)*W
	DO 30 I=1,N
 30	C(I)=Y(I)*Z
	RETURN
	END

load ch6_result_nozero.dat
load ch6_result_64.dat
load ch6_result_128.dat

figure(5)
hold on
subplot(3,2,1)
plot(ch6_result_nozero(:,2), ch6_result_nozero(:,3))
title('周期，无补零')
xlabel('frequency Hz')
ylabel('amplitude')

subplot(3,2,2)
plot(ch6_result_nozero(:,2), ch6_result_nozero(:,4))
title('非周期，无补零')
xlabel('frequency Hz')
ylabel('amplitude')


subplot(3,2,3)
plot(ch6_result_64(:,2), ch6_result_64(:,3))
title('周期，64个点')
xlabel('frequency Hz')
ylabel('amplitude')

subplot(3,2,4)
plot(ch6_result_64(:,2), ch6_result_64(:,4))
title('非周期，64个点')
xlabel('frequency Hz')
ylabel('amplitude')


subplot(3,2,5)
plot(ch6_result_128(:,2), ch6_result_128(:,3))
title('周期，128个点')
xlabel('frequency Hz')
ylabel('amplitude')

subplot(3,2,6)
plot(ch6_result_128(:,2), ch6_result_128(:,4))
title('非周期，128个点')
xlabel('frequency Hz')
ylabel('amplitude')

线性褶积与圆周褶积比较

program ch7_circle
real, parameter :: pi=3.1415926535
real, dimension(101) :: X=0., h_low=0., t
! real,dimension (101) :: h_low=0.
real,dimension (100+101-1) :: result_low=0., result_circle
real,dimension (201) :: X_c=0., h_r=0., t2, result_circle_r
real :: wc1=70.0, f ! 滤波器的滤波边界
integer :: i, dt=2
open(1, file='INPUT1.DAT', form='formatted', status='old')
open(2, file='ch7_result.dat',status='replace',action='write')
open(3, file='ch7_result_low.dat',status='replace',action='write')
open(4, file='ch7_result_circle.dat',status='replace',action='write')
read (1,*)(X(i), i=1,100)

! 理想低通滤波器

low_pass_fliter: do i=1,101
t(i) = dt*(i-100/2-1)/1000.
	if ((t(i).lt.0.0001) .and. (t(i).gt.-0.0001)) then
	h_low(i)=2*pi*wc1/pi
	else
	h_low(i)=sin(2*pi*wc1*t(i))/(pi*t(i))
	endif
end do low_pass_fliter

! 卷积
call CON(X,h_low,result_low,100,101,100+101-1)  ! result_low is the result of linear convolution
call CIRCON(X, h_low, result_circle, 101)

do i=1,101   ! 为圆周卷积补零
X_c(i) = X(i)
enddo

do i=1,201
t2(i) = dt*(i-200/2-1)/1000.
	if ((t2(i).lt.0.0001) .and. (t2(i).gt.-0.0001)) then
	h_r(i)=2*pi*wc1/pi
	else
	h_r(i)=sin(2*pi*wc1*t2(i))/(pi*t2(i))
	endif
end do
! 圆周卷积
call CIRCON(X_c, h_r, result_circle_r, 201)

do i=1,101
    write(2,300) i, (i-1.)*dt, X(i), result_low(i+50), result_circle(i)
end do
300 format(I4, 4(x, F15.4))


do i=1,101
    write(4,200) i, (i-1.)*dt, X_c(i), result_circle_r(i+100)
end do
200 format(I4, 3(x, F15.4))


end program ch7_circle



	SUBROUTINE CON(A,B,C,I,J,K)
	DIMENSION A(I),B(J),C(K)
	DO 1 K1=1,K
  1	C(K1)=0.0
	DO 2 I1=1,I
	DO 2 I2=1,J
	II=I1+I2-1
  2	C(II)=C(II)+A(I1)*B(I2)*0.004
	RETURN
	END


	SUBROUTINE CIRCON(X,H,Y,N)
	REAL X(N),H(N),Y(N)
	DO I=1,N
	Y(I)=0.0
	ENDDO
	DO I=1,N
	   DO J=1,N
	     IF(I-J.GT.0) Y(I)=Y(I)+X(J)*H(I-J)
	   ENDDO 
	ENDDO
	END

clear
load ch7_result.dat
load ch7_result_circle.dat
figure(5)
hold on
subplot(3,1,1)
plot(ch7_result(:,2), ch7_result(:,3))
title('源数据')
xlabel('T s')
ylabel('amplitude')

subplot(3,1,2)
plot(ch7_result(:,2), ch7_result(:,4))
title('线性褶积')
xlabel('T s')
ylabel('amplitude')


subplot(3,1,3)
plot(ch7_result(:,2), ch7_result(:,5))
title('圆周褶积')
xlabel('T s')
ylabel('amplitude')


figure(6)
hold on
subplot(3,1,1)
plot(ch7_result_circle(:,2), ch7_result_circle(:,3))
title('源数据')
xlabel('T s')
ylabel('amplitude')

subplot(3,1,2)
plot(ch7_result_circle(:,2), ch7_result_circle(:,4))
title('圆周褶积')
xlabel('T s')
ylabel('amplitude')

最小平方法求反射系数

program ch8_leasqu
implicit none
integer, parameter :: N=200, Nb=12, m=61
real, dimension(Nb) :: bn
real, dimension(N) :: ref_coe
real, dimension(m) :: rbb, a
real, dimension(N+Nb-1) :: X
real, dimension(N+Nb-1+m-1) :: ref_coe_inver
integer :: i
bn=[17.5, 12.5, 7.0, 0.0, -3.5, -7.0,-3.0, -1.0, 0.0, 1.4, 3.5, 2.0]
! 读取反射系数序列 ref_coe
open(1, file='INPUT8.DAT', form='formatted', status='old')
read(1,*)(ref_coe(i), i=1,200)
! 计算合成地震记录
call CON(bn, ref_coe, X, Nb, N, N+Nb-1)  

! 使用最小平方法求取反射系数序列
! 求子波自相关
call AUTCOR(Nb,m,bn,rbb)
! 解相应的托布利兹矩阵，求反子波a。
call PEO(m,rbb,a)
! 反子波与地震记录褶积，得到反射系数序列
call CON(a, X, ref_coe_inver, m, N+Nb-1, N+Nb-1+m-1)
open(3, file='ch8_input.dat', status='replace',action='write')
do i=1,N
    write(3,200) i, ref_coe(i)
enddo
200 format(I4, x, F15.4)
open(2, file='ch8_result.dat', status='replace',action='write')
do i=1,N+Nb-1+m-1
    write(2,100) i, ref_coe_inver(i)
enddo
100 format(I4, x, F15.6)

end program ch8_leasqu

	SUBROUTINE CON(A,B,C,I,J,K)
	DIMENSION A(I),B(J),C(K)
	DO 1 K1=1,K
  1	C(K1)=0.0
	DO 2 I1=1,I
	DO 2 I2=1,J
	II=I1+I2-1
  2	C(II)=C(II)+A(I1)*B(I2)*0.004
	RETURN
	END

! 解托布利兹矩阵，只有期望输出为脉冲时适用
	SUBROUTINE PEO(LR,R,A)
	DIMENSION R(LR),A(LR)
	V=R(1)
	A(1)=1.0
	IF(LR.EQ.1)RETURN
	DO 6 L=2,LR
	D=0.0
	L3=L-1
	DO 1 J=1,L3
	K=L-J+1
 1	D=D+A(J)*R(K)
	C=D/V
	IF(L.EQ.2)GOTO 4
	L1=(L-2)/2
	L2=L1+1
	IF(L2.LT.2)GOTO 3
	DO 2 J=2,L2
	HOLD=A(J)
	K=L-J+1
	A(J)=A(J)-C*A(K)
 2	A(K)=A(K)-C*HOLD
	IF(2*L1.EQ.L-2)GOTO 4
 3	LT3=L2+1
	A(LT3)=A(LT3)-C*A(LT3)
 4	A(L)=-C
 6	V=V-C*D
	RETURN
	END

! 自相关程序
! C自相关的数组 A输出
    SUBROUTINE AUTCOR(J,K,C,A)
	DIMENSION C(J),A(K)
	DO 7 N=1,K
	A(N)=0.0
	I=N
	DO 6 IN=I,J
	IL=IN-N+1
  6	A(N)=A(N)+C(IN)*C(IL)
  7	CONTINUE
	RETURN
	END

load ch8_input.dat
load ch8_result.dat

figure(1)
hold on
subplot(3,1,1)
plot(ch8_input(:,1), ch8_input(:,2))
title('原始反射系数序列')
subplot(3,1,2)
plot(ch8_result(:,1), ch8_result(:,2))
title('反演反射系数序列未截取')

subplot(3,1,3)
plot(ch8_result(1:200,1), ch8_result(1:200,2))
title('反演反射系数序列截取到200')

子波零相位变换

program ch9_zerophase
integer, parameter :: N=25
real, dimension(N) :: bo
complex, dimension(N,2) :: bo_comp


bo=[0.,1.,3.,5.,7.,9.,10.,8.,6.,4.,2.,0.,-1.,0.,1.,2.,3.,1.5,0.,-1.,0.,2.,1.,0.,0.]
do i=1,N
    bo_comp(i,:) = CMPLX(bo(i), 0.0)
enddo
! 对原始子波做傅氏变换到频率域
call DFT(N,bo_comp,1.0)
! 在频率域保持振幅谱不变，将虚部设置为0
do i=1,N
    bo_comp(i,:) = CMPLX(sqrt(real(bo_comp(i,1))**2+imag(bo_comp(i,2))**2), 0.0)
enddo
! 变换到时域
call DFT(N, bo_comp, -1.0)

open(1, file='ch9_result.dat', status='replace',action='write')
do i=1,N
    write(1,100) i, bo(i), real(bo_comp(i,1))
enddo
100 format (I4,2(x,F15.4))

end program ch9_zerophase


	SUBROUTINE DFT(N,C,S)
	COMPLEX C(N),Y(512),W
	IF(S.EQ.1.0)Z=1.0
	IF(S.EQ.-1.0)Z=1.0/FLOAT(N)
	DO 20 I=1,N
	Y(I)=(0.0,0.0)
	DO 20 J=1,N
	A=COS(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))
	B=SIN(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N))*(-S)
	W=CMPLX(A,B)
 20	Y(I)=Y(I)+C(J)*W
	DO 30 I=1,N
 30	C(I)=Y(I)*Z
	RETURN
	END

load ch9_result.dat

figure(1)
hold on
subplot(2,1,1)
plot(ch9_result(:,1), ch9_result(:,2))
title('原始子波')
xlabel('序号i，无意义')

subplot(2,1,2)
plot(ch9_result(:,1), ch9_result(:,3))
title('零相位子波')
xlabel('序号i，无意义')

子波最小相位转换

program ch10_leastphase
integer, parameter :: Nb=25, m=26
real, dimension(Nb) :: b
real, dimension(m) :: rbb, y, a0
real, dimension(Nb+m-1) :: by, g, grev
real :: b00,bysingle
integer :: i
real, dimension(2*Nb+m-2) :: b0

! 已知非最小相位子波
b=[0.,1.,2.,3.,1.5,0.,-1.,0.,1.,3.,5.,7.,9.,10.,8.,6.,4.,2.,0.,-1.,0.,2.,3.,2.,0.]
! 求子波自相关
call AUTCOR(Nb,m,b,rbb)
! 解矩阵，得到y(t)
call PEO(m,rbb,y)
!得到b0（0）,用b00表示
call CON(b,y,by,Nb,m,Nb+m-1)
bysingle=0.
do i = 1,Nb+m-1
    bysingle=bysingle+by(i)**2
enddo
b00=1./bysingle

! 得到 a0
do i=1,m
    a0(i) = b00*y(i)
enddo

! 得到g(t)
call CON(b,a0,g, Nb,m,Nb+m-1)

! 计算g(-t) 即grev
do i=1,Nb+m-1
    grev(i) = g(Nb+m-i)
enddo

! 计算b0
call CON(b,grev,b0,Nb,Nb+m-1,2*Nb+m-2)

open(2, file='ch10_result.dat', status='replace',action='write')
open(1, file='ch10_input.dat', status='replace',action='write')
do i=1,Nb
    write(1,100) i, b(i)
enddo
do i=1,2*Nb+m-2
    write(2,100) i, b0(i)
enddo
100 format(I4,x, F15.4)
end program ch10_leastphase


	SUBROUTINE CON(A,B,C,I,J,K)
	DIMENSION A(I),B(J),C(K)
	DO 1 K1=1,K
  1	C(K1)=0.0
	DO 2 I1=1,I
	DO 2 I2=1,J
	II=I1+I2-1
  2	C(II)=C(II)+A(I1)*B(I2)*0.004
	RETURN
	END

! 解托布利兹矩阵，只有期望输出为脉冲时适用
	SUBROUTINE PEO(LR,R,A)
	DIMENSION R(LR),A(LR)
	V=R(1)
	A(1)=1.0
	IF(LR.EQ.1)RETURN
	DO 6 L=2,LR
	D=0.0
	L3=L-1
	DO 1 J=1,L3
	K=L-J+1
 1	D=D+A(J)*R(K)
	C=D/V
	IF(L.EQ.2)GOTO 4
	L1=(L-2)/2
	L2=L1+1
	IF(L2.LT.2)GOTO 3
	DO 2 J=2,L2
	HOLD=A(J)
	K=L-J+1
	A(J)=A(J)-C*A(K)
 2	A(K)=A(K)-C*HOLD
	IF(2*L1.EQ.L-2)GOTO 4
 3	LT3=L2+1
	A(LT3)=A(LT3)-C*A(LT3)
 4	A(L)=-C
 6	V=V-C*D
	RETURN
	END

! 自相关程序
! C自相关的数组 A输出
    SUBROUTINE AUTCOR(J,K,C,A)
	DIMENSION C(J),A(K)
	DO 7 N=1,K
	A(N)=0.0
	I=N
	DO 6 IN=I,J
	IL=IN-N+1
  6	A(N)=A(N)+C(IN)*C(IL)
  7	CONTINUE
	RETURN
	END

load ch10_input.dat
load ch10_result.dat

figure(1)
hold on
subplot(2,1,1)
plot(ch10_input(:,1), ch10_input(:,2))
title('原始子波')
subplot(2,1,2)
plot(ch10_result(:,1), ch10_result(:,2))
title('最小相位子波未截取')

静校正

program ch11_static
real, dimension(1000) :: X
real, dimension(100) :: refer
real, dimension(10,100) :: Xorigin, X_static
integer, dimension(10,1) :: static
integer start, i, k, j, sta
integer, parameter :: M=50
real, dimension(10,M) :: rxy
open(1,file='SEIS.DAT', form='formatted',status='old')
read(1,*)(X(i),i=1,1000)
start=1
do j=1,10
    do k=1,100
        Xorigin(j,k)=X(start)
        start=start+1
    enddo
enddo

! 选择第1道作为参考道，形成参考道refer
do i=1,100
    refer(i) = Xorigin(1,i)
enddo

! refer与各道做互相关求取各道的相对剩余静校正量static
! 做互相关
do i=1,10
    call COAUTCOR(100,100,M,refer, Xorigin(i,:), rxy(i,:))
enddo

! 绘制互相关曲线
open(3, file='ch11_outTAO.dat', status='replace',action='write')

do j=1,M
    write(3,300)  j-1, rxy(1,j), rxy(2,j), rxy(3,j), rxy(4,j), rxy(5,j), rxy(6,j), rxy(7,j), rxy(8,j), rxy(9,j), rxy(10,j)
enddo
300 format(I4,10(x,F15.4))

! 由上面的互相关曲线得到了剩余静校正量：static
do i=1,10
    static(i,:)=maxloc(rxy(i,:))
enddo

! 进行静校正
! 初始化为0，避免生成错误数据
do i=1,10
    do j=1,100
        X_static(i,j)=0.
    enddo
enddo
! 进行静校正
do i=1,10
    sta=static(i,1)
    k=1
    do j=sta,100
        X_static(i,k)=Xorigin(i,j)
        k=k+1
    enddo
enddo
! 输出原始数据集
open(2, file='ch11_input.dat', status='replace',action='write')
do i=1,1000
    write(2,200)  i, X(i)
enddo
200 format(I4,x,F15.4)
! 输出静校正后的数据
open(4, file='ch11_outstatic.dat', status='replace',action='write')
do j=1,100
    write(4,400)  j, X_static(1,j), X_static(2,j), X_static(3,j), &
    X_static(4,j), X_static(5,j), X_static(6,j), X_static(7,j), X_static(8,j), X_static(9,j), X_static(10,j)
enddo
400 format(I4,10(x,F15.4))

end program ch11_static

! 互相关
	SUBROUTINE COAUTCOR(J,L,K,A,B,C)     
	DIMENSION A(J),B(L),C(K)
	DO 7 N=1,K
	C(N)=0.0
	IF (J.LE.L+1-N) THEN
		MIN=J
	ELSE 
		MIN=L+1-N
	ENDIF
 	DO 6 IA=1,MIN
	IB=IA+N-1
  6	C(N)=C(N)+A(IA)*B(IB)
  7	CONTINUE
	RETURN
	END

clear
load ch11_input.dat
load ch11_outTAO.dat
load ch11_outstatic.dat

figure(1)
hold on
sgtitle('原始地震道')
for i=1:10
    subplot(10,1,i)
    plot(1:100, ch11_input(i*100-99:i*100,2))
end

figure(2)
hold on
sgtitle('相关函数曲线')
for i=1:10 
    subplot(10,1,i)
    plot( ch11_outTAO(:,1), ch11_outTAO(:,i+1))
end

figure(3)
hold on
sgtitle('静校正后的地震数据')
for i=1:10 
    subplot(10,1,i)
    plot( ch11_outstatic(:,1), ch11_outstatic(:,i+1))
end