function [cloc] = CD_ALG(data,TH)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Algorithm to detect Back-to-Zero (BZ) type of clipped points on a given
% waveform data.
% Usage:
% [cloc] = CD_ALG(data,TH)
%  data:     Waveform data with or without BZ type of clipped points
%  TH  :     Threshold (e.g. 60% of the Observed Range)
%  cloc:     The locations of BZ type of clipped points on the data.
%            If cloc is empty, there is no BZ type of clipped point 
%            detected on the given waveform data.
%
%              Wenzheng Yang and Yehuda Ben-Zion, 2009
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If no threshold value is given, use a 60% of the maximum data amplitude
% as the threshold
if nargin <2 
  TH       = 0.6*max(data);
end

% Taper the first and the last N_tape points on the waveform data
N_tape     = 10;
temp       = data;
temp_begin = N_tape +1;
temp_end   = length(temp) - N_tape;
data1      = temp(temp_begin:temp_end);

% Criteria 1. Treat waveform data with maximum amplitude larger
%             than threshold as potentially clipped waveform.
%
if max(abs(data1))>TH
 criteria1 = 1;
else
 criteria1 = 0;   
end

% Criteria 2: For the potential clipped waveform, assume that points
%             with zero amplitude, between the first and last samples
%             with absolute value larger than 0.5 of the waveform peak
%             amplitude, are tentatively clipped.
%
clip_point = find(data1 ==0);
clip_point = clip_point + N_tape;

temp       = find((data>(max(data)*0.5)) | (data < (min(data)*0.5)));
if length(temp) >= 2
 criteria2 = clip_point>temp(1) & clip_point< temp(end);
else
 criteria2 = (clip_point==clip_point);
end

% Criteria 3. Require the clipped samples are bounded on both sides by
%             points with the same sign or another zero.
%
criteria3  = (sign(data(clip_point+1))+sign(data(clip_point-1)))~=0;

% Criteria 4. Require that the maximum or mininmum values of +/- 10
%             samples on either side of the candidate clipped points is
%             larger than 0.8 of the waveform peak amplitude.
%
N_clip     = length(clip_point);
criteria4  = zeros(N_clip,1);
for i      = 1:N_clip
    temp   =  data((clip_point(i)-N_tape):(clip_point(i)+N_tape));
    if max(temp)>0.8*max(data) | min(temp) < 0.8*min(data)
       criteria4(i)  = 1;
    end
end

cloc       = clip_point(criteria1 & criteria2 & ...
                        criteria3 & criteria4);

end