load('trouble22467_waveforms_12-14DIV.mat');

for chanNum = 1:64

clearvars -except timestamps chanNum
% chanNum = 12;
numUnitsDisplay = 4;

close all;

for ii = 1:3
    if ii==1
        load('trouble22467_12-14DIV_LineSort.SS','-mat'); plotShift = 0; label = 'Manual Sorting'; figure('units','normalized','outerposition',[0 0 1 1]);
    elseif ii==2
        load('trouble22467_12-14DIV_WaveletSort.SS','-mat'); plotShift = numUnitsDisplay*1; label = 'Wavelet Transform';
    else
        load('trouble22467_12-14DIV_PCASort.SS','-mat'); plotShift = numUnitsDisplay*2; label = 'Principal Component Analysis';
    end

% load('trouble22467_12-14DIV_DoublePeakSortNR.SS','-mat'); plotShift = 12;

%% Find Firing Rates
clear spkRate* tmp*

chanUnits = unique(SS.unit(SS.channel==chanNum));
chanUnits = chanUnits(chanUnits>0);

% NO UNITS ON CHANNEL
if isempty(chanUnits)
    continue;

% UNITS ON CHANNEL
else
    for i = 1:length(chanUnits)
    %     spkRateTotal(i,:) = length(SS.time(SS.unit==chanUnits(i)))/timestamps(2);
        spkRateTotal(i,:) = length(SS.time(SS.unit==chanUnits(i)  & SS.time<timestamps(2) & SS.time>timestamps(1)-3*3600))/(timestamps(2)-timestamps(1)+3600*3);
    end;

    % Sort the units in order of firing rate (highest FR first)
    tmpSort = [chanUnits spkRateTotal];
    tmpSort = sortrows(tmpSort,-2);
    chanUnits = tmpSort(:,1);
    spkRateTotal = tmpSort(:,2);

    for i = 1:length(chanUnits)
    %     spkRatePre(i,:) = length(SS.time(SS.unit==chanUnits(i) & SS.time<timestamps(1)))/timestamps(1);
        spkRatePre(i,:) = length(SS.time(SS.unit==chanUnits(i) & SS.time<timestamps(1) & SS.time>timestamps(1)-3*3600))/(3*3600);
        spkRatePost(i,:) = length(SS.time(SS.unit==chanUnits(i) & SS.time>timestamps(1)))/(timestamps(2)-timestamps(1));
    end
end

%% Plot Avg Waveforms

SS.CalculateAvgWaveform
SS.fs = 25000;
ymax = 35; %uV

for i = 1:numUnitsDisplay
    if length(chanUnits)<i
        break
    end;
    dSpkRate = spkRatePost(i)/spkRatePre(i)*100-100;
    if dSpkRate > 0
        col = [0.1 0.9 0.3];
        sym = '+';
    else
        col = [1 0 0];
        sym = '';
    end
    subplot(3,numUnitsDisplay,i+plotShift);
    confplot(1000/SS.fs.*(1:length(SS.avgwaveform.avg(:,chanUnits(i)))),...
        SS.avgwaveform.avg(:,chanUnits(i)),...
        SS.avgwaveform.std(:,chanUnits(i)),...
        SS.avgwaveform.std(:,chanUnits(i)),...
        'LineWidth',2,'Color','k');
    axis tight
    grid on;
    ylim([-ymax ymax]);
    text(2,ymax,[num2str(spkRatePre(i)),' Hz'],'Color','k','HorizontalAlignment','right','VerticalAlignment','top')
    text(2,ymax-5,[sym,num2str(dSpkRate),'%'],'Color',col,'HorizontalAlignment','right','VerticalAlignment','top')
    if i == 1
        ylabel([label,10,10,'Voltage (uV)'])
    else
        ylabel('Voltage (uV)')
    end
    xlabel('Time (ms)')
end

end

%% Save Figure
export_fig(['trouble22467_12-14DIV_spkSortCompareCustom2_chan' num2str(chanNum)],'-png','-transparent',gcf)

end