%% PORCC ALGORITHM AND FUNCTIONS
Fs = 500000; % Sampling frequency
FFT = 512; %FFT value for fft
Th1 = 0.9999; % threshold for HQ clicks
Th2 = 0.55; % threshold for LQ clicks
% Coefficients for Model 1 (HQ)
logitCoefHQ = [-12.857529903901488;1.480233893649301;-0.037748903932102];
% Coefficients for Model 2 (LQ)
logitCoefLQ = [4.359888452493062;0.873710953710921;-0.017386257791731;...
  -6.901099488779461;35.272929893479020;-1.443845172568886e-05;-0.142569799845068];

click_model = [0.00336253736901471;0.00336253736901471;-0.0201752242140883;...
    0;0.0437129857971912;-0.0302628363211324;-0.0672507473802942;...
    0.114326270546500;0.0336253736901471;-0.201752242140883;0.0941510463324119;...
    0.211839854247927;-0.238740153200044;-0.124413882653544;0.346341349008515;...
    -0.0235377615831030;-0.359791498484574;0.171489405819750;0.302628363211324;...
    -0.272365526890192;-0.228652541093000;0.353066423746545;0.137864032129603;...
    -0.396779409543736;-0.0403504484281765;0.413592096388809;-0.0538005979042354;...
    -0.427042245864868;0.178214480557780;0.386691797436692;-0.299265825842309;...
    -0.285815676366250;0.393416872174721;0.137864032129603;-0.400141946912750;...
    0.00336253736901471;0.346341349008515;-0.121051345284530;-0.232015078462015;...
    0.147951644236647;0.134501494760588;-0.141226569498618;-0.0605256726422648;...
    0.104238658439456;0.0302628363211324;-0.0672507473802942;-0.0269002989521177;...
    0.0504380605352207;0.0437129857971912;-0.0807008968563530;-0.00672507473802942;...
    0.0773383594873383;-0.0201752242140883;-0.0638882100112795;0.0504380605352207;...
    0.0369879110591618;-0.0672507473802942;-0.0100876121070441;0.0706132847493089;...
    -0.0168126868450736;-0.0571631352732501;0.0235377615831030;0.0538005979042354;...
    -0.0470755231662059;-0.0403504484281765;0.0504380605352207;0.0269002989521177;...
    -0.0571631352732501;-0.0100876121070441;0.0605256726422648;-0.00336253736901471;...
    -0.0437129857971912;0.0201752242140883;0.0235377615831030;-0.0168126868450736;...
    -0.0100876121070441;0.00336253736901471;0;0.0168126868450736;-0.0134501494760588;...
    -0.0269002989521177;0.0336253736901471;0.0302628363211324;-0.0571631352732501;...
    -0.0134501494760588;0.0807008968563530;-0.0168126868450736;-0.0807008968563530;...
    0.0538005979042354;0.0605256726422648;-0.0773383594873383;-0.0269002989521177;...
    0.0840634342253677;-0.0100876121070441;-0.0739758221183236;0.0369879110591618;...
    0.0470755231662059;-0.0470755231662059;-0.0168126868450736;0.0403504484281765;0;...
    -0.0269002989521177;0.00672507473802942;0.0134501494760588;-0.00672507473802942;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;...
    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];



for i=1:length(NewRandomAndPAM) %each pgdf file
  click = NewRandomAndPAM(i).wave; % the wave as extracted from PAMGUARD
  [PSD, f, CF, PF] = CentroidAndPeakFreq(click, FFT, Fs);
  [RMSBW] =  RMSBWD(PSD, f, CF); 
  Q = (CF/RMSBW)/1000; 
    if PF > 100000 && PF < 160000 && CF > 100000 && CF < 160000 && Q > 4
        % Møhl & Andersen, 1973
     duration = ClickDuration(click); 
     ratio = PF/CF;
     [XCoeff, ~] = xcorr(click, click_model); % click_model is below
        XC = max(XCoeff);
     BW =  powerbw(click, Fs)/1000;  
       
    % estimate probabilities for HQ and LQ
    ProbHQ = glmval(logitCoefHQ,[Q duration],'logit'); % Log Model 1
             % assign each clip to a category
            if ProbHQ >= Th1 
               Porps = 1;  % HQ click
            else 
        ProbLQ = glmval(logitCoefLQ,[Q duration ratio XC CF BW],'logit'); % Log Model 2
                if ProbLQ > Th2
                   Porps = 2; % LQ click
                else 
                   Porps = 3; % HF Noise
                end 
            end
    else 
     Porps = 3;
   end % if CF, PF, Q
end  % for loop

%% Functions 
function [PSD, f, CF, PF] = CentroidAndPeakFreq(click, FFT, Fs)
[PSD,f] = periodogram(click,[],FFT,Fs,'power');
%psd=abs(fft(click))
PSD = PSD./max(PSD); % normalize spectrum
CF = sum(f.*PSD.^2)./sum(PSD.^2);
[~,indfc] = max(PSD);
PF = f(indfc); 
end

function [BW] = RMSBWD(PSD, f, CF)
BW = (sqrt(sum((f-CF).^2.*PSD.^2 ) / sum(PSD.^2)))/1000; 
end 

function [duration] = ClickDuration(click)
% Click duration based on Madsen & Walhberg 2007 - 80%
    Ener = cumsum(click.^2); 
    Ener2_5 = Ener(end) * 0.1; % 1.5% of acumulated energy
    Ener97_5 = Ener(end) * 0.9;% 98.5% of acumulated energy
    ix = sum(Ener<=Ener2_5)+1; % index of where the 1.5% is
    iy = sum(Ener<=Ener97_5)+1; % index of where the 98.5% is
    duration = (iy-ix)*2; %in microseconds
    %new_click1 = click1a(ix1:iy1);
end