% 06/02/19 This script will loop through clipped wav files each containing 
%an echolocation click, and extract the instantaneous frequency and the
%duration based upon an established criterion of amplitude difference from 
%the peak of the signal.

%Path Input of directory of wav files
pathIn = 'C:\Program Files\MATLAB\R2018b\EchoPulseIn\';
%Set the path for output figures of the signal plus envelope, exponential curve, and start/end
%points
pathOut = 'C:\Program Files\MATLAB\R2018b\EchoPulseOut\';
myFiles = dir(fullfile(pathIn,'*.wav')); %gets all wav files in struct

%Preallocate memory for extracted parameters
n = length(myFiles);
%duration of signal
duration = nan(n,1);
%Steepness of the wavelet from peak to end duration (if desired)
Steepness = nan(n,1);
%Frequency at the peak amplitude of the click wavelet
PeakF = nan(n,1);
%Filename of the click that was analyzed
fname = cell(n,1);
%Voltage at the peak of the click wavelet
PeakV = nan(n,1);

%Loop through the wav files each containinng a click
for k = 1:n
   clear figure
   filename = [pathIn myFiles(k).name];
   [~, fname{k}, ext] = fileparts(filename);
   [y, fs] = audioread(filename);
   %Make sure it's the appropriate channel. Use y(:,3) for 3-channel
   %recordings, just y for normal single channel
   sig = y(:,1);
   %resample fix for 500 khz sampled bats. Comment out if fs from audioread
   %is sufficient.
   sig = resample(sig,1,2);
   fs = 250000;
   
   t = (0:length(sig)-1)/fs;
   t_vert = t.';

    %%%Get the envelope of the signal. The envelope func creates the same as the hilbert code piece below%%
    [up lo] = envelope(sig,100,'analytic');

    %Get the first three peaks of the signal from the envelope
    [pks, locs, w, p] = findpeaks(up,t_vert,'MinPeakHeight',.02,'MinPeakProminence',.001);
    figure(k); plot(t_vert,sig, '-o'); hold on; plot(t_vert,up); hold on; plot(locs, pks,'o');

    %findmaxpeak
    [M I] = max(pks);

    %Create vector of top peak and the two subsequent peaks
    toppeaks = [M, pks(I+1), pks(I+2)]; 
    %peak location vector of the three peaks
    toplocs = [locs(I), locs(I+1), locs(I+2)];
    %Fit exponential curve to the points and extract the steepness coeff. b
    f = fit(toplocs.',toppeaks.','exp1');
    plot(f,toplocs.',toppeaks.')
    Steepness(k) = f.b;

    %Set amplitude criterion to calculate duration. E.g. -10 is -10 db from either side of the
    %peak of the envelope of the signal to get start and end of the signal
    criterion = -10; 
    %Calculate what the Y value should be using the -10 criterion
    CriticalY = (10^(criterion/20))*M;
    %Use exponential equation and function coeffiecients from above to extract time for given y
    %A decreasing exponential equation is used for the slope from the peak
    endDur = (log(CriticalY/f.a))/f.b;
    plot(endDur,CriticalY,'*')

    %get the index of the peak point on the envelope array
    Mloc = find(up == M); 
    %Get first half of the envelope array
    firsthalf = up(1:Mloc);

    %Target the closest possible array point to CriticalY value from the
    %beginning of the signal to the peak of the envelope
    [c index] = min(abs(firsthalf - CriticalY));

    %Use interpolation to calulate x given Critical Y, and closest Y points
    %above and below
    if c < CriticalY
        startdur = t_vert(index)+(t_vert(index+1)- t_vert(index))*((CriticalY-up(index))/(up(index+1)-up(index)));
   
    elseif c > CriticalY
        startdur = t_vert(index)+(t_vert(index-1)+ t_vert(index))*((CriticalY-up(index-1))/(up(index)-up(index-1)));
 
    else 
        disp('Incorrect Index.')
    end
 
    plot(startdur,CriticalY, '*')

    duration(k) = endDur-startdur;

    fig = figure(k);
    saveas(fig, [pathOut fname{k}], 'jpeg')
    
    %Using more traditional instantaneous frequency methods
    %Find peak of the raw signal
    [Sigpks, Siglocs, Sigw, Sigp] = findpeaks(sig,t,'MinPeakHeight',.02,'MinPeakProminence',.001);
    [SigM SigI] = max(Sigpks);
    PeakV(k) = SigM;
    Peakloc = find(sig == SigM);
    hilsig = hilbert(sig);
    instfrq = diff(unwrap(angle(hilsig)))/((1/fs)*2*pi);
    %Get the frequency at the peak of the signal
    PeakF(k) = instfrq(Peakloc);
    fig2 = figure(k+n);
    
end
