Finding the shaping filter

The program shp() is used to find the shaping filter that reconstructs the HF from the LF for each interpolated trace. The routine is a modification of shaper() Claerbout (1992a) with contrunc() replacing contran().

# Shaping Filter (Modification of shaper [PVI])
# find the shaping filter for the cubed data in time domain
#
#  rr = yy - xx * ff
#%
subroutine  shp(n1,n2,nf,ny,yh,yl,x)
#
integer n1,n2,i1,i2,nf,ny,lag,iter,niter,i,nx
real yh(n1,n2),yl(n1,n2),x(nf,n2)
temporary real ff(nf),xx(n1),yy(ny),rr(ny)
temporary real res(ny,n2),df(nf),dr(ny),sf(nf),sr(ny)
#
lag=nf/2
niter=nf
nx=n1
do i2=1,n2 {
	call zero(n1,xx)
        call zero(ny,yy)
        call zero(ny,rr)
        call zero(nf,ff)
	do i1=1,n1 {
		xx(i1)=yl(i1,i2)**3
		yy(i1)=yh(i1,i2)
		}
do i=1,ny
	rr(i)=yy(i)
do iter= 0, niter {
        call contrunc( 1, 0,lag,nx,xx,  nf,df,  ny,rr)  # df=xx*rr
        call contrunc( 0, 0,lag,nx,xx,  nf,df,  ny,dr)  # dr=xx*df
        call cgstep( iter, nf,ff,df,sf,  ny,rr,dr,sr)   # rr=rr-dr; ff=ff+df
        }
#}
#do i2=1,n2 {
	do i1=1,nf {
		x(i1,i2)=ff(i1)
		}
}
call slice('filterxc.h',4,nf,1,ff) 
call slice('filter.h',4,nf,n2,x) 
return;end