From b578b0a09d13ade6d10dddcf44423e64b4022281 Mon Sep 17 00:00:00 2001 From: hcorzopola <55563675+hcorzopola@users.noreply.github.com> Date: Wed, 9 Oct 2024 12:47:42 -0500 Subject: [PATCH] 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 | 39 +++ .../2024soundfx/synthabma/SConstruct | 85 +++++ book/yilmaz/README | 1 + book/yilmaz/SConstruct | 16 + book/yilmaz/fourier1d/SConstruct | 3 + book/yilmaz/fourier1d/aliasing/SConstruct | 107 ++++++ .../yilmaz/fourier1d/freqfiltering/SConstruct | 183 ++++++++++ book/yilmaz/fourier1d/paper.tex | 84 +++++ book/yilmaz/fourier1d/phase/SConstruct | 324 ++++++++++++++++++ book/yilmaz/fourier1d/sinusoids/SConstruct | 77 +++++ book/yilmaz/intro.tex | 5 + 11 files changed, 924 insertions(+) create mode 100644 book/hackathon/2024soundfx/real/SConstruct create mode 100644 book/hackathon/2024soundfx/synthabma/SConstruct create mode 100644 book/yilmaz/README create mode 100644 book/yilmaz/SConstruct create mode 100644 book/yilmaz/fourier1d/SConstruct create mode 100644 book/yilmaz/fourier1d/aliasing/SConstruct create mode 100644 book/yilmaz/fourier1d/freqfiltering/SConstruct create mode 100644 book/yilmaz/fourier1d/paper.tex create mode 100644 book/yilmaz/fourier1d/phase/SConstruct create mode 100644 book/yilmaz/fourier1d/sinusoids/SConstruct create mode 100644 book/yilmaz/intro.tex diff --git a/book/hackathon/2024soundfx/real/SConstruct b/book/hackathon/2024soundfx/real/SConstruct new file mode 100644 index 000000000..604192f7b --- /dev/null +++ b/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() diff --git a/book/hackathon/2024soundfx/synthabma/SConstruct b/book/hackathon/2024soundfx/synthabma/SConstruct new file mode 100644 index 000000000..a1733e324 --- /dev/null +++ b/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() diff --git a/book/yilmaz/README b/book/yilmaz/README new file mode 100644 index 000000000..3b857e7c0 --- /dev/null +++ b/book/yilmaz/README @@ -0,0 +1 @@ +This document aims to reproduce the synthetic data examples presented on Öz Yilmaz' Seismic Data Analysis, Volume I. diff --git a/book/yilmaz/SConstruct b/book/yilmaz/SConstruct new file mode 100644 index 000000000..60220efe0 --- /dev/null +++ b/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') diff --git a/book/yilmaz/fourier1d/SConstruct b/book/yilmaz/fourier1d/SConstruct new file mode 100644 index 000000000..2f8ee2166 --- /dev/null +++ b/book/yilmaz/fourier1d/SConstruct @@ -0,0 +1,3 @@ +from rsf.tex import * + +End(color='all') diff --git a/book/yilmaz/fourier1d/aliasing/SConstruct b/book/yilmaz/fourier1d/aliasing/SConstruct new file mode 100644 index 000000000..43078ed53 --- /dev/null +++ b/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() diff --git a/book/yilmaz/fourier1d/freqfiltering/SConstruct b/book/yilmaz/fourier1d/freqfiltering/SConstruct new file mode 100644 index 000000000..552c82288 --- /dev/null +++ b/book/yilmaz/fourier1d/freqfiltering/SConstruct @@ -0,0 +1,183 @@ +from rsf.proj import * + +""" +This SCons script aims to reproduce Fig 1.1-21 to Fig 1.1-?, presented in Öz Yilmaz' Seismic Data Analysis, Vol 1. +""" + +# Get wavelet from summation of sinusoids +def wavelet(n): + wavelet = 'wavelet-%gHz'%(n) + Flow(wavelet,sinusoids[:n], + ''' + add ${SOURCES[1:%g]} | math output=input/%g + '''%(n,n)) + Plot(wavelet,'graph title="" screenratio=1 min2=-1 max2=1') + Result(wavelet,'graph title="" screenratio=2 min2=-1 max2=1') + + waveletFrame = wavelet+'-frame' + Flow(waveletFrame,[*sinusoids[:n],wavelet], + ''' + cat axis=2 d=1 o=1 ${SOURCES[1:%g]} | math output="input+x2" + '''%(n+1)) + Plot(waveletFrame,'graph') + Result(waveletFrame,'graph') + +# Get boxcar amplitude spectra +def get_boxcar(freq1,freq2): + wavelet = 'wavelet-boxcar-%gHz-%gHz'%(freq1,freq2) + fampspectra = 'fampspectra-boxcar-%gHz-%gHz'%(freq1,freq2) + fourier = 'fourier-boxcar-%gHz-%gHz'%(freq1,freq2) + fampspectra2 = 'fampspectra-'+wavelet + fampspectra3 = 'overlay-'+fampspectra2 + frame = 'famp-n-'+wavelet + ## A) Define a boxcar amplitude spectrum... + Flow(fampspectra,None, + ''' + math n1=257 d1=0.976562 o1=0 output="1" | + put label1="Frequency" unit1="Hz" | + cut max1=%g | cut min1=%g + '''%(freq1,freq2)) + Plot(fampspectra, + 'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=120 min2=0 max2=1.2') + Result(fampspectra, + 'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=120 min2=0 max2=1.2') + ## ...and a zero-phase spectrum + Flow(fourier+'-r',fampspectra, + ''' + math output="input" + n1=257 d1=0.976562 o1=0 + ''') + Flow(fourier+'-i',fampspectra, + ''' + math output="0" + n1=257 d1=0.976562 o1=0 + ''') + Flow(fourier,[fourier+'-r',fourier+'-i'], + ''' + cmplx $SOURCES + ''') + ## B) Apply inverse FFT and obtain a filter operator + Flow(wavelet,fourier,'fft1 inv=y | math output="input*257/%g" | put o1=-0.5 | rotate rot1=250'%(freq2-freq1)) # rotate shifts the wavelet to be centered at 0 + Plot(wavelet,'graph title="" screenratio=1 min1=-0.2 max1=0.2 min2=-1 max2=1') + Result(wavelet,'graph title="" screenratio=2 min2=-1 max2=1') + ## C) Truncate the operator, and D) Apply forward FFT and compute th amplitude spectrum of the truncated operator + Flow(fampspectra2,wavelet, + ''' + window min1=-0.2 max1=0.2 | fft1 | + math output="2*abs(input)/%g" | real | put label1="Frequency" unit1="Hz" + '''%(512/(freq2-freq1))) # normalized by 2/N + Plot(fampspectra2, + 'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=120 min2=0 max2=1.2') + Result(fampspectra2, + 'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=120 min2=0 max2=1.2') + ## Overlay amplitude spectra + Plot(fampspectra3,[fampspectra, fampspectra2],'Overlay') + Result(fampspectra3,[fampspectra, fampspectra2],'Overlay') + ## Plot wavelets and amplitude spectra + Plot(frame,[wavelet, fampspectra3],'OverUnderAniso') + Result(frame,[wavelet, fampspectra3],'OverUnderAniso') + +# Apply boxcar filter +def get_trapezoidal(freq1,freq2,freq3,freq4): + # Multiply wavelet amplitude spectra by + pass + +# Get wavelet's Fourier spectra +def get_fspectra(freq): + ## File names + wavelet = 'wavelet-%gHz'%(freq) + fourier = 'fourier-'+wavelet + argf = 'argf-'+wavelet + fmask = 'fmask-'+wavelet + fampspectra = 'fampspectra-'+wavelet + fphispectra = 'fphispectra-'+wavelet + fspectra = 'fspectra-'+wavelet + frame = 'famp-n-'+wavelet + ## Fourier transform + Flow('fourier-'+wavelet,wavelet,'fft1') + + ## Amplitude spectra + Flow(fampspectra,fourier, + ''' + math output="2*abs(input)/%g" | real | put label1="Frequency" unit1="Hz" + '''%(512/freq)) # normalized by 2/N + Plot(fampspectra, + 'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=120 min2=0 max2=1.2') + Result(fampspectra, + 'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=120 min2=0 max2=1.2') + + ## Phase spectra + """ + This one is trickier. On paper, phi=arg(F(w)). However, arg(z) is not well defined + for |z|=0. + Hence, we first compute arg(F), and then we mute the values for which |z|=0. + """ + Flow(argf,fourier, + ''' + math output="arg(input)*90/asin(1)" | real | put label1="Frequency" unit1="Hz" + ''') # multiplied by (90/arcsin(1)) to convert phi from radians to degrees + Flow(fmask,fampspectra, + ''' + mask min=%f | dd type=float + '''%(1/2)) # mask with non-zero values + Flow(fphispectra,[argf, fmask], + ''' + math mask=${SOURCES[1]} output="input*mask" + ''') + Plot(fphispectra, + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=120 min2=-180 max2=180 d2num=90') + Result(fphispectra, + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=120 min2=-180 max2=180 d2num=90') + + ## Plot Fourier spectra + Plot(fspectra,[fampspectra, fphispectra],'SideBySideAniso') + Result(fspectra,[fampspectra, fphispectra],'SideBySideAniso') + + ## Plot wavelets and amplitude spectra + Plot(frame,[wavelet, fampspectra],'OverUnderAniso') + Result(frame,[wavelet, fampspectra],'OverUnderAniso') + +# Compute sinusoids up to Nyquist frequency (i.e. 125 Hz) +nfreq = 125 +dt = 0.002 # sampling rate; should be 0.004, but we will double the sampling rate for the next example +sinusoids = [] +for i in range(nfreq): + isinusoid = 'sinusoid-%ghz'%(i+1) + # Produce sinusoid with frequency of 'i' Hz + Flow(isinusoid,None, + ''' + math n1=%g d1=%f o1=-0.5 + output="cos(x1*%f*2*(2*asin(1)))" | + put label1="Time" unit1="s" + '''%(1/dt+1,dt,i+1)) # pi is defined as 2*asin(1) + sinusoids.append(isinusoid) +# Result(isinusoid,'graph title="" screenratio=2'%(i+1)) # uncomment to visualize individual sinusoids + +# Fig 1.1-21 +wavelet(2) # 'wavelet-2Hz-frame' +wavelet(4) # 'wavelet-4Hz-frame' +wavelet(8) # 'wavelet-8Hz-frame' +wavelet(16) # 'wavelet-16Hz-frame' +wavelet(32) # 'wavelet-32Hz-frame' +Result('frame-1-1-21',['wavelet-2Hz-frame','wavelet-4Hz-frame','wavelet-8Hz-frame','wavelet-16Hz-frame','wavelet-32Hz-frame'],'SideBySideAniso') +# Fig 1.1-22 +wavelet(125) # 'wavelet-125Hz-frame' +Result('frame-1-1-22b',['wavelet-2Hz-frame','wavelet-4Hz-frame','wavelet-8Hz-frame','wavelet-16Hz-frame','wavelet-32Hz-frame','wavelet-125Hz-frame'],'SideBySideAniso') # own interpretation; not in the original + +# Fig 1.1-23 +wavelet(30) # 'wavelet-30Hz-frame' +get_fspectra(30)# 'famp-n-wavelet-30Hz' +wavelet(40) # 'wavelet-40Hz-frame' +get_fspectra(40)# 'famp-n-wavelet-40Hz' +wavelet(50) # 'wavelet-50Hz-frame' +get_fspectra(50)# 'famp-n-wavelet-50Hz' +wavelet(60) # 'wavelet-60Hz-frame' +get_fspectra(60)# 'famp-n-wavelet-60Hz' +wavelet(70) # 'wavelet-70Hz-frame' +get_fspectra(70)# 'famp-n-wavelet-70Hz' +Result('frame-1-1-23',['famp-n-wavelet-30Hz','famp-n-wavelet-40Hz','famp-n-wavelet-50Hz','famp-n-wavelet-60Hz','famp-n-wavelet-70Hz'],'SideBySideAniso') + +# Boxcars +get_boxcar(1,30) + +End() diff --git a/book/yilmaz/fourier1d/paper.tex b/book/yilmaz/fourier1d/paper.tex new file mode 100644 index 000000000..0f13a82c5 --- /dev/null +++ b/book/yilmaz/fourier1d/paper.tex @@ -0,0 +1,84 @@ +\title{Fundamentals of Signal Processing} + +\maketitle + +\begin{abstract} +This chapter reproduces some of the synthetic data examples presented in {\"O}z Yilmaz' Seismic Data Analysis, Chapter 1: Fundamentals of Signal Processing. +\end{abstract} + +\section{The 1-D Fourier Transform} +\inputdir{sinusoids} + +We reproduced the example used by Yilmaz to illustrate how a sinusoidal time function can be described just in terms of its frequency, peak amplitude, and phase. The original figures showed the recorded motion of hypothetical springs when a weight is attached to them. Figure~\ref{fig:frame1} and Figure~\ref{fig:frame3} show the motion for the same spring-weight setup, recorded with a delay of 20 ms for the latter. Figure~\ref{fig:frame2} shows a stiffer spring, which results in a shorter period and a smaller peak amplitude. + +The peak amplitudes, periods and lags of the three sinusoidal functions are: +\begin{itemize} + \item Figure~\ref{fig:frame1}: Amplitude - 0.8 u, Period - 0.080 s, Lag - 0 ms + \item Figure~\ref{fig:frame2}: Amplitude - 0.4 u, Period - 0.040 s, Lag - 0 ms + \item Figure~\ref{fig:frame3}: Amplitude - 0.8 u, Period - 0.080 s, Lag - 20 ms +\end{itemize} + +\multiplot{3}{frame1,frame2,frame3}{width=0.3\textwidth}{Sinusoids from FIG 1.1-1, and their amplitude and phase spectra.} + +\pagebreak + +\subsection{Frequency Aliasing} +\inputdir{aliasing} + +We reproduced the figures that Yilmaz used to exemplify aliasing: the consequence of sampling below the Nyquist rate. + +Figures~\ref{fig:frame1},~\ref{fig:frame2}, and ~\ref{fig:frame3} show sinusoidal functions with different frequencies: +\begin{itemize} + \item Figure~\ref{fig:frame1}: 25 Hz + \item Figure~\ref{fig:frame2}: 75 Hz + \item Figure~\ref{fig:frame3}: 150 Hz +\end{itemize} + +Figure~\ref{fig:frame4} shows the composition of two sinusoids with frequencies: +\begin{itemize} + \item 12.5 Hz + \item 75 Hz +\end{itemize} + +All the functions above are sampled at three different rates: +\begin{itemize} + \item 2 ms + \item 4 ms + \item 8 ms +\end{itemize} + +\plot{frame1}{width=\textwidth}{25-Hz sinusoid sampled at 2, 4, and 8 ms. As shown on FIG 1.1-7.} +\plot{frame2}{width=\textwidth}{75-Hz sinusoid sampled at 2, 4, and 8 ms. As shown on FIG 1.1-8.} +\plot{frame3}{width=\textwidth}{150-Hz sinusoid sampled at 2, 4, and 8 ms. As shown on FIG 1.1-9} +\plot{frame4}{width=\textwidth}{Sum of a 12.5-Hz and 75-Hz sinusoids sampled at 2, 4, and 8 ms. As shown on FIG 1.1-10} + +\pagebreak + +\subsection{Phase Considerations} +\inputdir{phase} + +The figures of this subsection show the effects of phase shifts on a time signal. + +Figures ~\ref{fig:wavelet-0ms-frame}, ~\ref{fig:wavelet-200ms-frame}, ~\ref{fig:wavelet-90deg-shift-frame}, and ~\ref{fig:wavelet-200ms-90deg-shift-frame} show the summation of sinusoids with frequencies ranging from 1 to 32 Hz to form a wavelet. + +Figures ~\ref{fig:figure-1-1-13}, ~\ref{fig:figure-1-1-15}, and ~\ref{fig:figure-1-1-18} show the effects of applying constant and linear phase shifts on zero-phase wavelets. + +\plot{wavelet-0ms-frame}{width=\textwidth}{Summation of a discrete number of sinusoids (1-32 Hz) with no phase-lag. As shown on FIG 1.1-11.} +\plot{wavelet-200ms-frame}{width=\textwidth}{Summation of a discrete number of sinusoids (1-32 Hz) with a constant -0.2 s time-lag. As shown on FIG 1.1-12.} +\plot{wavelet-90deg-shift-frame}{width=\textwidth}{Summation of a discrete number of sinusoids (1-32 Hz) with a constant 90 degree phase-shift. As shown on FIG 1.1-14.} +\plot{wavelet-200ms-90deg-shift-frame}{width=\textwidth}{Summation of a discrete number of sinusoids (1-32 Hz) with a constant -0.2 s time-lag and a constant 90 degree phase-shift. As shown on FIG 1.1-17.} + +\plot{figure-1-1-13}{width=\textwidth}{Linear phase shifts applied to shift a wavelet in time. As shown on FIG 1.1-13.} +\plot{figure-1-1-15}{width=\textwidth}{Constant phase shifts applied to a wavelet change its shape. As shown on FIG 1.1-15.} +\plot{figure-1-1-18}{width=\textwidth}{Effects of applyign both constant and linear phase shifts to a wavelet. As shown on FIG 1.1-18.} + +\pagebreak + +\subsection{Frequency Filtering} +\inputdir{freqfiltering} + +In Yilmaz' book, this subsection discusses some fundamental ideas about filtering in the frequency domain. However, the synthetic examples presented only illustrate how to produce synthetic wavelets from a summation of a finite set of zero-phase sinusoids, as shown in Figure~\ref{fig:frame-1-1-21}. The key idea that, the wavelets become more compact in time as we increase their frequency bandwith (i.e. the number of sinusoids in the summation). + +\plot{frame-1-1-21}{width=\textwidth}{Synthetic wavelet resulting from the summation of different number of zero-phase sinusoids with identical peak amplitudes. In all subfigures, the upper most trace is the resulting wavelet. As shown on FIG 1.1-21.} + +Fig 1.1-21 to 1.1-27 diff --git a/book/yilmaz/fourier1d/phase/SConstruct b/book/yilmaz/fourier1d/phase/SConstruct new file mode 100644 index 000000000..5e75de1dc --- /dev/null +++ b/book/yilmaz/fourier1d/phase/SConstruct @@ -0,0 +1,324 @@ +from rsf.proj import * + +""" +This SCons script aims to reproduce Fig 1.1-11 to Fig 1.1-18, presented in Öz Yilmaz' Seismic Data Analysis, Vol 1. + +Figures 1.1-11,12,14, and 17 show the summation of sinusoids with frequencies ranging from 1 to 32 Hz to form a wavelet. + +Figures 1.1-13,15, and 18 show the effects of applying constant and linear phase shifts on zero-phase wavelets. +""" + +# FUNCTIONS TO COMPUTE WAVELETS AS SINUSOID SUMATIONS ################################### +""" +The formulation of both linear and constant phase shifts is more elegant in the Fourier +domain. However, I am including these functions to illustrate what is happening in the +time domain. Reproducing Figs 1.1-11,12,14, and 17 with these functions also reduces +the number of Fourier transforms used. +""" +# Wavelets from constant lag sinusoid sumations +def constant_lag_wavelet(lag): + # Sinusoid parameters + dt = 0.004 # sampling rate + sinusoids = [] + for i in range(32): + isinusoid = 'sinusoid-%ghz-%gms'%(i+1,lag*1000) + # Produce sinusoid with frequency of 'i' Hz + Flow(isinusoid,None, + ''' + math n1=%g d1=%f o1=-0.5 + output="cos((x1+%f)*%f*2*(2*asin(1)))" | + put label1="Time" unit1="s" + '''%(1/dt+1,dt,lag,i+1)) # pi is defined as 2*asin(1) + sinusoids.append(isinusoid) +# Result(isinusoid,'graph title="" screenratio=2'%(i+1)) # uncomment to visualize individual sinusoids + + wavelet = 'wavelet-%gms'%(lag*1000) + Flow(wavelet,sinusoids, + ''' + add ${SOURCES[1:32]} | math output=input/32 + ''') + Plot(wavelet,'graph title="" screenratio=1 min2=-1 max2=1') + Result(wavelet,'graph title="" screenratio=2 min2=-1 max2=1') + + waveletFrame = wavelet+'-frame' + Flow(waveletFrame,[*sinusoids,wavelet], + ''' + cat axis=2 d=1 o=1 ${SOURCES[1:33]} | math output="input+x2" + ''') + Plot(waveletFrame,'graph title=""') + Result(waveletFrame,'graph title=""') + +# Wavelets from constant phase shifted sinusoid sumations +def constant_phase_shift_wavelet(shift): + # Sinusoid parameters + dt = 0.004 # sampling rate + ratio = 360/shift # get ratio between full cycle and shift + lags = [1/((i+1)*ratio) for i in range(32)] # get lag in seconds <- (1/(cycles per second*fraction of full cycle)) + sinusoids = [] + for i in range(32): + isinusoid = 'sinusoid-%ghz-%gdeg-shift'%(i+1,shift) + # Produce sinusoid with frequency of 'i' Hz + Flow(isinusoid,None, + ''' + math n1=%g d1=%f o1=-0.5 + output="cos((x1+%f)*%f*2*(2*asin(1)))" | + put label1="Time" unit1="s" + '''%(1/dt+1,dt,lags[i],i+1)) # pi is defined as 2*asin(1) + sinusoids.append(isinusoid) +# Result(isinusoid,'graph title="" screenratio=2'%(i+1)) # uncomment to visualize individual sinusoids + + wavelet = 'wavelet-%gdeg-shift'%(shift) + Flow(wavelet,sinusoids, + ''' + add ${SOURCES[1:32]} | math output=input/32 + ''') + Plot(wavelet,'graph title="" screenratio=1 min2=-1 max2=1') + Result(wavelet,'graph title="" screenratio=2 min2=-1 max2=1') + + waveletFrame = wavelet+'-frame' + Flow(waveletFrame,[*sinusoids,wavelet], + ''' + cat axis=2 d=1 o=1 ${SOURCES[1:33]} | math output="input+x2" + ''') + Plot(waveletFrame,'graph title=""') + Result(waveletFrame,'graph title=""') + +# Wavelets from constant phase shifted sinusoid sumations +def mixed_phase_shift_wavelet(shift,lag): + # Sinusoid parameters + dt = 0.004 # sampling rate + ratio = 360/shift # get ratio between full cycle and shift + lags = [1/((i+1)*ratio)+lag for i in range(32)] # get lag in seconds <- (1/(cycles per second*fraction of full cycle)) + sinusoids = [] + for i in range(32): + isinusoid = 'sinusoid-%ghz-%gms-%gdeg-shift'%(i+1,lag*1000,shift) + # Produce sinusoid with frequency of 'i' Hz + Flow(isinusoid,None, + ''' + math n1=%g d1=%f o1=-0.5 + output="cos((x1+%f)*%f*2*(2*asin(1)))" | + put label1="Time" unit1="s" + '''%(1/dt+1,dt,lags[i],i+1)) # pi is defined as 2*asin(1) + sinusoids.append(isinusoid) +# Result(isinusoid,'graph title="" screenratio=2'%(i+1)) # uncomment to visualize individual sinusoids + + wavelet = 'wavelet-%gms-%gdeg-shift'%(lag*1000,shift) + Flow(wavelet,sinusoids, + ''' + add ${SOURCES[1:32]} | math output=input/32 + ''') + Plot(wavelet,'graph title="" screenratio=1 min2=-1 max2=1') + Result(wavelet,'graph title="" screenratio=2 min2=-1 max2=1') + + waveletFrame = wavelet+'-frame' + Flow(waveletFrame,[*sinusoids,wavelet], + ''' + cat axis=2 d=1 o=1 ${SOURCES[1:33]} | math output="input+x2" + ''') + Plot(waveletFrame,'graph title=""') + Result(waveletFrame,'graph title=""') +######################################################################################### + +# FUNCTIONS TO APPLY PHASE SHIFTS ON THE FOURIER DOMAIN ################################# +""" +The functions below apply phase shifts on the Fourier domain. +Results can be improved. +""" +def get_fspectra(): + ## Fourier transform + Flow('fourier','wavelet-0ms','fft1') + + ## Amplitude spectra + Flow('fampspectra','fourier', + ''' + math output="2*abs(input)/%g" | real | put label1="Frequency" unit1="Hz" + '''%100) # normalized by 2/N + Plot('fampspectra', + 'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=125 min2=0 max2=1') + Result('fampspectra', + 'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=125 min2=0 max2=1') + + ## Phase spectra + """ + This one is trickier. On paper, phi=arg(F(w)). However, arg(z) is not well defined + for |z|=0. + Hence, we first compute arg(F), and then we mute the values for which |z|=0. + """ + Flow('argf','fourier', + ''' + math output="arg(input)*90/asin(1)" | real | put label1="Frequency" unit1="Hz" + ''') # multiplied by (90/arcsin(1)) to convert phi from radians to degrees + Flow('fmask','fampspectra', + ''' + mask min=%f | dd type=float + '''%(1/2)) # mask with non-zero values + Flow('fphispectra','argf fmask', + ''' + math mask=${SOURCES[1]} output="input*mask" + ''') + Plot('fphispectra', + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-180 max2=180 d2num=90') + Result('fphispectra', + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-180 max2=180 d2num=90') + + Plot('fspectra','fampspectra fphispectra','SideBySideAniso') + Result('fspectra','fampspectra fphispectra','SideBySideAniso') + +def get_fphispectraFrames(): + # 2 pi frame + Plot('fphispectra-2pi','fphispectra', + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-360 max2=360 d2num=90') + Plot('frame-2pi',['fphispectra-2pi','wavelet-0ms'],'OverUnderAniso') + Result('frame-2pi',['fphispectra-2pi','wavelet-0ms'],'OverUnderAniso') + # 10 pi frame + Plot('fphispectra-10pi','fphispectra', + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-1800 max2=1800 d2num=90') + Plot('frame-10pi',['fphispectra-10pi','wavelet-0ms'],'OverUnderAniso') + Result('frame-10pi',['fphispectra-10pi','wavelet-0ms'],'OverUnderAniso') + # 20 pi frame + Plot('fphispectra-20pi','fphispectra', + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-3600 max2=3600 d2num=90') + Plot('frame-20pi',['fphispectra-20pi','wavelet-0ms'],'OverUnderAniso') + Result('frame-20pi',['fphispectra-20pi','wavelet-0ms'],'OverUnderAniso') + +# Linear phase shift +def linear_phase_shift(slope): + m = 125/360 + Flow('fphispectra-%gslope-shift'%(slope*m),'fphispectra', + ''' + math output="input+(%f*x1)" + '''%slope) + Plot('fphispectra-%gslope-shift'%(slope*m), + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-3600 max2=3600 d2num=1800') + Result('fphispectra-%gslope-shift'%(slope*m), + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-3600 max2=3600 d2num=1800') + + Flow('fourier-%gslope-shift-r'%(slope*m),'fampspectra fphispectra-%gslope-shift'%(slope*m), + ''' + math a=${SOURCES[0]} phi=${SOURCES[1]} output="a*cos(phi*asin(1)/90)" + n1=129 d1=0.976562 o1=0 + ''') + + Flow('fourier-%gslope-shift-i'%(slope*m),'fampspectra fphispectra-%gslope-shift'%(slope*m), + ''' + math a=${SOURCES[0]} phi=${SOURCES[1]} output="a*sin(phi*asin(1)/90)" + n1=129 d1=0.976562 o1=0 + ''') + + Flow('fourier-%gslope-shift'%(slope*m),'fourier-%gslope-shift-r fourier-%gslope-shift-i'%(slope*m,slope*m), + ''' + cmplx $SOURCES + ''') + + Flow('wavelet-%gslope-shift-fshifted'%(slope*m),'fourier-%gslope-shift'%(slope*m),'fft1 inv=y | math output="input*129"') + Plot('wavelet-%gslope-shift-fshifted'%(slope*m),'graph title="" screenratio=1 min2=-1 max2=1') + Result('wavelet-%gslope-shift-fshifted'%(slope*m),'graph title="" screenratio=2 min2=-1 max2=1') + + Plot('frame-%gslope-shift'%(slope*m),['fphispectra-%gslope-shift'%(slope*m),'wavelet-%gslope-shift-fshifted'%(slope*m)],'OverUnderAniso') + Result('frame-%gslope-shift'%(slope*m),['fphispectra-%gslope-shift'%(slope*m),'wavelet-%gslope-shift-fshifted'%(slope*m)],'OverUnderAniso') + +# Constant phase shift +def constant_phase_shift(shift): + Flow('fphispectra-%gdeg-shift'%shift,'fphispectra', + ''' + math output="input+(%f)" + '''%shift) + Plot('fphispectra-%gdeg-shift'%shift, + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-360 max2=360 d2num=90') + Result('fphispectra-%gdeg-shift'%shift, + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-360 max2=360 d2num=90') + + Flow('fourier-%gdeg-shift-r'%shift,'fampspectra fphispectra-%gdeg-shift'%shift, + ''' + math a=${SOURCES[0]} phi=${SOURCES[1]} output="a*cos(phi*asin(1)/90)" + n1=129 d1=0.976562 o1=0 + ''') + + Flow('fourier-%gdeg-shift-i'%shift,'fampspectra fphispectra-%gdeg-shift'%shift, + ''' + math a=${SOURCES[0]} phi=${SOURCES[1]} output="a*sin(phi*asin(1)/90)" + n1=129 d1=0.976562 o1=0 + ''') + + Flow('fourier-%gdeg-shift'%shift,'fourier-%gdeg-shift-r fourier-%gdeg-shift-i'%(shift,shift), + ''' + cmplx $SOURCES + ''') + + Flow('wavelet-%gdeg-shift-fshifted'%shift,'fourier-%gdeg-shift'%shift,'fft1 inv=y | math output="input*129"') + Plot('wavelet-%gdeg-shift-fshifted'%shift,'graph title="" screenratio=1 min2=-1 max2=1') + Result('wavelet-%gdeg-shift-fshifted'%shift,'graph title="" screenratio=2 min2=-1 max2=1') + + Plot('frame-%gdeg-shift'%shift,['fphispectra-%gdeg-shift'%shift,'wavelet-%gdeg-shift-fshifted'%shift],'OverUnderAniso') + Result('frame-%gdeg-shift'%shift,['fphispectra-%gdeg-shift'%shift,'wavelet-%gdeg-shift-fshifted'%shift],'OverUnderAniso') + +# Mixed phase shift +def mixed_phase_shift(shift,slope): + m = 125/360 + Flow('fphispectra-%gslope-%gdeg-shift'%(slope*m,shift),'fphispectra', + ''' + math output="input+(%f*x1)+(%f)" + '''%(slope,shift)) + Plot('fphispectra-%gslope-%gdeg-shift'%(slope*m,shift), + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-1800 max2=1800 d2num=1800') + Result('fphispectra-%gslope-%gdeg-shift'%(slope*m,shift), + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-1800 max2=1800 d2num=1800') + + Flow('fourier-%gslope-%gdeg-shift-r'%(slope*m,shift),'fampspectra fphispectra-%gslope-%gdeg-shift'%(slope*m,shift), + ''' + math a=${SOURCES[0]} phi=${SOURCES[1]} output="a*cos(phi*asin(1)/90)" + n1=129 d1=0.976562 o1=0 + ''') + + Flow('fourier-%gslope-%gdeg-shift-i'%(slope*m,shift),'fampspectra fphispectra-%gslope-%gdeg-shift'%(slope*m,shift), + ''' + math a=${SOURCES[0]} phi=${SOURCES[1]} output="a*sin(phi*asin(1)/90)" + n1=129 d1=0.976562 o1=0 + ''') + + Flow('fourier-%gslope-%gdeg-shift'%(slope*m,shift),'fourier-%gslope-%gdeg-shift-r fourier-%gslope-%gdeg-shift-i'%(slope*m,shift,slope*m,shift), + ''' + cmplx $SOURCES + ''') + + Flow('wavelet-%gslope-%gdeg-shift-fshifted'%(slope*m,shift),'fourier-%gslope-%gdeg-shift'%(slope*m,shift),'fft1 inv=y | math output="input*129"') + Plot('wavelet-%gslope-%gdeg-shift-fshifted'%(slope*m,shift),'graph title="" screenratio=1 min2=-1 max2=1') + Result('wavelet-%gslope-%gdeg-shift-fshifted'%(slope*m,shift),'graph title="" screenratio=2 min2=-1 max2=1') + + Plot('frame-%gslope-%gdeg-shift'%(slope*m,shift),['fphispectra-%gslope-%gdeg-shift'%(slope*m,shift),'wavelet-%gslope-%gdeg-shift-fshifted'%(slope*m,shift)],'OverUnderAniso') + Result('frame-%gslope-%gdeg-shift'%(slope*m,shift),['fphispectra-%gslope-%gdeg-shift'%(slope*m,shift),'wavelet-%gslope-%gdeg-shift-fshifted'%(slope*m,shift)],'OverUnderAniso') + +######################################################################################### + +# Wavelets from sinusoid sumations +constant_lag_wavelet(0) # Fig 1.1-11; 'wavelet-0ms-frame' +constant_lag_wavelet(0.2) # Fig 1.1-12; 'wavelet-200ms-frame' +constant_phase_shift_wavelet(90) # Fig 1.1-14; 'wavelet-90deg-shift-frame' +mixed_phase_shift_wavelet(90,0.2) # Fig 1.1-17; 'wavelet-200ms-90deg-shift-frame' + +# Recovering Amplitude & Phase spectra through FFT +get_fspectra() + +# Phase shifts applied on the Fourier domain +## Fig 1.1-13 +linear_phase_shift(1*360/125) # Fig 1.1-13b; 'frame-1slope-shift' +linear_phase_shift(2*360/125) # Fig 1.1-13c; 'frame-2slope-shift' +linear_phase_shift(4*360/125) # Fig 1.1-13d; 'frame-4slope-shift' +linear_phase_shift(8*360/125) # Fig 1.1-13e; 'frame-8slope-shift' +## Fig 1.1-15 +constant_phase_shift(90) # Fig 1.1-15b; 'frame-90deg-shift' +constant_phase_shift(180) # Fig 1.1-15c; 'frame-180deg-shift' +constant_phase_shift(270) # Fig 1.1-15d; 'frame-270deg-shift' +constant_phase_shift(360) # Fig 1.1-15e; 'frame-360deg-shift' +## Fig 1.1-18 +mixed_phase_shift(90,4*360/125) # Fig 1.1-18b; 'frame-4slope-90deg-shift' +mixed_phase_shift(-90,-4*360/125) # Fig 1.1-18c; 'frame--4slope--90deg-shift' + +# Final figures of phase shifts applied on the Fourier domain +## Produce copies of 'fphispectra' with different "zoom" +get_fphispectraFrames() +Result('figure-1-1-13',['frame-20pi','frame-1slope-shift','frame-2slope-shift','frame-4slope-shift','frame-8slope-shift'],'SideBySideAniso') +Result('figure-1-1-15',['frame-2pi','frame-90deg-shift','frame-180deg-shift','frame-270deg-shift','frame-360deg-shift'],'SideBySideAniso') +Result('figure-1-1-18',['frame-20pi','frame-4slope-90deg-shift','frame--4slope--90deg-shift'],'SideBySideAniso') + +End() diff --git a/book/yilmaz/fourier1d/sinusoids/SConstruct b/book/yilmaz/fourier1d/sinusoids/SConstruct new file mode 100644 index 000000000..44adbd3ce --- /dev/null +++ b/book/yilmaz/fourier1d/sinusoids/SConstruct @@ -0,0 +1,77 @@ +from rsf.proj import * + +""" +This SCons script aims to reproduce the subfigures of Fig 1.1-1, presented in Öz Yilmaz' Seismic Data Analysis, Vol 1. + +We produce three sinusoidal functions with different peak amplitudes, periods and lags: + * 1: Amplitude - 0.8 u, Period - 0.080 s, Lag - 0 ms + * 2: Amplitude - 0.4 u, Period - 0.040 s, Lag - 0 ms + * 3: Amplitude - 0.8 u, Period - 0.080 s, Lag - 20 ms + +The original figure showed the recorded motion of hypothetical springs when a weight is attached to them. Sinusoids 1 and 3 represent the same spring-weight setup, but recorded with a delay of 20 ms. Sinusoid 2 represent a stiffer spring, which results in a shorter period and a smaller peak amplitude. + +This hypothetical experiment was used by Yilmaz to exemplify how a sinusoidal time function can be described just in terms of its frequency, peak amplitude, and phase. +""" + +# Sinusoid parameters +amplitude = [0.8,0.4,0.8] # peak amplitudes of the sinusoids +frequency = [12.5,25,12.5] # frequencies in Hz +lag = [0,0,0.020] # lags in seconds + +for i in range(3): + # Produce sinusoid + Flow('sinusoid%g'%(i+1),None, + ''' + math n1=1001 d1=0.001 o1=0 + output="%f*cos((x1-%f)*%f*2*(2*asin(1)))" | + put label1="Time" unit1="s" + '''%(amplitude[i],lag[i],frequency[i])) # pi is defined as 2*asin(1) + Plot('sinusoid%g'%(i+1), + 'graph title="Sinusoid %g" screenratio=0.75 min1=0 max1=1 min2=-1 max2=1'%(i+1)) + Result('sinusoid%g'%(i+1), + 'graph title="Sinusoid %g" screenratio=0.5 min1=0 max1=1 min2=-1 max2=1'%(i+1)) + + # Recovering Amplitude & Phase spectra through FFT + ## Fourier transform + fftwindow = 800 # we window 10 full cycles + Flow('fourier%g'%(i+1),'sinusoid%g'%(i+1),'window n1=%g | fft1'%fftwindow) + + ## Amplitude spectra + Flow('fampspectra%g'%(i+1),'fourier%g'%(i+1), + ''' + math output="2*abs(input)/%g" | real | put label1="Frequency" unit1="Hz" + '''%fftwindow) # normalized by 2/N + Plot('fampspectra%g'%(i+1), + 'graph title="Amplitude Spectrum" screenratio=1 min1=0 max1=125 min2=0 max2=1') + Result('fampspectra%g'%(i+1), + 'graph title="Amplitude Spectrum %g" screenratio=1 min1=0 max1=125 min2=0 max2=1'%(i+1)) + + ## Phase spectra + """ + This one is trickier. On paper, phi=arg(F(w)). However, arg(z) is not well defined + for |z|=0. + Hence, we first compute arg(F), and then we mute the values for which |z|=0. + """ + Flow('argf%g'%(i+1),'fourier%g'%(i+1), + ''' + math output="arg(input)*90/asin(1)" | real | put label1="Frequency" unit1="Hz" + ''') # multiplied by (90/arcsin(1)) to convert phi from radians to degrees + Flow('fmask%g'%(i+1),'fampspectra%g'%(i+1), + ''' + mask min=%f | dd type=float + '''%(amplitude[i]/2)) # mask with non-zero values + Flow('fphispectra%g'%(i+1),'argf%g fmask%g'%(i+1,i+1), + ''' + math mask=${SOURCES[1]} output="input*mask" + ''') + Plot('fphispectra%g'%(i+1), + 'graph title="Phase Spectrum" screenratio=1 min1=0 max1=125 min2=-180 max2=180 d2num=90') + Result('fphispectra%g'%(i+1), + 'graph title="Phase Spectrum %g" screenratio=1 min1=0 max1=125 min2=-180 max2=180 d2num=90'%(i+1)) + + Plot('fspectra%g'%(i+1),'fampspectra%g fphispectra%g'%(i+1,i+1),'SideBySideAniso') + Result('fspectra%g'%(i+1),'fampspectra%g fphispectra%g'%(i+1,i+1),'SideBySideAniso') + + Result('frame%g'%(i+1),'sinusoid%g fspectra%g'%(i+1,i+1),'TwoRows') + +End() diff --git a/book/yilmaz/intro.tex b/book/yilmaz/intro.tex new file mode 100644 index 000000000..33429cae7 --- /dev/null +++ b/book/yilmaz/intro.tex @@ -0,0 +1,5 @@ +\title{Introduction} + +This is a placeholder file for an \emph{Introduction} to the final document. + +The aim of the final version of this document is to reproduce the synthetic figures presented by {\"O}z Yilmaz in his book Seismic Data Analysis, Chapter 1 Fundamentals of Signal Processing.