spiking.m

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

function [f,o] = spiking(d,NF,mu);  
%SPIKING: Spiking deconvolution  using Levinson's recursion
%
%  [f,o] = spiking(d,NF,mu)
%
%  IN     d:  data (a column trace)
%         NF: lenght of the spiking operator
%         mu: prewhitening in percentage  
%
%  OUT    f:  the filter
%         o:  the ouput or convolution of the data with 
%             the filter (adjusted to the length of the
%             input data) 
%
%  Note: we assume a minimum phase wavelet
%
%  SeismicLab
%  Version 1
%
%  written by M.D.Sacchi, last modified December 10,  1998.
%  [e-mail] 
%
%  Copyright (C) 1998 Seismic Processing and Imaging Group
%                     Department of Physics
%                     The University of Alberta
%
%

NF = NF - 1;
R = xcorr(d,d,NF);

r = R(NF+1:2*NF+1);
r(1) = r(1)*(1+mu/100.);

[f]=levinson(r,NF);

if nargout == 2
o = conv(f,d);          
o = o(1:length(d)); %adjust lenght of output
end

return