Homework 3

Next: Interpolation after coordinate transformation Up: Homework 3 Previous: Theoretical part

# Fourier compression

In this exercise, we will use a depth slice selected from a 3-D seismic volume and shown in Figure 1 (Hall, 2007). Notice a channel structure.

data
Figure 1.
Seismic depth slice with a channel structure.

fft
Figure 2.
Absolute value of the Fourier transform of the seismic slice from Figure 1. The circle inside shows a window selected for compression.

The goal of your assignment is to find a compressed representation of the data in the Fourier transform domain. Figure 2 shows the Fourier transform of the data from Figure 1. We can see that most of the energy gets concentrated near the center (zero frequency).

There are two alternative ways to compress data in the Fourier domain:

• One approach is to select a range of frequencies that contain the most important information. An advantage of this approach is the ability to subsample the original data by transforming back from a windowed range of frequencies. The results from this method are shown in Figure 3.
• Another approach is to zero all Fourier coefficients below a certain threshold value, regardless of which frequencies they represent. The results from this method are shown in Figure 4. Figure 5 shows the selected threshold plotted against the histogram of Fourier coefficients.

sig,cut
Figure 3.
Data separated into signal (a) and noise (b) by applying Fourier compression with windowing.

thr,noi
Figure 4.
Data separated into signal (a) and noise (b) by applying Fourier compression with thresholding.

hist
Figure 5.
Normalized histogram of Fourier coefficients (by absolute value). The vertical line shows the selected threshold.

1. Change directory to hw3/fourier.
2. Run
```scons view
```
to reproduce the figures on your screen.
3. Modify the SConstruct file to decrease the size of the window so that the noise level increases in Figure 3b. How do you measure the noise level? Find a level that you find negligibly small.
4. Modify the SConstruct file to increase the threshold value so that the compression achieves the same quality as in the previous case. The noise level in Figure 4b should match that in Figure 3b.
5. Compare the number of nonzero Fourier coefficients in both cases. Which method achieves a better compression?
6. EXTRA CREDIT for finding a way for a better compression of the data in the Fourier domain. Your data reconstruction should have the same noise level, yet the number of non-zero coefficients in the Fourier domain should be smaller.

 ```from rsf.proj import * # Get data ########## Fetch('horizon.asc','hall') # Convert format Flow('data','horizon.asc', ''' echo in=\$SOURCE data_format=ascii_float n1=3 n2=57036 | dd form=native | window n1=1 f1=-1 | add add=-65 | put n2=291 o2=35.031 d2=0.01 label2=y unit2=km n1=196 o1=33.139 d1=0.01 label1=x unit1=km | costaper nw1=25 nw2=25 ''') # Display def plot(title): return ''' grey color=j title="%s" transp=y yreverse=n clip=14 ''' % title Result('data',plot('Horizon')) # 2-D Fourier Transform ####################### Flow('fft','data', 'rtoc | fft3 axis=1 pad=1 | fft3 axis=2 pad=1') Plot('fft', ''' math output="abs(input)" | real | grey title="Fourier Transform" allpos=y screenratio=1 ''') # A. Compression by Windowing ############################# cut = 8 # !!! CHANGE ME !!! # Create a frame Flow('frame','fft','real | math output="sqrt(x1*x1+x2*x2)" ') Plot('frame', ''' contour nc=1 c0=%g plotfat=5 plotcol=3 wantaxis=n wanttitle=n screenratio=1 ''' % cut) Result('fft','fft frame','Overlay') # Cut a hole Flow('fcut','frame fft', ''' mask max=%g | dd type=float | rtoc | mul \${SOURCES[1]} ''' % cut) # Inverse FFT Flow('sig','fcut', 'fft3 axis=2 inv=y | fft3 axis=1 inv=y | real') Result('sig',plot('Compression Signal')) Flow('cut','data sig','add scale=1,-1 \${SOURCES[1]}') Result('cut',plot('Compression Noise') + 'color=I') # B. Compression by Thresholding ################################ thr = 80 # !!! CHANGE ME !!! # Plot histogram Plot('hist','fft', ''' math output="abs(input)" | real | histogram o1=0 d1=%g n1=101 | dd type=float | scale axis=1 | graph title="Scaled Histogram" pad1=n label1= unit1= label2= unit2= ''' % (0.03*thr)) Flow('line.asc',None, 'echo 0 0 0 1 n1=4 data_format=ascii_float in=\$TARGET') Plot('line','line.asc', ''' dd type=complex form=native | graph min1=-1 max1=2 plotcol=5 wantaxis=n wanttitle=n dash=1 ''') Result('hist','hist line','Overlay') # Thresholding Flow('fthr','fft','thr thr=%g' % thr) # Inverse FFT Flow('thr','fthr', 'fft3 axis=2 inv=y | fft3 axis=1 inv=y | real') Result('thr',plot('Compression Signal')) # Subtract from Data Flow('noi','data thr','add scale=1,-1 \${SOURCES[1]}') Result('noi',plot('Compression Noise') + 'color=I') End() ```

 Homework 3

Next: Interpolation after coordinate transformation Up: Homework 3 Previous: Theoretical part

2014-10-02