% Walrus pulse analysis
% Magnus Wahlberg Apr 2015
% Script edited to include a normalization step for the entire file (CR/JM July 2020)
% Script edited to add 90% energy and duration (CR/JS July 2020)
% Script edited to add  export to .csv function (CR/JS July 2020)
% Be sure to add the correct sampling rate on line 12!

clear all
close all

file = 'walrus clap pulses normalized.wav'; % enter file name here
fs = 96e3; %confirm sampling rate needs to be entered here for each file
sig = audioread(file);
sig = sig/max(abs(sig));

thr = 0.4; % set threshold value to detect pulses and avoid outliers

i=1;
c = 0;
jump = 3e3;
while i< length(sig)-jump-1,
    if sig(i) > thr,
        c=c+1;
        % next line edited from bank(:,c) = sig(i-500:i+1999);
        bank(:,c) = sig(i-500:i+1999)/max(abs(sig(i-500:i+1999)));
        samp(c) = i;
        i=i+jump;
    end
    i=i+1;
end
ICI = diff(samp)/fs;
nc=c;

for c=1:nc,
   en = cumsum(bank(:,c).^2);
   en = 100*en/max(en);
   mi = min(find(en>2.5)); % 95% energy content from 2.5-97.5 %
   ma = min(find(en>97.5)); % 95% energy content from 2.5-97.5 %
   E(:,c) = en;
   En(c) = 20*log10( sqrt( sum(bank(mi:ma,c).^2)));
   
   mi2 = min(find(en>5)); % 90% energy content from 5-95 %
   ma2 = min(find(en>95)); % 90% energy content from 5-95 %
   E(:,c) = en;
   En90(c) = 20*log10( sqrt( sum(bank(mi2:ma2,c).^2)));
   % JS added 90% energy calculations 7/1/2020
   
   dur95(c) = (ma-mi)*1e3/fs;
   dur90(c) = (ma2-mi2)*1e3/fs;
   % added the 90% duration variable 7/20/2020
   hbank(:,c)=abs(hilbert(bank(:,c)));
   peak = min(find(hbank(:,c)==max(hbank(:,c))));
   peakvalue(c)=hbank(peak,c);
   vec10 = find(hbank(:,c)>max(hbank(:,c))*10^(-10/20));
   dur10(c) = (max(vec10)-min(vec10))*1e3/fs;
   vec3 = find(hbank(:,c)>max(hbank(:,c))*10^(-3/20));
   dur3(c) = (max(vec3)-min(vec3))*1e3/fs;
   t=[0:length(bank)-1]'*1e3/fs;
   t0 = sum(t.*bank(:,c).^2)/sum(bank(:,c).^2);
   trms(c) = sqrt (sum( (t-t0).^2.*bank(:,c).^2)/sum(bank(:,c).^2) );
   f = abs(fft(bank(:,c).*hanning(length(bank))))';
   f = f(1:length(f)/2)/max(f);
   fbank(:,c) = 20*log10(f);
   fx = [0:length(fbank)-1]*fs / (2e3 * length(fbank));
   pf(c) = fx(min(find(fbank(:,c)==max(fbank(:,c)))));
   bw3int = fx(find(fbank(:,c)>-3));
   bw3(c) = max(bw3int)-min(bw3int);
   bw10int = fx(find(fbank(:,c)>-10));
   bw10(c) = max(bw10int)-min(bw10int);
   centf(c) = sum(fx.*f.^2) / sum(f.^2);
   rmsbw(c) = sqrt(sum((fx-centf(c)).^2.*f.^2)/sum(f.^2));
end

template = find(peakvalue==max(peakvalue));
for i=1:max(c),
    corrcoeff(i) = max(xcorr(bank(:,i),bank(:,template)))/max(xcorr(bank(:,template),bank(:,template)));
% figure(2); plot(xcorr(bank(:,i),bank(:,template)));
% pause
end

subplot(3,1,1); plot(t,bank)
subplot(3,1,2); plot(t,E)
ylabel('Cumulative energy (%)')
xlabel('Time(ms)')



subplot(3,1,3); plot(fx,fbank)
set(1,'Color',[1 1 1])
disp(['Inter-signal interval: ' num2str(ICI) ' s'])
disp(['95% energy duration: ' num2str(dur95) ' ms'])
disp(['90% energy duration: ' num2str(dur90) ' ms'])
disp(['-3 dB duration: ' num2str(dur3) ' ms'])
disp(['-10 dB duration: ' num2str(dur10) ' ms'])
disp(['RMS duration: ' num2str(trms) ' ms'])
disp(['95% Energy: ' num2str(En) ' dB re 1 µPa2s'])
disp(['90% Energy: ' num2str(En90) ' dB re 1 µPa2s'])
% added variable on 7/1/2020
disp(['Peak frequency: ' num2str(pf) ' kHz'])
disp(['Centroid frequency: ' num2str(centf) ' kHz'])
disp(['RMS BW: ' num2str(rmsbw) ' kHz'])
disp(['-3dB Band width: ' num2str(bw3) ' kHz'])
disp(['-10dB Band width: ' num2str(bw10) ' kHz'])

output = [dur95; dur90; dur3; dur10; trms; En; En90; pf; centf; rmsbw; bw3; bw10];
names={'95%_energy_duration_(dur95)', '90%_energy_duration_(dur90)', '-3dB_duration_(dur3)', '-10dB_duration_(dur10)', 'RMS_duration_(trms)', '95%_Energy_(En)', '90%_Energy_(En90)', 'Peak_Frequency_(pf)', 'Centroid_Frequency_(centf)', 'RMS_BW_(rmsbw)', '-3dB_Bandwidth_(bw3)', '-10_dB_Bandwidth_(bw10)', 'Inter-signal_interval_(ICI)'};
delete output.csv
FID=fopen('output.csv','w');
fprintf(FID,'You are so smart!\n');
for i=1:size(output,1)
    temp=[names{i} ' ' num2str(output(i,:))];
    fprintf(FID,'%s ',names{i});
    fprintf(FID,'%f ',output(i,:));
    fprintf(FID,'\n');
end

names={'Inter-signal_interval_(ICI)'};

for i=1:size(ICI,1)
    temp=[names{i} ' ' num2str(ICI(i,:))];
    fprintf(FID,'%s ',names{i});
    fprintf(FID,'%f ',ICI(i,:));
    fprintf(FID,'\n');
end
fclose(FID);

% remember that the output file will get overwritten each time you
% run the code unless you save it and rename it.


