diff --git a/README.md b/README.md new file mode 100644 index 0000000..343e903 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/dataset/Data_Eval_E_1.mat b/dataset/Data_Eval_E_1.mat new file mode 100644 index 0000000..ec791de Binary files /dev/null and b/dataset/Data_Eval_E_1.mat differ diff --git a/dataset/Data_Eval_E_2.mat b/dataset/Data_Eval_E_2.mat new file mode 100644 index 0000000..f5786b1 Binary files /dev/null and b/dataset/Data_Eval_E_2.mat differ diff --git a/dataset/Data_Eval_E_3.mat b/dataset/Data_Eval_E_3.mat new file mode 100644 index 0000000..d113603 Binary files /dev/null and b/dataset/Data_Eval_E_3.mat differ diff --git a/dataset/Data_Eval_E_4.mat b/dataset/Data_Eval_E_4.mat new file mode 100644 index 0000000..4e35024 Binary files /dev/null and b/dataset/Data_Eval_E_4.mat differ diff --git a/dataset/Data_Test_1.mat b/dataset/Data_Test_1.mat new file mode 100644 index 0000000..ff4a3e7 Binary files /dev/null and b/dataset/Data_Test_1.mat differ diff --git a/dataset/Data_Test_2.mat b/dataset/Data_Test_2.mat new file mode 100644 index 0000000..eac9be2 Binary files /dev/null and b/dataset/Data_Test_2.mat differ diff --git a/dataset/Data_Test_3.mat b/dataset/Data_Test_3.mat new file mode 100644 index 0000000..401c2e2 Binary files /dev/null and b/dataset/Data_Test_3.mat differ diff --git a/dataset/Data_Test_4.mat b/dataset/Data_Test_4.mat new file mode 100644 index 0000000..628bfb0 Binary files /dev/null and b/dataset/Data_Test_4.mat differ diff --git a/dataset/Data_Test_5.mat b/dataset/Data_Test_5.mat new file mode 100644 index 0000000..d367f8f Binary files /dev/null and b/dataset/Data_Test_5.mat differ diff --git a/dataset/Data_Test_6.mat b/dataset/Data_Test_6.mat new file mode 100644 index 0000000..c28386c Binary files /dev/null and b/dataset/Data_Test_6.mat differ diff --git a/dataset/Data_Test_7.mat b/dataset/Data_Test_7.mat new file mode 100644 index 0000000..b2fa6ad Binary files /dev/null and b/dataset/Data_Test_7.mat differ diff --git a/dataset/Data_Test_8.mat b/dataset/Data_Test_8.mat new file mode 100644 index 0000000..8bc5ad5 Binary files /dev/null and b/dataset/Data_Test_8.mat differ diff --git a/spike_sorting.m b/spike_sorting.m new file mode 100644 index 0000000..1c62a56 --- /dev/null +++ b/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