authbiomedicaltechnologyassignment2018/spike_sorting.m

%% AUTHOR[1] : Apostolos Fanakis (8261)
%% EMAIL[1]  : apostolof@auth.gr
%% AUTHOR[2] : Charalampos Papadiakos (8302)
%% EMAIL[2]  : charaldp@ece.auth.gr
%% AUTHOR[3] : Hlektra Mitsi ()
%% EMAIL[3]  :
%% $DATE     : 28-December-2018 12:45:00 $
%% $Revision : 1.00 $
%% DEVELOPED : 9.0.0.341360 (R2016a)
%% FILENAME  : spike_sorting.m
%%
%% =================================================================================================
%% S.1
datasetMedians = zeros(8);
datasetFactors = zeros(8);

for fileIndex=1:8
    fprintf('Loading test dataset no. %d\n', fileIndex);
    filename = sprintf('dataset\\Data_Test_%d.mat', fileIndex);
    Dataset = load(filename);
    data = double(Dataset.data);
    
    %% Q.1.1
    figure();
    plot(data(1:10000));
    xlim([0, 10000]);
    title(['First 10000 samples of dataset #' num2str(fileIndex)]);
    xlabel('Sample #');
    ylabel('Trivial Unit'); %TODO: Is this mVolts?
    drawnow;
    
    %% Q.1.2
    dataMedian = median(abs(data)/0.6745);
    datasetMedians(fileIndex) = dataMedian;
    
    thresholdFactorInitValue = 2; % k starting value
    thresholdFactorEndValue = 14; % k ending value
    thresholdFactorStep = 0.01; % k jumping step
    numberOfFactors = length(thresholdFactorInitValue:thresholdFactorStep:thresholdFactorEndValue);
    numberOfSpikesPerFactor = zeros(numberOfFactors);
    
    parfor factorIteration=1:numberOfFactors % runs for each k
        % builds threshold
        thresholdFactor = thresholdFactorInitValue + (factorIteration - 1) * thresholdFactorStep;
        threshold = thresholdFactor * dataMedian;
        
        % calculates number of spikes
        sample = 1;
        while sample <= length(data)
            if data(sample) >= threshold
                % spike found
                numberOfSpikesPerFactor(factorIteration) = numberOfSpikesPerFactor(factorIteration) + 1;
                
                % skips cheking until values are below threshold again
                while sample <= length(data)
                    sample = sample + 1;
                    if (data(sample) <= threshold)
                        break;
                    end
                end
            end
            sample = sample + 1;
        end
    end
    
    figure();
    % trims zeros
    numberOfSpikesTrimmed = numberOfSpikesPerFactor(1:find(numberOfSpikesPerFactor,1,'last'));
    endValue = thresholdFactorInitValue + thresholdFactorStep * (length(numberOfSpikesTrimmed) - 1);
    plot(thresholdFactorInitValue:thresholdFactorStep:endValue, numberOfSpikesTrimmed);
    title(['Number of spikes for different values of k for dataset #' num2str(fileIndex)]);
    xlabel('Threshold factor (k)');
    ylabel('Number of spikes');
    hold on;
    plot([thresholdFactorInitValue endValue], [Dataset.spikeNum, Dataset.spikeNum]);
    xlim([thresholdFactorInitValue endValue]);
    drawnow;
    hold off;
    
    % finds dataset's theshold factor k that produces the closest number of
    % spikes ot the ground truth
    [minValue, closestIndex] = min(abs(numberOfSpikesTrimmed-Dataset.spikeNum));
    datasetFactors(fileIndex) = thresholdFactorInitValue + (closestIndex - 1) * thresholdFactorStep;
    
    clear dataset
    clear data
end
fprintf('\n');

%% Q.1.3
figure();
plot(datasetMedians, datasetFactors, 'o');
title('Polynomial curve fitting on median-threshold factor value pairs');
xlabel('Dataset median');
ylabel('Threshold factor');
hold on;
empiricalRule = polyfit(datasetMedians, datasetFactors, 8);
visualizationX = linspace(0, 0.5, 50);
visualizationY = polyval(empiricalRule, visualizationX);
plot(visualizationX, visualizationY);
hold off

%% =================================================================================================
%% S.2
clearvars = {'closestIndex' 'datasetFactors' 'datasetMedians' 'endValue' 'minValue' 'numberOfFactors' ...
    'numberOfSpikesPerFactor' 'numberOfSpikesTrimmed' 'thresholdFactorEndValue' 'thresholdFactorInitValue' ...
    'thresholdFactorStep' 'visualizationX' 'visualizationY'};
clear(clearvars{:})
clear clearvars

for fileIndex=1:4
    fprintf('Loading evaluation dataset no. %d \n', fileIndex);
    filename = sprintf('dataset\\Data_Eval_E_%d.mat', fileIndex);
    Dataset = load(filename);
    data = double(Dataset.data);
    
    %% Q.2.1 and Q.2.2
    dataMedian = median(abs(data)/0.6745);
    
    factorEstimation = polyval(empiricalRule, dataMedian);
    threshold = factorEstimation * dataMedian;
    numberOfSpikes = 0;
    spikesTimesEst(2500) = 0;
    spikesEst(2500, 64) = 0;
    
    % calculates number of spikes
    spikeStartIndex = 1;
    spikeEndIndex = 1;
    sample = 1;
    while sample <= length(data)
        if data(sample) >= threshold
            % spike found
            numberOfSpikes = numberOfSpikes + 1;
            spikeStartIndex = sample;
            
            % skips cheking until values are below threshold again
            while sample <= length(data)
                sample = sample + 1;
                if (data(sample) <= threshold)
                    spikeEndIndex = sample;
                    [~, minIndex] = min(data(spikeStartIndex:spikeEndIndex));
                    [~, maxIndex] = max(data(spikeStartIndex:spikeEndIndex));
                    firstIndex = min([minIndex maxIndex]);
                    % Q.2.1
                    spikesTimesEst(numberOfSpikes) = firstIndex;
                    % Q.2.2
                    spikesEst(numberOfSpikes, :) = data(spikeStartIndex-1+firstIndex-34:spikeStartIndex-1+firstIndex+29);
                    break;
                end
            end
        end
        sample = sample + 1;
    end
    
    fprintf('%d spikes found for dataset #%d\n', numberOfSpikes, fileIndex);
    fprintf('actual number of spikes = %d\n', length(Dataset.spikeTimes));
    fprintf('diff = %d\n\n', abs(length(Dataset.spikeTimes) - numberOfSpikes));
        
    figure();
    hold on;
    for spike=1:numberOfSpikes
        plot(1:64, spikesEst(spike, :));
    end
    drawnow;
    hold off;
end