[SP_maxa, index_maxa]=peak_index(SP, center_bin, pts_per_bin, num_pts_per_pk_intrl)

function [SP_maxa, index_maxa]=peak_index(SP, center_bin, pts_per_bin, num_pts_per_pk_intrl)
% % peak_index: Finds the highest data point centered about center_bin
% %
% % Syntax:
% %
% % [SP_maxa, index_maxa]=peak_index(SP, center_bin, pts_per_bin, num_pts_per_pk_intrl);
% % 
% % ********************************************************************
% %
% % Description
% %
% % This program finds the value adn index of the maximum absolute value 
% % of the data array SP within pts_per_bin number of data points of the
% % center_bin.  The center_bin is the bin number of the peak.  
% % 
% % 
% % ********************************************************************
% %
% % Input Variables
% % 
% % SP is the data array
% % 
% % center_bin is the bin where the center of the impulese occurs. 
% % 
% % pts_per_bin is the number of points in each bin
% % 
% % num_pts_per_pk_intrl is the number of datat points in each peak interval 
% %  
% % 
% % ********************************************************************
% % 
% % Output Variables
% % 
% % SP_maxa is an array of the absolute values of the 
% % amplitude of each peak.   
% %
% % index_maxa is an array of the indices of the peaks.
% % 
% % ********************************************************************
% 
% 
% Example='1';
% 
% SP=randn(1, 100000);  % input data array one-dimensional
% peak_ix=20250;        % choose an index for the impulsive peak
% SP(peak_ix)=100;      % for this example create a point impulse
% pts_per_bin=500;      % the pts_per_bin is used to find the highest peak 
%                       % within a range of indices fo the data array
% center_bin=ceil(peak_ix/pts_per_bin);
%                       % the center_bin is the bin number for choosing the 
%                       % range of indices to search
%                       % center_bin can be a one-dimensional array
%                       
% [SP_maxa index_maxa]=peak_index(SP, center_bin, pts_per_bin);
% 
% 
% 
% % ********************************************************************
% %
% % peak_index is written by Edward L. Zechmann 
% %    
% %     date 26 December    2007
% %    
% % modified 27 December    2007    Added comments and an example
% %    
% % modified 19 September   2008    Updated Comments
% %   
% % modified 11 October     2010    Improved Error handling.  
% %                                 Updated Comments
% %
% % modified 11 February    2011    Updated to new method of finding peaks.
% %                                 Updated Comments.
% %     
% %     
% %     
% % ********************************************************************
% % 
% % Please feel free to modify this code.
% % 

if nargin < 1 || isempty(SP) || ~isnumeric(SP)
    SP=randn(4, 100000);
end

[m1 n1]=size(SP);

if m1 > n1
    SP=SP';
    [m1 n1]=size(SP);
end

if nargin < 2 || isempty(center_bin) || ~isnumeric(center_bin)
    
    center_bin=zeros(m1, 1);
    for e1=1:m1;
        [buf buf_ix]=max(SP(e1, :));
        center_bin(e1,1)=buf_ix;
    end
    
end


if nargin < 3 || isempty(pts_per_bin) || ~isnumeric(pts_per_bin)
    pts_per_bin=500;
end

if nargin < 4 || isempty(num_pts_per_pk_intrl) || ~isnumeric(num_pts_per_pk_intrl)
    num_pts_per_pk_intrl=pts_per_bin;
end




nc1=length(center_bin);

SP_maxa=zeros(nc1, 1);
index_maxa=zeros(nc1, 1);

low_bin=floor(num_pts_per_pk_intrl);
high_bin=floor(num_pts_per_pk_intrl);

for e1=1:nc1;

    index1=pts_per_bin*(center_bin(e1)-1)+1-low_bin;
    index2=pts_per_bin*(center_bin(e1)-1)+1+high_bin;

    if index1 < 1
        index1=1;
        index2=index1+low_bin+high_bin;
    end

    if index2 > n1
        index2=n1;
    end

    % index3 is the local index of the peak
    [SP_max index3]=max(abs(SP(1, index1:index2  )));
    if ~isempty(SP_max)
        SP_maxa(e1)=SP_max;
        index_maxa(e1)=index1-1+index3;
    end

end

index_maxa=index_maxa(index_maxa>=1);