```from rsf.proj import * import math from rsf.recipes.beg import server as private nt=501 dt = 0.008 nf = 300 ot=0 df = 1/(nt*dt) wf = 2*math.pi def grey(title,bias,other): return ''' transp | grey title="%s" bias=%g color=j scalebar=n screenratio=0.8 barwidth=0.2 crowd1=0.75 crowd2=0.3 wherexlabel=b allpos=y label2=Time unit2=s label1=Frequency unit1=Hz %s ''' % (title,bias,other) ########################### # Model 0 Ricker wavelets ########################### ## Flow('ricker',None, ## ''' ## spike n1=512 d1=1 o1=0 k1=128,256 nsp=2 mag=1,0.5 | ## ricker1 frequency=0.2 ## ''') ## Flow('rtft rbasis','ricker', ## ''' ## rtft rect=5 nw=257 w0=0 dw=0.002 ## verb=y niter=100 basis=\${TARGETS[1]} ## ''') ## Result('rtft', ## ''' ## stack axis=3 norm=n | ## grey transp=n yreverse=n color=j ## title="LTF transform spectra (smooth=7)" ## crowd1=0.75 crowd2=0.25 ## labelfat=3 font=2 titlefat=3 parallel2=n ## ''') ## Flow('inv','rtft','rtft inv=y') ## Result('inv','inv ricker', ## ''' ## cat axis=3 \${SOURCES[1]} | ## graph title=Comparison ## ''') ## Flow('ltft0','ricker','ltft rect=7 verb=n') ## Result('ltft0', ## ''' ## window n1=%d | ## math output="abs(input)" | real | ## grey transp=n yreverse=n color=j ## title="LTF transform spectra (smooth=7)" ## crowd1=0.75 crowd2=0.25 ## labelfat=3 font=2 titlefat=3 parallel2=n ## ''' % nt) ########################### # Model 1 crossing chirps ########################### Flow('cchirps',None, ''' spike n1=%d d1=1 o1=0 | math output="cos(%g*(10+x1/7)*x1/%d)+cos(%g*(%d/2.8-x1/6.0)*x1/%d)" | pad end1=11 ''' % (nt,2*math.pi,nt,2*math.pi,nt,nt)) Result('cchirps', ''' window n1=%d | graph title=Input labelfat=3 font=2 titlefat=3 crowd1=0.75 crowd2=0.25 parallel2=n wherexlabel=t wheretitle=b plotfat=3 '''% nt) ## nf0 = 257 ## df = 0.00195312 ## wf = 2*math.pi ## #---cos basis ## Flow('cbasis','cchirps', ## ''' ## spray axis=2 n=%d d=%g o=0 | ## math output="cos(%g*x2*x1)" ## ''' % (nf0,df,wf)) ## Result('cbasis', ## ''' ## grey transp=n yreverse=n color=j ## title="Cosine bases" label2=Frequency unit2=Hz ## screenratio=0.8 crowd1=0.75 crowd2=0.3 ## labelfat=3 font=2 titlefat=3 parallel2=n ## ''') ## #---sin basis ## Flow('sbasis','cchirps', ## ''' ## spray axis=2 n=%d d=%g o=0 | math output="sin(%g*x2*x1)" ## ''' % (nf0,df,wf)) ## Result('sbasis', ## ''' ## grey transp=n yreverse=n color=j ## title="Sine bases" label2=Frequency unit2=Hz ## screenratio=0.8 crowd1=0.75 crowd2=0.3 ## labelfat=3 font=2 titlefat=3 parallel2=n ## ''') ## # 'window j2=4 | wiggle poly=y title="Sine bases" ') # S transform Flow('st1','cchirps','st') Result('st1', ''' window n1=%d | math output="abs(input)" | real | grey transp=n yreverse=n color=j title="S transform spectra" crowd1=0.75 crowd2=0.25 labelfat=3 font=2 titlefat=3 parallel2=n ''' % nt) ## Flow('ist1','st1','st inv=y') ## Result('ist1', ## ''' ## window n1=%d | ## graph title="Inverse S transform" label2= unit2= ## labelfat=3 font=2 titlefat=3 ## screenratio=0.8 crowd1=0.75 crowd2=0.3 parallel2=n ## '''% nt) ## Flow('diff2','cchirps ist1','add scale=1,-1 \${SOURCES[1]}') ## Result('diff2', ## ''' ## window n1=%d | ## graph title="Difference between input and IST" label2= unit2= ## labelfat=3 font=2 titlefat=3 ## screenratio=0.8 crowd1=0.75 crowd2=0.3 ## min2=-0.05 max2=0.05 parallel2=n ## ''' % nt) # Local prediction filter ## Flow('tfreq1','cchirps','timefreq rect=20') ## Result('tfreq1', ## ''' ## grey color=j wanttitle="Time-frequency spectra (lpf)" ## ''') # test sfltft Flow('ltft1 scbasis','cchirps','ltft rect=7 verb=n basis=\${TARGETS[1]}') ## Result('sbasis','scbasis', ## ''' ## window n2=257 | ## grey transp=n yreverse=n color=j ## title="Sine bases" label2=Frequency unit2=Hz ## crowd1=0.75 crowd2=0.2 ## labelfat=3 font=2 titlefat=3 parallel2=n ## ''') ## Result('cbasis','scbasis', ## ''' ## window f2=257 | put o2=0 | ## grey transp=n yreverse=n color=j ## title="Cosine bases" label2=Frequency unit2=Hz ## crowd1=0.75 crowd2=0.2 ## labelfat=3 font=2 titlefat=3 parallel2=n ## ''') Result('ltft1', ''' window n1=%d | math output="abs(input)" | real | grey transp=n yreverse=n color=j title="LTF decomposition spectra (smooth=7)" crowd1=0.75 crowd2=0.25 labelfat=3 font=2 titlefat=3 parallel2=n ''' % nt) ## Flow('iltft1','ltft1','ltft inv=y') ## Result('iltft1', ## ''' ## window n1=%d | ## graph title="Inverse LTF transform" label2= unit2= ## labelfat=3 font=2 titlefat=3 ## screenratio=0.8 crowd1=0.75 crowd2=0.3 ## ''' % nt) ## Flow('diff1','cchirps iltft1','add scale=1,-1 \${SOURCES[1]}') ## Result('diff1', ## ''' ## window n1=%d | ## graph title="Difference between input and ILTFT" label2= unit2= ## labelfat=3 font=2 titlefat=3 ## screenratio=0.8 crowd1=0.75 crowd2=0.3 ## min2=-0.05 max2=0.05 parallel2=n ## ''' % nt) # choose parameters Flow('ltft12','cchirps','ltft rect=14 verb=n') Result('ltft12', ''' window n1=%d | math output="abs(input)" | real | grey transp=n yreverse=n color=j title="LTF decomposition spectra (smooth=14)" crowd1=0.75 crowd2=0.25 labelfat=3 font=2 titlefat=3 parallel2=n ''' % nt) ## Flow('iltft12','ltft12','ltft inv=y') ## Result('iltft12', ## ''' ## window n1=%d | ## graph title="Inverse LTF transform" label2= unit2= ## labelfat=3 font=2 titlefat=3 ## screenratio=0.8 crowd1=0.75 crowd2=0.3 ## ''' % nt) ## Flow('diff12','cchirps iltft12','add scale=1,-1 \${SOURCES[1]}') ## Result('diff12', ## ''' ## window n1=%d | ## graph title="Difference between input and ILTFT" label2= unit2= ## labelfat=3 font=2 titlefat=3 ## screenratio=0.8 crowd1=0.75 crowd2=0.3 ## min2=-0.05 max2=0.05 parallel2=n ## ''' % nt) ## Flow('ltft13','cchirps','ltft rect=11 verb=y') ## Result('ltft13', ## ''' ## window n1=%d | ## math output="abs(input)" | real | ## grey transp=n yreverse=n color=j ## title="LTF transform spectra (weak)" ## crowd1=0.75 crowd2=0.25 ## labelfat=3 font=2 titlefat=3 parallel2=n ## ''' % nt) ## Flow('iltft13','ltft13','ltft inv=y') ## Result('iltft13', ## ''' ## window n1=%d | ## graph title="Inverse LTF transform" label2= unit2= ## labelfat=3 font=2 titlefat=3 ## screenratio=0.8 crowd1=0.75 crowd2=0.3 ## ''' % nt) ## Flow('diff13','cchirps iltft13','add scale=1,-1 \${SOURCES[1]}') ## Result('diff13', ## ''' ## window n1=%d | ## graph title="Difference between input and ILTFT" label2= unit2= ## labelfat=3 font=2 titlefat=3 ## screenratio=0.8 crowd1=0.75 crowd2=0.3 ## min2=-0.05 max2=0.05 parallel2=n ## ''' % nt) ## ########################### ## # Model 2 three cosines ## ########################### ## nlevel0=0 ## nf0=200 ## Flow('cosines',None, ## ''' ## math n1=%d o1=0 d1=%g n2=1 ## output="cos(%g*(10*x1))+cos(%g*(20*x1))+cos(%g*(30*x1))" | ## put label1=Time unit1=s label2=Amplitude ## ''' % (nt,dt,wf,wf,wf)) ## Flow('spikes',None, ## ''' ## spike n1=%d o1=0 d1=%g n2=1 nsp=2 mag=10,10 k1=%d,%d ## ''' % (nt,dt,nt/2,nt/2+40)) ## Flow('cspikes','cosines spikes','add \${SOURCES[1]} | noise var=%g' % (nlevel0)) ## Result('cspikes', ## ''' ## graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude ## label1=Time unit1=s crowd2=0.3 wanttitle=n ## ''') ## Flow('ltft2','cspikes','ltft rect=6 verb=n') ## Result('ltft2', ## ''' ## math output="abs(input)" | real | ## '''+ grey('LTFT',0.01,'yreverse=n wheretitle=t \ ## labelfat=3 font=2 titlefat=3')) ## Flow('iltft2','ltft2','ltft inv=y') ## Result('iltft2', ## ''' ## graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude ## label1=Time unit1=s crowd2=0.3 title="Inverse LTF transform" ## ''') ## # Guochang's method ## #---cos basis ## Flow('bcos-0','cspikes', ## 'spray axis=2 n=%d d=%g o=0| math output="cos(%g*x2*x1)"' % (nf0,df,wf)) ## #---sin basis ## Flow('bsin1-0','cspikes', ## ''' ## spray axis=2 n=%d d=%g o=0| math output="sin(%g*x2*x1)"| window f2=1 ## ''' % (nf0,df,wf)) ## #-----similarity vs projection ## Flow('mask','cspikes', ## 'noise rep=y type=n seed=2008 | mask min=0') ## #Flow('zero','s-0 mask','headercut mask=\${SOURCES[1]}') ## Flow('bsincos-0','bsin1-0 bcos-0','cat axis=2 \${SOURCES[1]}') ## Flow('s-0s2','cspikes mask', ## ''' ## spray axis=2 n=%d d=%g o=0 | transp | ## headercut mask=\${SOURCES[1]} | transp ## ''' % (2*nf0-1, df)) ## Result('s-0s2', ## ''' ## window n2=1 | ## graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude ## label1=Time unit1=s crowd2=0.3 wanttitle=n ## ''') ## Flow('project-0 pred','bsincos-0 s-0s2', ## 'lpf match=\${SOURCES[1]} niter=100 rect1=20 pred=\${TARGETS[1]}') ## Flow('proj-0-old','project-0', ## ''' ## pad beg2=1 | put n2=%d n3=2 | math output=input*input | ## stack axis=3 norm=n | math output="sqrt(input)" ## ''' % nf0) ## Result('proj-0-old',grey('',0.2,'')) ## Flow('proj-0','cspikes','timefreq nw=200 dw=0.249501 rect=15') ## Result('proj-0',grey('',0.2,'')) ## Result('pred', ## ''' ## stack | ## graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude ## label1=Time unit1=s crowd2=0.3 wanttitle=n ## ''') ## ########################### ## # Model 3 frequency-vary sinusoid ## ########################### ## Flow('simple',None, ## ''' ## math n1=251 o1=0 d1=0.004 n2=1 output="cos(%g*(3.0*x1*x1+5)*x1)" ## label1=Time unit1=s ## ''' % (2*wf)) ## Result('simple', ## ''' ## graph yreverse=y transp=y wanttitle=n ## screenratio=2.5 labelsz=6 ## ''') ## Flow('ltft3','simple','ltft rect=6 verb=n') ## Result('ltft3', ## ''' ## math output="abs(input)" | real | ## grey yreverse=y transp=y wanttitle=n color=j ## screenratio=2.5 labelsz=6 ## ''') ## Flow('iltft3','ltft3','ltft inv=y') ## Result('iltft3', ## ''' ## graph yreverse=y transp=y title=Inverse ## screenratio=2.5 labelsz=6 ## ''') ## #---cos basis ## Flow('bcos1','simple', ## ''' ## spray axis=2 n=20 d=2 o=0 label=Frequency unit=Hz | ## math output="cos(%g*x2*x1)" | window f2=1 ## ''' % wf) ## Flow('bcos','simple bcos1', ## ''' ## spray axis=2 n=1 d=2 o=0 label=Frequency unit=Hz | ## math output=0.5 | cat axis=2 \${SOURCES[1]} ## ''') ## Result('bcos', ## ''' ## wiggle yreverse=y transp=y wanttitle=n ## poly=y screenratio=2 labelsz=6 ## ''') ## #---sin basis ## Flow('bsin','simple', ## ''' ## spray axis=2 n=20 d=2 o=0 label=Frequency unit=Hz | ## math output="sin(%g*x2*x1)"| window f2=1 ## ''' % wf) ## Result('bsin', ## ''' ## wiggle yreverse=y transp=y wanttitle=n poly=y ## screenratio=2 wantaxis1=n labelsz=6 ## ''') ## ########################### ## # Model 4 chirp signal: ## # cos(g(f)*t+a)+cos(g(f)*t+b) ## ########################### ## nlevel1=0 ## dt1 = 0.004 ## Flow('tchirps',None, ## ''' ## math n1=%d o1=0 d1=%g n2=1 ## output="cos(%g*(3.0*x1*x1+5)*x1)+cos(%g*(8*x1*x1+19)*x1)" ## ''' % (nt, dt1, wf, wf)) ## Result('tchirps', ## ''' ## graph pad=n screenratio=0.8 crowd1=0.75 ## label2=Amplitude label1=Time unit1=s crowd2=0.3 wanttitle=n ## ''') ## Flow('ltft4','tchirps','ltft rect=6 verb=n') ## Result('ltft4', ## ''' ## math output="abs(input)" | real | ## ''' + grey('',0.02,'')) ## Flow('iltft4','ltft4','ltft inv=y') ## Result('iltft4', ## ''' ## put unit2= | ## graph pad=n screenratio=0.8 crowd1=0.75 ## label2=Amplitude label1=Time unit1=s crowd2=0.3 title=Inverse ## ''') ## # S transform ## Flow('st2','tchirps','pad end1=11 | st') ## Result('st2', ## ''' ## window n1=501 | ## math output="abs(input)" | real | ## ''' + grey('S transform',0.02,'wheretitle=t')) ## Flow('ist2','st2','st inv=y') ## Result('ist2', ## ''' ## window n1=501 | ## graph pad=n screenratio=0.8 crowd1=0.75 ## label2=Amplitude unit2= label1=Time unit1=s crowd2=0.3 ## title="Inverse S transform" ## ''') ## # Guochang's method ## #---cos basis ## Flow('bcos-1','tchirps', ## ''' ## spray axis=2 n=%d d=%g o=0 | ## math output="cos(%g*x2*x1)" ## ''' % (nf,df,wf)) ## #---sin basis ## Flow('bsin1-1','tchirps', ## ''' ## spray axis=2 n=%d d=%g o=0 | ## math output="sin(%g*x2*x1)"| ## window f2=1 ## ''' % (nf,df,wf)) ## #-----similarity vs projection ## Flow('bsincos-1','bsin1-1 bcos-1','cat axis=2 \${SOURCES[1]}') ## Flow('s-1s2','tchirps','spray axis=2 n=%d d=%g o=0' % (2*nf-1, df)) ## Flow('project-1','bsincos-1 s-1s2', ## 'lpf match=\${SOURCES[1]} niter=30 rect1=23 ') ## Flow('proj-1-old','project-1', ## ''' ## pad beg2=1 | put n2=%d n3=2 | ## math output=input*input | ## stack axis=3 norm=n | ## math output="sqrt(input)" ## ''' % nf) ## Flow('proj-1','tchirps','timefreq nw=300 dw=0.249501 rect=18') ## Result('proj-1',grey('',0.2,'')) ## ########################### ## # Model 5 synthetic trace ## # by ricker wavelet ## ########################### ## nlevel2=0 ## dt2=0.002 ## nf2=500 ## Flow('syn1',None, ## ''' ## spike n1=%d d1=%g o1=0 nsp=4 k1=%g,%g,%g,%g mag=1,1,1,1 | ## ricker1 frequency=50 ## ''' % (nt,dt2,0.2*nt,0.4*nt,0.605*nt,0.8*nt)) ## Flow('syn2',None, ## ''' ## spike n1=%d d1=%g o1=0 nsp=2 k1=%g,%g mag=1,1 | ## ricker1 frequency=40 ## ''' % (nt,dt2,0.6*nt,0.79*nt)) ## Flow('syn3',None, ## ''' ## spike n1=%d d1=%g o1=0 nsp=3 k1=%g,%g,%g mag=1,1,1 | ## ricker1 frequency=30 ## ''' % (nt,dt2,0.21*nt,0.59*nt,0.81*nt)) ## Flow('syn4',None, ## ''' ## spike n1=%d d1=%g o1=0 nsp=1 k1=%g mag=2 | ## ricker1 frequency=20 ## ''' % (nt,dt2,0.43*nt)) ## Flow('syn','syn1 syn2 syn3 syn4', ## ''' ## add \${SOURCES[1:4]} | ## noise var=%g ## ''' % (nlevel2)) ## Flow('secs','syn syn1 syn2 syn3 syn4', ## ''' ## cat \${SOURCES[1:5]} axis=2 ## ''') ## Result('secs', ## ''' ## dots gaineach=n ## labels="Sum\^ \_\^ \_=:50Hz\^ ## \_\^ \_+:40Hz\^ \_\^ \_+:30Hz\^ \_\^ \_+:20Hz" ## wanttitle=n titles= ## ''') ## Result('syn', ## ''' ## graph pad=n screenratio=0.8 crowd1=0.75 ## label2=Amplitude label1=Time unit1=s ## crowd2=0.3 wanttitle=n pad=n ## ''') ## Flow('ltft5','syn','ltft rect=6 verb=n nw=500 dw=0.249501') ## Result('ltft5', ## ''' ## math output="abs(input)" | real | ## ''' + grey('',0.0,'')) ## Flow('iltft5','ltft5','ltft inv=y') ## Result('iltft5', ## ''' ## put unit2= | ## graph pad=n screenratio=0.8 crowd1=0.75 ## label2=Amplitude label1=Time unit1=s ## crowd2=0.3 title=Inverse pad=n ## ''') ## # S transform ## Flow('st3','syn','pad end1=11 | st fhi=124.7505') ## Result('st3', ## ''' ## math output="abs(input)" | real | ## ''' + grey('S transform',0.0,'wheretitle=t')) ## Flow('ist3','st3','st inv=y') ## Result('ist3', ## ''' ## window n1=501 | put unit2= | ## graph pad=n screenratio=0.8 crowd1=0.75 ## label2=Amplitude label1=Time unit1=s ## crowd2=0.3 title="Inverse ST" pad=n ## ''') ## #---STFT with 30ms and 100ms Hamming window ## for case in ('s-2sh','s-2lg'): ## Fetch(case+'.rsf','guochang',private) ## Result(case,grey('',0,'') ) ## # Guochang's method ## #---cos basis ## Flow('bcos-2','syn', ## ''' ## spray axis=2 n=%d d=%g o=0 | ## math output="cos(%g*x2*x1)" ## ''' % (nf2,df,wf)) ## #---sin basis ## Flow('bsin1-2','syn', ## ''' ## spray axis=2 n=%d d=%g o=0 | ## math output="sin(%g*x2*x1)"| ## window f2=1 ## ''' % (nf2,df,wf)) ## #-----similarity vs projection ## Flow('bsincos-2','bsin1-2 bcos-2','cat axis=2 \${SOURCES[1]}') ## Flow('s-2s2','syn','spray axis=2 n=%d d=%g o=0' % (2*nf2-1, df)) ## Flow('project-2','bsincos-2 s-2s2','lpf match=\${SOURCES[1]} niter=80 rect1=8') ## Flow('proj-2-old','project-2', ## ''' ## pad beg2=1 | put n2=%d n3=2 | ## math output=input*input | ## stack axis=3 norm=n | ## math output="sqrt(input)" ## ''' % nf2) ## Flow('proj-2','syn','timefreq nw=500 dw=0.249501 rect=8 niter=80') ## Result('proj-2',grey('',0.,'')) ## ########################### ## # Model 6 random reflectivity ## # with ricker frequency f=at+b ## ########################### ## nf6=600 ## nt6=2000 ## dt6=0.001 ## nlevel6=0 ## Flow('refl',None, ## ''' ## math n1=%d d1=%g o1=%g output=0 | noise var=1 ## ''' % (nt6,dt6,ot)) ## Flow('rfreq','refl','math output="25*x1*x1+15"') ## Flow('rphase','refl','math output=0') ## Flow('random','refl rfreq rphase', ## ''' ## ricker2 tfreq=\${SOURCES[1]} tphase=\${SOURCES[2]} ## norm=n hiborder=100 | ## noise var=%g | ## agc rect1=100 ## ''' % nlevel6) ## Result('refl', ## ''' ## graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude ## label1=Time unit1=s crowd2=0.3 wanttitle=n ## ''') ## Result('random', ## ''' ## graph pad=n screenratio=0.8 crowd1=0.75 ## label2=Amplitude label1=Time unit1=s crowd2=0.3 wanttitle=n ## ''') ## Flow('ltft6','random','ltft rect=20 verb=n nw=800 dw=0.249501') ## Result('ltft6', ## ''' ## math output="abs(input)" | real | ## ''' + grey('',0.0,'')) ## Flow('iltft6','ltft6','ltft inv=y') ## Result('iltft6', ## ''' ## put unit2= | ## graph pad=n screenratio=0.8 crowd1=0.75 ## label2=Amplitude label1=Time unit1=s crowd2=0.3 title=Inverse ## ''') ## # Guochang's method ## #---cos basis ## Flow('bcos-6','random', ## ''' ## spray axis=2 n=%d d=%g o=0 | ## math output="cos(%g*x2*x1)" ## ''' % (nf6,df,wf)) ## #---sin basis ## Flow('bsin1-6','random', ## ''' ## spray axis=2 n=%d d=%g o=0 | ## math output="sin(%g*x2*x1)"| ## window f2=1 ## ''' % (nf6,df,wf)) ## Flow('bsincos-6','bsin1-6 bcos-6','cat axis=2 \${SOURCES[1]}') ## Flow('s-6s','random','spray axis=2 n=%d d=%g o=0' % (2*nf6-1, df)) ## Flow('project-6','bsincos-6 s-6s', ## 'lpf match=\${SOURCES[1]} niter=100 rect1=90 ') ## Flow('proj-6-old','project-6', ## ''' ## pad beg2=1 | put n2=%d n3=2 | ## math output=input*input | ## stack axis=3 norm=n | ## math output="sqrt(input)" ## ''' % nf6) ## Flow('proj-6','random','timefreq nw=800 dw=0.249501 rect=90') ## Result('proj-6','window max2=149 | ' + grey('',0.2,'')) ## #---test local frequency ## Flow('numer-6','proj-6',' math output="input*x2" | stack axis=2') ## Flow('denom-6','proj-6','stack axis=2') ## Flow('lcf-6','numer-6 denom-6','math x=\${SOURCES[1]} output="input/(x+0.01)"') ## Flow('theo-6','lcf-6','math output="25*x1*x1+15" ') ## Flow('flcf-6','random','iphase rect1=90 order=10 hertz=y') ## Result('lcf-6','lcf-6 theo-6 flcf-6', ## ''' ## cat axis=2 \${SOURCES[1]} \${SOURCES[2]}| ## graph pad=n screenratio=0.8 yreverse=y min2=0 max2=149 ## label2=Frequency unit2=Hz plotcol=6,5,3 ## crowd1=0.75 label1=Time unit1=s crowd2=0.3 wanttitle=n Xdash=0,5,1 ## ''') ## ########################### ## # Model 7 The whole 4D data ## ########################### ## nf9=601 ## df9=0.25 ## #-----read data from segy file ## Fetch('ln472_old.sgy','guochang',private) ## Flow('old','ln472_old.sgy', ## ''' ## segyread bfile=/dev/null hfile=/dev/null read=data | ## put n1=751 n2=201 o1=0 d1=0.002 02=1 d2=0.01 label1="Time" unit1="s" ## label2="X" unit2="km" | window min1=0 ## ''') ## Result('old','old', ## 'grey wanttitle=n label2="Distance" label1=Time unit1=s') ## Flow('rltft','old','ltft rect=5 verb=y nw=601 dw=0.25 niter=50') ## Result('rltft', ## ''' ## math output="abs(input)" | real | ## transp plane=23 | ## byte allpos=y gainpanel=40 pclip=100 | ## grey3 color=j frame1=615 frame3=40 frame2=76 label1=Time flat=y ## unit1=s label2=Distance label3=Frequency unit3=Hz ## point1=0.8 point2=0.8 wanttitle=n ## ''') ## Flow('riltft','rltft','ltft inv=y verb=y') ## Result('riltft', ## ''' ## grey title="Inverse LTF" label2="Distance" ## unit2=km label1=Time unit1=s ## ''') ## # Dominant frequency ## Flow('numer','rltft', ## ''' ## math output="abs(input)" | real | ## math output="input*x2"| stack ## ''') ## Flow('denom','rltft', ## ''' ## math output="abs(input)" | real | ## stack ## ''') ## Flow('lf','numer denom', ## 'math den=\${SOURCES[1]} output="input/(den+0.001)"') ## Result('lf', ## ''' ## grey allpos=y color=j label1=Time scalebar=y ## barlabel=Frequency barunit=Hz pclip=99. ## unit1=s label3=Frequency label2=Distance unit3=Hz wanttitle=n ## ''' ) ## #---cos basis ## Flow('bcos-9','old', ## ''' ## spray axis=3 n=%d d=%g o=0 | ## math output="cos(%g*x3*x1)" ## ''' % (nf9,df9,wf)) ## #---sin basis ## Flow('bsin1-9','old', ## ''' ## spray axis=3 n=%d d=%g o=0 | ## math output="sin(%g*x3*x1)"| ## window f3=1 ## ''' % (nf9,df9,wf)) ## #-----projection ## Flow('bsincos-9','bsin1-9 bcos-9','cat axis=3 \${SOURCES[1]}') ## Flow('olds','old','spray axis=3 n=%d d=%g o=0' % (2*nf9-1, df9)) ## #Flow('allproject-9','olds bsincos-9', ## # 'lpf match=\${SOURCES[1]} niter=30 rect1=14 ') ## win9 = [] ## for f9 in range(0,1200,100): ## winold9 = 'winold-%d' % f9 ## winbsincos9 = 'winsincos-9-%d' %f9 ## winproject9 = 'winproject-9-%d' % f9 ## win9.append(winproject9) ## Flow(winbsincos9,'bsincos-9','window f3=%d n3=100' % f9) ## Flow(winold9,'olds','window f3=%d n3=100' % f9) ## Flow(winproject9,[winbsincos9, winold9], ## 'lpf match=\${SOURCES[1]} niter=30 rect1=14') ## Flow('project-9',win9, ## ''' ## cat \${SOURCES[1:%d]} axis=3 | put o3=0 d3=%g label3=Frequency Unit3=Hz | ## math output="input*input" ## ''' % (len(win9),df9)) ## Flow('proj-9-old','project-9', ## ''' ## pad beg3=1 |pad end3=1 | put n3=%d n4=2 | math output=input*input | ## stack axis=4 norm=n | math output="sqrt(input)" ## ''' % nf9) ## Flow('proj-9','old','timefreq rect=14 nw=601 dw=0.25') ## Result('proj-9', ## ''' ## transp plane=23 | ## byte allpos=y gainpanel=40 pclip=100 | ## grey3 color=j frame1=615 frame3=40 frame2=76 label1=Time flat=y ## unit1=s label2=Distance label3=Frequency unit3=Hz ## point1=0.8 point2=0.8 wanttitle=n ## ''') ## #---test local frequency ## Flow('sproj-9','proj-9','smooth rect1=10 rect2=10 ') ## Flow('numer-9','proj-9','math output="input*x2" | stack') ## Flow('denom-9','proj-9','stack axis=2') ## Flow('glcf-9','numer-9 denom-9','add mode=d \${SOURCES[1]}') ## Flow('lcf-9','numer-9 denom-9','math x=\${SOURCES[1]} output="input/(x+0.001)"') ## Result('lcf-9', ## ''' ## grey allpos=y color=j label1=Time scalebar=y ## barlabel=Frequency barunit=Hz pclip=99. ## unit1=s label3=Frequency label2=Distance unit3=Hz wanttitle=n ## ''' ) ## ############################### End()```