The program Pstri() performs the interpolation on small windows of data which are then patched back together. The shaping filter is found by shpsafe() which is similar to shp(). The actual interpolation is performed by sip(). The decimated crosswell data set has dimensions n1=512, n2=158. The interpolation in this paper is performed by setting w1=32, w2=16 and k1=32, k2=64. This makes 2048 windows. The shaping filter is 32 points long (same length as the windows).
# Pstri: calculate shaping filter in windows and patch back together. # n1, w1, k1 (data dimension, window dimension, # of windows) # n2, w2, j2 #% end of self doc. integer n1,n2,nf,ntr, w1,w2, k1,k2 real d2 from history: integer n1, n2 from history: real d2 from par: integer ntr=1 from par: integer nf=32 from par: integer w1=32 , w2=16 from par: integer k1=64 , k2=128 allocate: real all(n1,n2), low(n1,n2) allocate: real shapeout(k1*nf,k2*w2), shape(nf,w2) allocate: real tempall(w1,w2), templow(w1,w2) allocate: real temptest(w1,(ntr+1)*w2),test(n1,(ntr+1)*n2) allocate: real final(n1,(ntr+1)*n2) allocate: real windwt(w1,(ntr+1)*w2),wallwt(n1,(ntr+1)*n2)subroutine Pstri(n1,n2,d2,nf,ntr,all,low,shapeout,tempall,templow, w1,w2, k1,k2,shape,test,temptest,windwt,wallwt,final) integer nf,ny,ntr, i1,n1, i2,n2, w1,j1,w2,j2, k1,k2, outfd,output, np real all(n1,n2), low(n1,n2),d2 real shapeout(k1*nf,k2*w2), shape(nf,w2) real tempall(w1,w2), templow(w1,w2) real temptest(w1,(ntr+1)*w2),test(n1,(ntr+1)*n2) real final(n1,(ntr+1)*n2) real windwt(w1,(ntr+1)*w2),wallwt(n1,(ntr+1)*n2) ny=w1+nf-1 np=(ntr+1)*w2 outfd=output() call putch('n1','i',n1) call putch('n2','i',(ntr+1)*n2) call putch('d2','f',d2/(ntr+1)) call putch('title','s','Interpolation done in patches') call hclose()
call sreed('in', all, 4*n1*n2 ) call sreed('inlow', low, 4*n1*n2 )
call cinloiwt(1,1,1,1,windwt,w1,np) call mkwallwt(k1,k2, windwt, w1,np, wallwt,n1,(ntr+1)*n2)
do i1=1,k1{ do i2=1,k2{ call patch( 0,0, i1,i2, k1,k2, all,n1,n2, tempall,w1,w2) call patch( 0,0, i1,i2, k1,k2, low,n1,n2, templow,w1,w2) call shpsafe(w1,w2,nf,ny,tempall,templow,shape) call patch( 1,1, i1,i2, k1,k2, shapeout,k1*nf,k2*w2, shape,nf,w2) call sip(w1,w2,tempall,templow,np,shape,ntr,nf,ny,temptest) call diag(0,0, windwt,w1*np,temptest,temptest) call patch( 1,1, i1,i2, k1,k2, test,n1,(ntr+1)*n2, temptest,w1,np) }} call diag(0,1,wallwt,n1*(ntr+1)*n2,test,final)
call srite('out', final, 4*n1*(ntr+1)*n2) return; end