```from rsf.proj import * # Download input trace trace = 'trace_il1190_xl1155.trace' Fetch(trace,'data', server='https://raw.githubusercontent.com', top='seg/tutorials/master/1410_Phase') # Convert format to RSF Flow('trace',trace, ''' echo n1=1503 in=\$SOURCE data_format=ascii_float o1=0 d1=0.004 label1=Time unit1=s | dd form=native | window f1=2 ''') # Apply Hilbert transform Flow('hilb','trace','envelope hilb=y') # Compute envelope and phase Flow('envp','trace hilb','math h=\${SOURCES[1]} output="sqrt(input*input+h*h)" ') Flow('envm','envp','scale dscale=-1') Flow('phase','trace hilb','math h=\${SOURCES[1]} output="h&input" ') Plot('trace','trace hilb envp envm', ''' cat \${SOURCES[1:4]} axis=2 | window f1=500 n1=201 | graph title="Trace Amplitudes" dash=0,1,0,0 plotcol=5,5,4,4 ''') Plot('phase', ''' window f1=500 n1=201 | graph title="Instantaneous Phase" plotcol=3 label2=Degree unit2=radians ''') Result('trace','trace phase','OverUnderAniso') # Repeat on 2-D section sgy = 'penobscot_xl1155.sgy' Fetch(sgy,'data', server='https://raw.githubusercontent.com', top='seg/tutorials/master/1410_Phase') Flow('data tdata',sgy,'segyread tfile=\${TARGETS[1]}') Flow('hdata','data','envelope hilb=y') Flow('edata','data hdata','math h=\${SOURCES[1]} output="sqrt(input*input+h*h)" ') Flow('pdata','data hdata','math h=\${SOURCES[1]} output="h&input" ') Plot('data','window min1=2 max1=3 | grey title="Original" scalebar=y') Plot('hdata','window min1=2 max1=3 | grey title="Hilbert" scalebar=y') Plot('edata','window min1=2 max1=3 | grey title="Envelope" scalebar=y color=j allpos=y') Plot('pdata','window min1=2 max1=3 | grey title="Phase" scalebar=y color=j') Result('data','data hdata edata pdata','TwoRows') # Plot envelope peaks Plot('env','trace envp envm', ''' cat \${SOURCES[1:3]} axis=2 | window f1=500 n1=201 | graph title="Envelope" plotcol=5,4,4 grid=y min1=2 max1=2.8 min2=-7200 max2=7200 ''') Flow('peaks','envp','max1') Flow('meaks','peaks','math output="conj(input)" ') Plot('peaks','peaks meaks', ''' cat axis=1 \${SOURCES[1]} | graph wanttitle=n wantaxis=n symbol=o plotcol=3 symbolsz=6 min1=2 max1=2.8 min2=-7200 max2=7200 ''') Plot('peaks2','peaks meaks', ''' cat axis=2 \${SOURCES[1]} | transp | graph wanttitle=n wantaxis=n plotcol=3 min1=2 max1=2.8 min2=-7200 max2=7200 ''') Plot('envpeaks','env peaks peaks2','Overlay') Flow('xpeaks','peaks','real') Flow('phasepeaks','phase xpeaks','inttest1 coord=\${SOURCES[1]} interp=spline nw=4') Flow('cpeaks','xpeaks','rtoc') Plot('ppeaks','xpeaks phasepeaks cpeaks', ''' cmplx \${SOURCES[1]} | cat axis=2 \${SOURCES[2]} | transp | graph wanttitle=n wantaxis=n plotcol=3 min1=2 max1=2.8 min2=-3.5 max2=3.5 ''') Plot('ppeaks2','xpeaks phasepeaks', ''' cmplx \${SOURCES[1]} | graph title=Phase symbol=o plotcol=3 symbolsz=6 grid=y min1=2 max1=2.8 min2=-3.5 max2=3.5 ''') # ideal phase from math import pi Flow('iphase0','phasepeaks', ''' mask min=%g max=%g | dd type=float | math p=\$SOURCE output="(1-input)*p" ''' % (-pi/2,pi/2)) Flow('iphase1','iphase0', ''' mask min=%g max=%g | dd type=float | math p=\$SOURCE output="%g*input+(1-input)*p" ''' % (-3*pi/2,-pi/2,-pi)) Flow('iphase2','iphase1', ''' mask min=%g max=%g | dd type=float | math p=\$SOURCE output="%g*input+(1-input)*p" ''' % (pi/2,3*pi/2,pi)) Plot('ipeaks','xpeaks iphase2 cpeaks', ''' cmplx \${SOURCES[1]} | cat axis=2 \${SOURCES[2]} | transp | graph wanttitle=n wantaxis=n plotcol=2 plotfat=3 min1=2 max1=2.8 min2=-3.5 max2=3.5 ''') Plot('ipeaks2','xpeaks iphase2', ''' cmplx \${SOURCES[1]} | graph title=Phase symbol=o plotcol=2 symbolsz=8 grid=y min1=2 max1=2.8 min2=-3.5 max2=3.5 plotfat=3 ''') Plot('phasepeaks','ipeaks ipeaks2 ppeaks ppeaks2','Overlay') Result('phasepeaks','envpeaks phasepeaks','OverUnderAniso') # Compute phase error for seismic section Flow('dpeaks','edata','max1 | real') Flow('ppeaks','pdata dpeaks','inttest1 coord=\${SOURCES[1]} interp=spline nw=4 same=n') Flow('ipeaks0','ppeaks', ''' mask min=%g max=%g | dd type=float | math p=\$SOURCE output="(1-input)*p" ''' % (-pi/2,pi/2)) Flow('ipeaks1','ipeaks0', ''' mask min=%g max=%g | dd type=float | math p=\$SOURCE output="%g*input+(1-input)*p" ''' % (-3*pi/2,-pi/2,-pi)) Flow('ipeaks2','ipeaks1', ''' mask min=%g max=%g | dd type=float | math p=\$SOURCE output="%g*input+(1-input)*p" ''' % (pi/2,3*pi/2,pi)) Flow('epeaks bar','ipeaks2 ppeaks', 'math p=\${SOURCES[1]} output="abs(p-input)" | byte allpos=y bar=\${TARGETS[1]} pclip=100') Result('epeaks','dpeaks epeaks', ''' rtoc | math output="input+I*x2" | graph symbol='.' min1=2 max1=3 min2=0 max2=600 title="Absolute Phase Error at Envelope Peaks" transp=y yreverse=y depth=\${SOURCES[1]} color=j scalebar=y ''') # Zero-phase correction? Flow('skew','trace', 'pad beg1=250 end1=250 | localskew rect=250 a0=-90 da=3 na=121 | window n1=1501 f1=250') Plot('skew', ''' grey color=j allpos=y title="Inverse Local Skewness" label2="Phase Rotation" Xunit2="\^o\_" ''') # Pick & Slice Flow('pick','skew','pick vel0=-45 rect1=50') Plot('pick','graph plotcol=7 plotfat=5 transp=y yreverse=y min2=-90 max2=270 wanttitle=n pad=n wantaxis=n') Result('skew','skew pick','Overlay') # Phase rotation Flow('trace90','trace pick','phaserot phase=\${SOURCES[1]}') Flow('trace0','trace90','envelope order=100 hilb=y phase=-90') Result('trace2','trace trace0', ''' cat axis=2 \${SOURCES[1]} | dots labels="Input:Output" yreverse=y ''') End()```

 sfdd sfwindow sfenvelope sfmath sfscale sfcat sfgraph sfsegyread sfgrey sfmax1 sftransp sfreal sfinttest1 sfrtoc sfcmplx sfmask sfbyte sfpad sflocalskew sfpick sfphaserot sfdots