linear_events.m

著名的seismiclab的代码 是地震学研究人员必备的工具

function [d] = linear_events(dt,f0,tmax,h,tau,p,amp);
%LINEAR_EVENTS: A program to generate data containing linear
%               events.
%
%  [d] = linear_events(dt,f0,tmax,h,tau,p,amp);
%
%  IN   dt: sampling interval in secs
%       f0: central freq. of a Ricker wavelet in Hz
%       tmax: maximun time of the simulation in secs
%       h: vector of desire offsets in meters
%       tau,p,amp: vectors of intercept, ray paramter 
%                  and amplitude of each linear event
%                  (p is in sec/m and tau in secs)
%
%  OUT  d: data that consist of a superposition of linear
%          events
%
%  Example: Generation of two linear events 
%
%           dt = 4./1000; f0=20; tmax=1.; 
%           h=1:5:100; tau=[0.1,0.6]; p=[0.01,-0.005]; amp=[2.,-1];
%           [d] = linear_events(dt,f0,tmax,h,tau,p,amp);
%           wigb(d);       
%
%  SeismicLab
%  Version 1
%
%  written by M.D.Sacchi
%  [e-mail]
%
%  Copyright (C) 2002 Signal Analysis and Imaging Group
%                     Department of Physics
%                     The University of Alberta
%

 nt = floor(tmax/dt)+1;

 nfft = 2*(2^nextpow2(nt));

 n_events = length(tau);

 nh = length(h);

 wavelet = ricker(f0,dt); 

 nw = length(wavelet);

 W = fft(wavelet,nfft);

 D = zeros(nfft,nh);
 i = sqrt(-1);

 delay = dt*(floor(nw/2)+1);

 for ifreq=1:nfft/2+1
  w = 2.*pi*(ifreq-1)/nfft/dt;
   for k=1:n_events
    Shift = exp(-i*w*(tau(k)+h*p(k)-delay));
   D(ifreq,:) = D(ifreq,:) +amp(k)* W(ifreq)*Shift;
  end
 end

 for ifreq=2:nfft/2
  D(nfft+2-ifreq,:) = conj(D(ifreq,:));
 end 

 d = ifft(D,[],1);
 d = real(d(1:nt,:));

 return;