%% Model of seismic source waveform .
%% In this illustration, the waveform is obtained by summation of multiple cosine functions
%% This summation is called the Fourier transform.
%% ALthough this transform can be performed by simply calling "fft" 
%% (Fast fourier transform, or FFT), we perform it here in the "slow" way - 
%% by direct summation of cosine harmonics.
%% In this way, we will see how the FFT works.

%% Parameters:
%%  t - array of times in milliseconds at which the values of the waveform are output 
%%      (these times should normally contain time t0)

%% Return values:
%%  wave - array of wavelet values at times given in 't'

function [wave] = wavelet_ft(t)
  
  freq = 0.02:0.0005:0.06;   % frequencies of cosine functions from 20 to 60 Hz
         
    % initialize the output array with zeros 
         
  wave = zeros(1,length(t)); 
	 
    % sum cosine functions to obtain an impulsive shape (symmetric in time)
% #	  x - source-receiver distance (can be an scalar or array, 

  for k=1:length(freq)
    wave = wave + cos(2*pi*freq(k)*t );
  end
   
    % scale the resulting waveform so that its amplitude at t=0 equals 1.0
    
  wave = wave/length(freq);
  
    % multiply the waveform by half-period sine function 
    % to ensure it equals zeros at both ends of t
    
  a = pi/(t(end) - t(1));   % factor equal pi when multiplied by the time range
  
  wave = wave.*sin(a*(t-t(1)));
  
end