adding folder with reproduced figures from Yilmaz' (#311)

* added SConstruct with real data example

* added SConstruct with Ray Abma's synthetic data example

* Colab Notebook for seismic trace interpolation using machine learning

* added google colab notebook for seismic trace interpolation using machine learning

* adding folder with reproduced figures from Oz Yilmaz' Seismic Analysis

book/hackathon/2024soundfx/real/SConstruct

@@ -0,0 +1,39 @@
+from rsf.proj import *
+
+# Fetch & plot dataset
+Fetch('elf-stk2.rsf','masha')
+Flow('elf','elf-stk2','dd form=native | put unit1=s')
+Result('elf','grey')
+
+# Decimate
+Flow('decimated','elf','window j2=4 | put d2=53.3333')
+Result('decimated','grey')
+
+#Result('original_n_decimated',['elf','decimated'],'SideBySideAniso')
+
+# Zoom-in
+Flow('original_w','elf','window max1=1.5 n2=600 | costaper nw1=40')
+Flow('decimated_w','decimated','window max1=1.5 n2=150 | costaper nw1=40')
+Result('original_w','grey title="Original data"')
+Result('decimated_w','grey title="Decimated data"')
+
+#Result('original_n_decimated_w',['elf_w','decimated_w'],'SideBySideAniso')
+
+# T-X interpolation
+Flow(['pad_w','mask_w'],'decimated_w','lpad jump=4 mask=${TARGETS[1]}')
+Flow(['pef_w','lag_w'],'decimated_w','lpef lag=${TARGETS[1]} a=17,5 jump=4')
+Flow('interpolated_w',['pad_w','mask_w','pef_w'],
+     'miss padin=4 filt=${SOURCES[2]} mask=${SOURCES[1]} prec=n')
+Result('interpolated_w','grey')
+
+Flow('pef_aw','pad_w','apef a=15,4 jump=4 rect1=20 rect2=5 niter=200 verb=y')
+Flow('interpolated_aw',['pad_w','pef_aw','mask_w'],'miss4 filt=${SOURCES[1]} mask=${SOURCES[2]} verb=y')
+Result('interpolated_aw','grey')
+
+## Interpolation error
+
+# F-X interpolation
+
+## Interpolation error
+
+End()

book/hackathon/2024soundfx/synthabma/SConstruct

@@ -0,0 +1,85 @@
+from rsf.proj import *
+from rsf.recipes.beg import server as private
+
+###########################################################################
+
+for dat in ('curve','linear','random','marm'):
+    input = 'input.%s.segy' % dat
+    Fetch(input,'ray',private)
+    Flow([dat,'./A'+dat,'./B'+dat],input,
+         '''
+         segyread tape=$SOURCE read=d hfile=${TARGETS[1]} bfile=${TARGETS[2]}
+         ''',stdin=0)
+
+Flow('linear2','linear','window min1=0.5 max1=2.7 | bandpass fhi=60')
+Plot('linear','linear2',
+     '''
+     grey yreverse=y transp=y poly=y label2=Position 
+     title=Input
+     ''')
+Plot('jlinear','linear2',
+     '''
+     window n2=11 f2=23 n1=150 min1=1.35 |
+     put d1=1 o1=675 label1=Sample unit1= |
+     wiggle yreverse=y transp=y poly=y label2=Position wherexlabel=t
+     title=Input wheretitle=b clip=0.0451806 labelsz=5. titlesz=7
+     labelfat=2 font=2 titlefat=2 screenratio=1.2
+     ''')
+
+# nonstationary PWD
+Flow('lindip','linear2','twodip2 order=3 nj1=4 nj2=4 eps=10 gauss=n')
+Flow('lindip2','lindip',
+     'transp | spline n1=240 o1=0 d1=0.25 | transp')
+
+Flow('lin4 linones4','linear2','lpad jump=4 mask=${TARGETS[1]}')
+Flow('lindeal','lin4 lindip2 linones4',
+     'planemis2 dip=${SOURCES[1]} mask=${SOURCES[2]} order=3 verb=y')
+Plot('lindeal','grey yreverse=y transp=y poly=y title=Interpolated')
+
+Result('linear-deal','linear lindeal','SideBySideAniso')
+
+# Stationary T-X PEFs
+Flow('lpef llag','linear2','lpef lag=${TARGETS[1]} a=10,4 jump=4')
+Flow('lscov','lpad lmask lpef',
+     'miss padin=4 filt=${SOURCES[2]} mask=${SOURCES[1]} prec=n')
+Plot('lscov',
+     'grey yreverse=y transp=y poly=y title="Stationary PEF"')
+
+Result('linear-scomp','linear lscov','SideBySideAniso')
+
+# Nonstationary T-X PEFs
+Flow('lpad lmask','linear2','lpad jump=2 mask=${TARGETS[1]}')
+Flow('lapef','lpad','apef a=15,4 jump=2 rect1=20 rect2=5 niter=200 verb=y')
+Flow('lacov','lpad lapef lmask',
+     'miss4 filt=${SOURCES[1]} mask=${SOURCES[2]} verb=y')
+Plot('lacov',
+     '''
+     grey yreverse=y transp=y poly=y label2=Position
+     title="Adaptive PEF"
+     ''')
+
+Plot('jlacov','lacov',
+     '''
+     window n2=22 f2=46 n1=150 min1=1.35 |
+     put d1=1 o1=675 label1=Sample unit1= |
+     wiggle yreverse=y transp=y poly=y label2=Position wherexlabel=t
+     title="Adaptive PEF" wheretitle=b clip=0.0225903 labelsz=5. titlesz=7
+     labelfat=2 font=2 titlefat=2 screenratio=1.2
+     ''')
+
+Result('linear-comp','linear lacov','SideBySideAniso')
+
+# Stationary F-X PEFs
+Flow('lfx','lpad',
+     '''
+     spitz norm=n verb=y
+     ''')
+Plot('lfx',
+     '''
+     grey yreverse=y transp=y poly=y label2=Position
+     title="Adaptive PEF"
+     ''')
+
+Result('linear-fxcomp','linear lfx','SideBySideAniso')
+
+End()

book/yilmaz/README

@@ -0,0 +1 @@
+This document aims to reproduce the synthetic data examples presented on Öz Yilmaz' Seismic Data Analysis, Volume I.

book/yilmaz/SConstruct

@@ -0,0 +1,16 @@
+from rsf.book import *
+from rsf.tex import Paper
+
+chapters=Split(
+'''
+fourier1d
+''')
+
+Book(chapters,
+     author='TCCS Apocrypha',
+     title="Reproducing Oz Yilmaz' Seismic Data Analysis")
+
+Paper('intro')
+
+End(options='book',
+   use='graphicx,color,amsmath,amssymb,amsbsy,hyperref,listings,framed,tabularx')

book/yilmaz/fourier1d/SConstruct

@@ -0,0 +1,3 @@
+from rsf.tex import *
+
+End(color='all')

book/yilmaz/fourier1d/aliasing/SConstruct

@@ -0,0 +1,107 @@
+from rsf.proj import *
+
+"""
+This SCons script aims to reproduce Fig 1.1-7 to Fig 1.1-10, presented in Öz Yilmaz' Seismic Data Analysis, Vol 1.
+
+The first three frames show sinusoidal functions with different frequencies:
+    * 1: Frequency -  25  Hz
+    * 2: Frequency -  75  Hz
+    * 3: Frequency - 150  Hz
+
+The fourth frame shows the composition of two sinusoids with frequencies:
+    * 1: Frequency - 12.5 Hz
+    * 2: Frequency - 75   Hz
+
+All the functions above are sampled at three different rates:
+    * 1: Sampling rate - 2 ms
+    * 2: Sampling rate - 4 ms
+    * 3: Sampling rate - 8 ms
+
+The original figures exemplified aliasing: the consequence of sampling below the Nyquist rate.
+
+*For some reason, the original figures appear to be "smoother" than they should.
+"""
+
+# Sinusoid parameters
+frequency = [25,75,150] 			# frequencies in Hz
+samplingrate = [0.002,0.004,0.008] 		# sampling rates in seconds
+
+for i in range(3):
+    for j in range(3):
+        jrate = samplingrate[j]
+        # Produce sinusoid
+        Flow('sinusoid%g-%gms'%(i+1,jrate*1000),None,
+             '''
+             math n1=%g d1=%f o1=0
+             output="sin(x1*%f*2*(2*asin(1)))" |
+             put label1="Time" unit1="s"
+             '''%(1/jrate+1,jrate,frequency[i])) # pi is defined as 2*asin(1)
+        Plot('sinusoid%g-%gms'%(i+1,jrate*1000),
+             'graph title="Sinusoid %g, sampling rate = %g ms" screenratio=0.75 min1=0 max1=1 min2=-1 max2=1'%(i+1,jrate*1000))
+        Result('sinusoid%g-%gms'%(i+1,jrate*1000),
+               'graph title="Sinusoid %g, sampling rate = %g ms" screenratio=0.5 min1=0 max1=1 min2=-1 max2=1'%(i+1,jrate*1000))
+
+        # Recovering Amplitude spectra through FFT
+        ## Fourier transform
+        fftwindow = int(0.4/jrate) # we window N full cycles
+        Flow('fourier%g-%gms'%(i+1,jrate*1000),'sinusoid%g-%gms'%(i+1,jrate*1000),'window n1=%g | fft1'%fftwindow)
+
+        ## Amplitude spectra
+        Flow('fampspectra%g-%gms'%(i+1,jrate*1000),'fourier%g-%gms'%(i+1,jrate*1000),
+            '''
+            math output="2*abs(input)/%g" | real | put label1="Frequency" unit1="Hz"
+            '''%fftwindow) # normalized by 2/N
+        Plot('fampspectra%g-%gms'%(i+1,jrate*1000),
+             'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=250 min2=0 max2=1')
+        Result('fampspectra%g-%gms'%(i+1,jrate*1000),
+               'graph title="Amplitude Spectrum %g" screenratio=1 min1=0 max1=250 min2=0 max2=1'%(i+1))
+
+        ## Plot side by side
+        Plot('subplot%g-%gms'%(i+1,jrate*1000),'sinusoid%g-%gms fampspectra%g-%gms'%(i+1,jrate*1000,i+1,jrate*1000),
+             'SideBySideAniso')
+        Result('subplot%g-%gms'%(i+1,jrate*1000),'sinusoid%g-%gms fampspectra%g-%gms'%(i+1,jrate*1000,i+1,jrate*1000),
+               'SideBySideAniso')
+
+    Plot('frame%g'%(i+1),['subplot%g-%gms'%(i+1,samplingrate[j]*1000) for j in range(3)],'OverUnderAniso')
+    Result('frame%g'%(i+1),['subplot%g-%gms'%(i+1,samplingrate[j]*1000) for j in range(3)],'OverUnderAniso')
+
+# Composite signal
+for j in range(3):
+    jrate = samplingrate[j]
+    # Produce sinusoid
+    Flow('sinusoid4-%gms'%(jrate*1000),None,
+         '''
+         math n1=%g d1=%f o1=0
+         output="sin(x1*12.5*2*(2*asin(1)))+sin(x1*75*2*(2*asin(1)))" |
+         put label1="Time" unit1="s"
+         '''%(1/jrate+1,jrate)) # pi is defined as 2*asin(1)
+    Plot('sinusoid4-%gms'%(jrate*1000),
+         'graph title="Sinusoid %g, sampling rate = %g ms" screenratio=0.75 min1=0 max1=1 min2=-1 max2=1'%(i+1,jrate*1000))
+    Result('sinusoid4-%gms'%(jrate*1000),
+           'graph title="Sinusoid %g, sampling rate = %g ms" screenratio=0.5 min1=0 max1=1 min2=-1 max2=1'%(i+1,jrate*1000))
+
+    # Recovering Amplitude spectrum through FFT
+    ## Fourier transform
+    fftwindow = int(0.4/jrate) # we window N full cycles
+    Flow('fourier4-%gms'%(jrate*1000),'sinusoid4-%gms'%(jrate*1000),'window n1=%g | fft1'%fftwindow)
+
+    ## Amplitude spectrum
+    Flow('fampspectra4-%gms'%(jrate*1000),'fourier4-%gms'%(jrate*1000),
+         '''
+         math output="2*abs(input)/%g" | real | put label1="Frequency" unit1="Hz"
+         '''%fftwindow) # normalized by 2/N
+    Plot('fampspectra4-%gms'%(jrate*1000),
+         'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=250 min2=0 max2=1')
+    Result('fampspectra4-%gms'%(jrate*1000),
+           'graph title="Amplitude Spectrum %g" screenratio=1 min1=0 max1=250 min2=0 max2=1'%(i+1))
+
+    ## Plot side by side
+    Plot('subplot4-%gms'%(jrate*1000),'sinusoid4-%gms fampspectra4-%gms'%(jrate*1000,jrate*1000),
+         'SideBySideAniso')
+    Result('subplot4-%gms'%(jrate*1000),'sinusoid4-%gms fampspectra4-%gms'%(jrate*1000,jrate*1000),
+           'SideBySideAniso')
+
+Plot('frame4',['subplot4-%gms'%(samplingrate[j]*1000) for j in range(3)],'OverUnderAniso')
+Result('frame4',['subplot4-%gms'%(samplingrate[j]*1000) for j in range(3)],'OverUnderAniso')
+
+End()