Init

README.md

@@ -0,0 +1,36 @@
+# Biomedical Technology Assignment 2018, AUTH
+> Spike sorting of filtered LFP recordings dataset using machine learning
+
+## Table of Contents
+
+- [Clone](#Clone)
+- [Execution](#execution)
+- [Support](#support)
+- [License](#license)
+
+---
+
+## Clone
+
+Clone this repo to your local machine using `https://gitlab.com/Apostolof/biomedicaltechnologyassignment2018.git`
+
+---
+
+## Execution
+
+Scripts were written and tested in Matlab R2016a (v9.0.0.341360). Any newer version of the software should also be able to execute the scripts without problems.
+
+---
+
+## Support
+
+Reach out to us:
+
+- [apostolof's email](mailto:apotwohd@gmail.com "apotwohd@gmail.com")
+- [charaldp's email](mailto:charaldp@ece.auth.gr "charaldp@ece.auth.gr")
+
+---
+
+## License
+
+[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://gitlab.com/Apostolof/biomedicaltechnologyassignment2018/blob/master/LICENSE)

spike_sorting.m

@@ -0,0 +1,163 @@
+%% 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
+clear
+
+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, data}
+end
+
+%% Q.1.3
+figure();
+title('Polynomial curve fitting on median-threshold factor value pairs');
+xlabel('Dataset median');
+ylabel('Threshold factor');
+plot(datasetMedians, datasetFactors, 'o');
+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
+
+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;
+            
+            % 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(firstIndex-34:firstIndex+29);
+                    break;
+                end
+            end
+        end
+        sample = sample + 1;
+    end
+    
+    fprintf('%d spikes found for dataset #%d\n', numberOfSpikes, fileIndex);
+    fprintf('actial 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