[peak_indices, gp, Pa_peak, Lp_peak]=peak_ix(SP, Fs, threshold, make_plot, ignore_threshold)

function [peak_indices, gp, Pa_peak, Lp_peak]=peak_ix(SP, Fs, threshold, make_plot, ignore_threshold)
% % peak_ix: Locates peaks for acoustic impulsive noise
% %
% % Syntax:
% %
% % [peak_indices, Pa_peak, Lp_peak]=peak_ix(SP, Fs, threshold, make_plot, ignore_threshold);
% %
% % ********************************************************************
% %
% % Description
% %
% % Locates peaks for acoustic impulsive noise
% % 
% % 
% % 
% % 
% % ********************************************************************
% %
% % Input Variables
% %
% % SP=randn(100000, 4);% Is the array of sound pressure data in Pa.
% %                     % should have size [num_samples num channels ];
% %                     % assumes the length of the tiem record is greater 
% %                     % than the number of channels and tranposes SP 
% %                     % if necessary
% %                     % default value is randn(100000, 4); 
% %
% % Fs=100000;          % sampling rate in Hz.  
% %                     % default value is Fs=100000; % Hz.
% %
% % threshold=100;      % 
% %                     % default value is threshold=100; % dB.
% %
% % make_plot=1;        % 1 Plots the processing data of the peaks.  
% %                     % Otherwise no plots are made. 
% %                     % 
% %                     % default is make_plot=0;  
% %
% % ignore_threshold=0; % 1 surpresses the threshold and other requirements
% %                     % for selecting peaks
% %                     % ignore_threshold=1; % generally will find more peaks.
% %                     % 
% %                     % 0 is the default.  Returns typical number of
% %                     % peaks. 
% %                     % 
% %                     % default is ignore_threshold=0;
% %
% %
% % ********************************************************************
% % 
% % Output Variables
% % 
% % peak_indices is a cell array of indices of the peaks for each channel 
% %
% % gp is the array of threshold values for finding peaks.  
% %
% % Pa_peak (Pa) is a cell array of peak pressures  of the peaks for each channel 
% %
% % Lp_peak (dB) is a cell array of the peak levels  of the peaks for each channel 
% %
% % ********************************************************************
%
%
% Example='1';
%
% tau=0.01;         % seconds decay constant for modelling the ringing
% td=1;             % seconds time duration of the test signal
% A2=20;            % Pa initial ampitude of the of the sin wave
% A1=1;             % Pa amplitude of the background gaussian noise
% Fs=100000;        % Hz sampling rate of the impulsive noise
% fc=1000;          % Hz ringing frequency of sin wave
% delay=0.1;        % seconds time delay between impulsive noise 
%                   % and test signal recording
% rt=0.0005;        % seconds time to reach peak value from back ground
%                   % noise level
% 
% [SP, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2, rt);
% threshold=80;
% make_plot=1;
% ignore_threshold=0;
% [peak_indices, Pa_peak, Lp_peak]=peak_ix(SP, Fs, threshold, make_plot, ignore_threshold);
% 
% 
% 
% % ********************************************************************
% %
% % Subprograms
% %
% % 
% % 
% % List of Dependent Subprograms for 
% % peak_ix
% % 
% % FEX ID# is the File ID on the Matlab Central File Exchange
% % 
% % 
% % Program Name   Author   FEX ID#
% %  1) convert_double		Edward L. Zechmann			
% %  2) dB_to_Pa		Edward L. Zechmann			
% %  3) estimatenoise		John D'Errico		16683	
% %  4) findjobj		Yair M. Altman		14317	
% %  5) func_threshold		Jing Tian		10462	
% %  6) LMSloc		Alexandros Leontitsis		801	
% %  7) maximize		Oliver Woodford		25471	
% %  8) moving		Aslak Grinsted		8251	
% %  9) Pa_to_dB		Edward L. Zechmann			
% % 10) peak_threshhold_function2		Edward L. Zechmann			
% % 11) peakfinder		Nate Yoder		25500	
% % 12) sub_mean		Edward L. Zechmann			
% % 13) wsmooth		Damien Garcia		NA	
% %
% %
% % 
% % ********************************************************************
% %
% % peak_threshhold_function is written by Edward L. Zechmann
% %
% %     date  8 October     2010
% %
% % modified 11 October     2010    Updated Comments
% %
% % modified 11 February    2011    Updated to new method of finding peaks.
% %                                 Updated Comments.
% %
% %
% %
% %
% %
% % ********************************************************************
% %
% % Please feel free to modify this code.
% %
% % See Also: Impulsive_Noise_Meter, plot_peaks, local_peaks
% %




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

% make sure that input data is double precision.
[SP]=convert_double(SP);

% Make sure the data has the channels and data array in the
% correct order, transpose if necessary.
[num_samples num_chans]=size(SP);

if num_samples < num_chans
    SP=SP';
    [num_samples num_chans]=size(SP);
end

if nargin < 2 || isempty(Fs) || ~isnumeric(Fs)
    Fs=100000;
end

[Fs]=convert_double(Fs);


if (nargin < 3 || isempty(threshold)) || ~isnumeric(threshold)
    threshold=100;
end

if (nargin < 4 || isempty(make_plot)) || ~isnumeric(make_plot)
    make_plot=0;
end

if (nargin < 5 || isempty(ignore_threshold)) || ~isnumeric(ignore_threshold)
    ignore_threshold=0;
end


Lt=dB_to_Pa(threshold);

% %
pts_per_peak=ceil(Fs/1000);



% % compute a threshold function for impoulsive noise
[gp]=peak_threshhold_function2(SP, Fs, pts_per_peak, ignore_threshold);

peak_indices=cell(num_chans, 1);
Pa_peak=cell(num_chans, 1);
Lp_peak=cell(num_chans, 1);


for e1=1:num_chans;

    if isequal(ignore_threshold, 1)

        bb2=LMSloc(gp(e1, :));

        if length(bb2 > gp(e1, :)) > ceil(0.01*num_samples)
            bb2=LMSloc(gp( e1, gp(e1, :) < bb2));
        end
        
    else
        
        bb2=LMSloc(gp(e1, :));
        
        if length(bb2 > gp(e1, :)) > ceil(0.1*num_samples)
            bb2=LMSloc(gp( e1, gp(e1, :) < bb2));
        end

    end
    
    
    [peak_loc]=peakfinder(gp(e1, :), bb2*1.1, +1);
    
    if ~isequal(ignore_threshold, 1)
        thre = func_threshold(gp( e1,:));
    
        peak_loc=peak_loc(gp(e1, peak_loc) > thre);
    end
    
    num_peaks=length(peak_loc);
    peak_ix_buf=zeros(num_peaks, 1);

    for e2=1:num_peaks;
        ix=floor(pts_per_peak*(peak_loc(e2)-30)):ceil(pts_per_peak*(peak_loc(e2)+30));
        ix(ix<1)=1;
        ix(ix>num_samples)=num_samples;
        [buf, buf_ix]=max(abs(SP(ix, e1)));
        peak_ix_buf(e2)=ix(1)-1+buf_ix(1);
    end

    peak_loc=unique(peak_ix_buf);
    
    if ~isequal(ignore_threshold, 1)
        peak_loc=peak_loc((abs(SP(peak_loc, e1)) > Lt)); 
    end
    
    peak_indices{e1, 1}=peak_loc;
    
    % Calculate peak pressure and peak levels
    Pa_peak{e1, 1}=SP(peak_loc, e1);
    Lp_peak{e1, 1}=Pa_to_dB(SP(peak_loc, e1));

end



if isequal(make_plot, 1)

    figure(1);
    delete(1);
    figure(1);
    for e1=1:num_chans;

        peak_loc=peak_indices{e1, 1};

        subplot(num_chans, 1, e1);
        plot((-1+(1:num_samples))/Fs, SP(: , e1), '-b', 'linewidth', 1);
        hold on;
        plot((peak_loc-1)./Fs, SP(peak_loc , e1), 'linestyle', 'none', 'marker', 'o', 'markeredgecolor', 'k');

    end

    fh1=figure(1);
    maximize(fh1);
        
end