seismic convolution
The convolution process uses the phase and amplitude obtained earlier, so if you don't want to look at the sweep wavelet, the above steps are not necessary.

Convolution of a Wavelet with the Reflection Coefficients
With the wavelet and reflection coefficient strings both in the time domain, we are now able to create a synthetic seismogram. The process is called convolution. In the time domain, the following steps are needed:
   1. reverse the order of the coefficients of the wavelet.
   2. multiply each of the W values of the reversed wavelet with the first W values of the reflection coefficient string and sum the results; this is the first value of the output trace. (W is the number of samples in the wavelet).
   3. shift the wavelet down the coefficient string by one sample and multiply/sum as in step 2; this is the second output sample.
   4. repeat until the wavelet has been shifted to the end of the reflection coefficient string.
   5. to prevent unwanted end effects, the reflection coefficient string must be padded with zeros at both ends before the convolution is started; the length of the pad is equal to the wavelet length.

NUMERICAL EXAMPLE:

Notice how the tail end of the first reflection interferes with the beginning of the second reflection at data point 6.

This technique, while simple to program, is computer intensive. The number of multiplications equals N squared, where N is the length of the padded reflection coefficient string.

A more efficient technique uses the Fast Fourier Transform, which converts the time series wavelet and reflection coefficient strings to their equivalent frequency domain amplitude and phase spectra. Convolution of the two involves the following:
   1. calculate phase and amplitude of wavelet, using FFT algorithm.
   2. calculate phase and amplitude of reflection coefficients, using FFT algorithm.
   3. multiply the two amplitude spectra and add the two phase spectra.
   4. convert these spectra to the time domain using the inverse FFT algorithm.

This method uses N * log2(N) multiplications instead of N squared. For example, a 5000 point transform takes 61,438 multiplications, while the direct approach takes 25,000,000, a speed improvement of over 600 times. Thus a computer which did the FFT in one minute would take 6.7 hours to compute longhand. For a million data points, the long method would take four weeks compared to a one minute FFT.

A sample program used for performing the FFT is shown below. As written, the program uses single precision numbers, and results may be a bit noisy. Real programs use double precision for the important array variables.

1 REM **********************************************************************
3 REM FAST FOURIER TRANSFORM PROGRAM - FORWARD
5 REM ********* Author: E.R. Crain, P.Eng. 29 Oct 1987 *******************
23 REM **** CAUTION: GWBASIC sets true = NEGATIVE!!!!!!!!!!!!!!!!
24 REM only on arithmetic operators - line 4070 needs fixes if true = POSITIVE
25 REM **** CAUTION: GWBASIC expects angles in radians
40 OPTION BASE 1 'array indexes start at 1
41 PI = 3.141593 'single precision
42 WVLTLENGTH = 256 'data points
43 SAMPRATE = 1 'milliseconds
44 SWEEPLENGTH = 1 'seconds
45 INITPHASE = 0 'radians
80 FFTLENGTH = WVLTLENGTH/2
83 L = LOG(FFTLENGTH) / LOG(2) 'number of cycles in FFT
84 M = FFTLENGTH 'abbreviated name for use in FFT
85 FREQINT = 1000 / (2^(L+1)) / SAMPRATE 'output frequency sample rate from FFT
86 DIM R [M], I[M], WVLT[WVLTLENGTH] 'Real and Imag arrays for FFT
4430 P1 = 0: IF FFTFLAG$ = "IFT" THEN P1 = 1: GOTO 4800 'Inverse Transform
4440 PRINT "FFT Loop 1": K = 0: FOR J = 1 TO M-1: I = 2 'FFT Starts Here
4455 IF K < M/I THEN GOTO 4470
4460 K = K - M / I: I = I + I: GOTO 4455
4470 K = K + M / I: IF K <= J THEN GOTO 4510
4480 A = R[J+1]: R[J+1] = R[K+1]: R[K+1] = A 'Binary Sort
4490 A = I[J+1]: I[J+1] = I[K+1]: I[K+1] = A
4510 NEXT J: G = .5: P = 1
4511 PRINT "FFT Loop 2": FOR I = 1 TO L: G = G + G: C = 1: E = 0
4520 Q = ((1-P)/2)^.5 * (1-2*P1)
4522 IF I = 1 THEN P = -1 ELSE P = ((1+P)/2)^.5
4525 FOR R = 1 TO G
4530 FOR J = R TO M STEP G+G: K = J + G
4540 A = C*R[K] + E*I[K]
4550 B = E*R[K] - C*I[K]
4560 R[K] = R[J] - A
4570 I[K] = I[J] + B
4580 R[J] = R[J] + A
4590 I[J] = I[J] - B
4600 NEXT J: A = E*P + C*Q: C = C*P - E*Q: E = A: NEXT R: NEXT I
4610 IF FFTFLAG$ = "IFT" THEN RETURN 'End of IFT
4801 REM **********************************************************************
4803 REM FAST FOURIER TRANSFORM PROGRAM - INVERSE
4805 REM ********* Author: E.R. Crain, P.Eng. 29 Oct 1987 *******************
4810 PRINT "FFT Loop 3" 'FFT Continues
4820 A = 4 * ATN(1) / M: P = COS(A) 'Start IFT
4830 Q = (1 - 2*P1) * SIN(A): A = R[1]: R[1] = A+I[1]: I[1] = A-I[1]
4835 IF FFTFLAG$ = "IFT" THEN GOTO 4850
4840 R[1] = R[1]/2: I[1] = I[1]/2
4850 C = 1 - 2*P1: E = 0: FOR J = 2 TO M/2: A = E*P + C*Q: C = C*P - E*Q:
E = A: K = M - J + 2
4855 A = R[J] + R[K]
4860 B = (I[J] + I[K])*C - (R[J] - R[K])*E: U = I[J] - I[K]
4865 V = (I[J] + I[K])*E + (R[J] - R[K])*C: R[J] = (A+B)/2: I[J] = (U-V)/2
4870 R[K] = (A-B)/2: I[K] = -(U+V)/2: NEXT J
4875 I[M/2+1] = -I[M/2+1]: IF FFTFLAG$ = "IFT" THEN GOTO 4440 'IFT Continues^
4880 FOR I = 1 TO M: R[I] = R[I]/128: I[I] = I[I]/128: NEXT I
4885 RETURN 'End of FFT
4901 REM **********************************************************************
4903 REM RECREATE WAVELET FROM PHASE AND AMPLITUDE
4905 REM ********* Author: E.R. Crain, P.Eng. 29 Oct 1987 *******************
4907 PRINT "Recreating Wavelet"
4910 FOR I = 1 TO FFTLENGTH: R[I] = A[I]*COS (P[I]*PI/180): I[I] = P[I]*SIN(P[I]*PI/180): NEXT I
4915 FFTFLAG$ = "IFT": GOSUB 4400 'Do IFT on Wavelet
4920 PRINT "FFT Loop 4R": RECSCALE = 5000*SAMPRATE
4925 FOR I = 1 TO FFTLENGTH: WVLT[2*I-1] = -I[I]/RECSCALE: NEXT I
4926 FOR I = 1 TO FFTLENGTH: WVLT[2*I] = -R[I]/RECSCALE: NEXT I
4927 WVLT[1] = 0
4930 PRINT "Plotting Recreated Wavelet"
4935 VIEW (300,10)-(500,110)
4940 WINDOW (-20,-2)-(128,2)
4945 LINE (0,0)-(0,0),0
4950 FOR I = 1 TO LENGTH: LINE -((I-1) * SAMPRATE, WVLT[I]),7:NEXT I
5000 REM **********************************************************************