gp=peak_threshhold_function2(SP, Fs, pts_per_bin, no_smoothing)

function gp=peak_threshhold_function2(SP, Fs, pts_per_bin, no_smoothing)
% % peak_threshhold_function2: Calculates a threshold for finding peaks
% %
% % Syntax:
% %
% % gp=peak_threshhold_function2(SP, Fs, pts_per_bin, no_smoothing);
% %
% % ********************************************************************
% %
% % Description
% %
% % In complex noise, continuous noise and impulsive noise are combined.
% % To make it easier to identify impulsive peaks.  The time record is
% % processed in a few steps.
% %
% % Step  Description
% %
% % 1)    Break up time record into bins and select the peak from each bin
% % 2)    Apply smoothing.
% % 3)    Subtract the mean value.
% %
% % This process makes it easier to apply a threshold to identify peaks.
% % The impulses are typically significantly greater than the
% % remaining noise.
% %
% %
% %
% % ********************************************************************
% %
% % Input Variables
% %
% % SP is the array of sound pressure data in Pa.
% %      should have size [num channels length_time_record]
% %      assumes the length of the time record is greater than the number
% %      of channels and tranposes SP if necessary
% %      default value is randn(4, 100000).
% %
% % Fs sampling rate in Hz.  default value is 100000 Hz.
% %
% % pts_per_bin is the number of points in the sub bins for speedings up
% %      the process of analyzing the data array.  The pts_per_bin is set
% %      so that there are 100 bins peak peak interval.
% %      The default value is pts_per_bin=floor(Fs/100);.
% %
% % no_smoothing=0;     % 1 does not apply any smothign to the threshold
% %                     % function.
% %                     %
% %                     % Otherwise smoothign is applied to the function
% %                     %
% %                     % The default value is no_smoothing=0;
% %
% % ********************************************************************
% %
% % Output Variables
% %
% % gp is a matrix of size [num_channels pts_per_bin].  gp is formed by
% %      processing SP.  The purpose is to maximize the likelyhood of
% %      detecting a peak while minimizing teh likelihood of falsely
% %      detecting a peak.
% %      gp makes identifying impulsive peaks much easier than
% %      applying a threshold on SP.
% %      The highest values of gp correspond to impulsive noise peaks.
% %
% % ********************************************************************
%
%
% Example
%
% Fs=100000; fc=1000; td=1; tau=0.01; delay=0.1; A1=5; A2=20;
% [SP, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
% pts_per_bin=floor(Fs/100);
%
% gp=peak_threshhold_function(SP, Fs, pts_per_bin);
%
%
% % ********************************************************************
% %
% % Subprograms
% %
% % 
% % List of Dependent Subprograms for 
% % peak_threshhold_function2
% % 
% % FEX ID# is the File ID on the Matlab Central File Exchange
% % 
% % 
% % Program Name   Author   FEX ID#
% % 1) convert_double		Edward L. Zechmann			
% % 2) estimatenoise		John D'Errico		16683	
% % 3) moving		Aslak Grinsted		8251	
% % 4) sub_mean		Edward L. Zechmann			
% % 5) wsmooth		Damien Garcia		NA	
% %
% %
% % ********************************************************************
% %
% % peak_threshhold_function is written by Edward L. Zechmann
% %
% %     date     August     2008
% %
% % modified  16 August     2008    Updated Comments
% %
% % modified 10 September   2008    Added exhaust_cycle finding option
% %                                 Updated Comemnts
% %
% % modified 11 November    2008    Fixed a bug with the exhaust noise
% %                                 feature and with bin size.
% %
% % modified 14 April       2009    Alterred algorithm greatly sped up
% %                                 progam.
% %
% % modified  6 October     2009    Updated comments
% %
% % modified  8 October     2010    Removed exhaust cycle option.
% %
% % 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(4, 100000);
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.
[m1 n1]=size(SP);

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

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

[Fs]=convert_double(Fs);


if nargin < 3 || isempty(Fs) || ~isnumeric(Fs)
    pts_per_bin=floor(Fs/100);

    if pts_per_bin < 1
        pts_per_bin=1;
    end

    if pts_per_bin > n1
        pts_per_bin=n1;
    end
end

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




% % ********************************************************************
%
% STEP 1 Break up time record into bins and select the peak from each bin
%
% Now the peaks are of primary interest.
%
% Break up the processed time record into 100 bins per second
% and calculate some statistics for each bin.
% The most important metric is the maximum peak in each bin.
n_bins=floor(n1/pts_per_bin);
gp=zeros(m1, n_bins);

for e1=1:n_bins;
    for e2=1:m1;

        gp( e2, e1)=max(abs(SP( e2, (e1-1)*pts_per_bin+(1:pts_per_bin) )));

    end
end



gp(gp<0)=0;


% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% STEP 2 Apply smoothing.
%

if ~isequal(no_smoothing, 1)
    [buf, gp]=sub_mean(abs((gp')'), Fs/pts_per_bin, 20);
end

for e1=1:m1;
    [gp(e1, :)] = wsmooth(gp(e1, :), 2);
end

% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% STEP 3 Subtract Mean Value.
%

gp=gp-min(gp, [], 2)*ones(1, size(gp, 2));