diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..16b9076 Binary files /dev/null and b/.DS_Store differ diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/.DS_Store b/HW_1_PNS_EE379K-385V_Neural_Eng/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/HW_1_PNS_EE379K-385V_Neural_Eng/.DS_Store differ diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/2010_Raspopovic_Cuff_Decoding.pdf b/HW_1_PNS_EE379K-385V_Neural_Eng/2010_Raspopovic_Cuff_Decoding.pdf new file mode 100644 index 0000000..0a37418 Binary files /dev/null and b/HW_1_PNS_EE379K-385V_Neural_Eng/2010_Raspopovic_Cuff_Decoding.pdf differ diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/c1_dataVis.m b/HW_1_PNS_EE379K-385V_Neural_Eng/c1_dataVis.m new file mode 100644 index 0000000..604a4aa --- /dev/null +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/c1_dataVis.m @@ -0,0 +1,57 @@ +%% +clear all +close all +clc + +%% Import Data +load('data.mat'); + +%% Example: Plot the raw signal +signal=Flex.signal; +labels=Flex.trigger; +TRIG = gettrigger(labels,0.5); +TRIGend = gettrigger(-labels,-0.5); + +figure('units','normalized','Position',[0.1,0.1,0.7,0.4]) +plot((1:length(signal))./fs,zscore(signal)); +hold on; +plot((1:length(signal))./fs,zscore(labels),'y'); +stem(TRIG./fs,ones(length(TRIG),1)*max(zscore(labels)),'Color','g'); +stem(TRIGend./fs,ones(length(TRIG),1)*max(zscore(labels)),'Color','r'); +grid on; grid minor; +xlim([0,length(signal)./fs]) +xlabel('Time (s)') +ylabel('Amplitude (uV)') +title('Raw VF signal with labels for stimulation periods') + +%% Example: PSD estimates +figure('units','normalized','Position',[0.1,0.1,0.5,0.5]) +[rows_act,cols_act,values_act] = find(labels>0); +[rows_rest1,cols_rest,values_rest] = find(labels==0); +notOfInterest = signal(rows_rest1); +signalOfInterest=signal(rows_act); +h = spectrum.welch; % creates the Welch spectrum estimator +SOIf=psd(h,signalOfInterest,'Fs',fs); % calculates and plot the one sided PSD +plot(SOIf); % Plot the one-sided PSD. +temp =get(gca); +temp.Children(1).Color = 'b'; + +% %% Bandpass Filtering +% +% fc1 = 0; % first cutoff frequency in Hz +% fc2 = 100; % second cutoff frequency in Hz +% +% % normalize the frequencies +% Wp = [fc1 fc2]*2/fs; + +% Build a Butterworth bandpass filter of 4th order +% check the "butter" function in matlab + +% Filter data of both classes with a non-causal filter +% Hint: use "filtfilt" function in MATLAB +% filteredSignal = ; + +%% + + + diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.m b/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.m new file mode 100644 index 0000000..b960386 --- /dev/null +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.m @@ -0,0 +1,35 @@ + +filteredSignal = ; % bandapass filtered signal +label = ; % labels of stimulus locations + +WSize = ; % window size in s +Olap = ; % overlap percentage + +%% Extracting Features over overlapping windows + +WSize = floor(WSize*fs); % length of each data frame, 30ms +nOlap = floor(Olap*WSize); % overlap of successive frames, half of WSize +hop = WSize-nOlap; % amount to advance for next data frame +nx = length(signal); % length of input vector +len = fix((nx - (WSize-hop))/hop); %length of output vector = total frames + +% preallocate outputs for speed +[MAV_feature, VAR_feature, featureLabels] = deal(zeros(1,len)); + +Rise1 = gettrigger(label,0.5); % gets the starting points of stimulations +Fall1 = gettrigger(-label,-0.5); % gets the ending points of stimulations + +for i = 1:len + segment = filteredSignal(((i-1)*hop+1):((i-1)*hop+WSize)); + MAV_feature(i) = ; + VAR_feature(i) = ; + + % re-build the label vector to match it with the feature vector + featureLabels(i) = sum(arrayfun(@(t) ((i-1)*hop+1) >= Rise1(t) && ((i-1)*hop+WSize) <= Fall1(t), 1:length(Rise1))); +end + +%% Plotting the features +% Note: when plotting the features, scale the featureLabels to the max of +% the feature values for proper visualization + + diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/c3_classification.m b/HW_1_PNS_EE379K-385V_Neural_Eng/c3_classification.m new file mode 100644 index 0000000..e12fb04 --- /dev/null +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/c3_classification.m @@ -0,0 +1,107 @@ +close all;clc; +%% +% Inputs: +% -------- +% MAVClass1: the features of the VF case (stimulus and rest features) +% MAVClass2: the features of the Pinch case (stimulus and rest features) +% TriggerClass1: labels for VF features (stimulus or rest label) +% TriggerClass2: labels for Pinch features (stimulus or rest label) + +% Build the datasets +MAV_class1 = MAVClass1(find(TriggerClass1==1)); +MAV_rest1 = MAVClass1(find(TriggerClass1==0)); + +VAR_class1 = VARClass1(find(TriggerClass1==1)); +VAR_rest1 = VARClass1(find(TriggerClass1==0)); + +MAV_class2 = MAVClass2(find(TriggerClass2==1)); +MAV_rest2 = MAVClass2(find(TriggerClass2==0)); + +VAR_class2 = VARClass2(find(TriggerClass2==1)); +VAR_rest2 = VARClass2(find(TriggerClass2==0)); + +% Concantenate the rest classes +MAV_rest = [MAV_rest1 MAV_rest2]; +VAR_rest = [VAR_rest1 VAR_rest2]; + + +%% +% Class1 vs Rest dataset +MAV_Data_Class1vsRest = [MAV_class1 MAV_rest]; +MAV_Labels_Class1vsRest = [ones(1,length(MAV_class1)) 2*ones(1,length(MAV_rest))]; + +VAR_Data_Class1vsRest = [VAR_class1 VAR_rest]; +VAR_Labels_Class1vsRest = MAV_Labels_Class1vsRest; + +% Class2 vs Rest dataset +MAV_Data_Class2vsRest = [MAV_class2 MAV_rest]; +MAV_Labels_Class2vsRest = [ones(1,length(MAV_class2)) 2*ones(1,length(MAV_rest))]; + +VAR_Data_Class2vsRest = [VAR_class2 VAR_rest]; +VAR_Labels_Class2vsRest = MAV_Labels_Class2vsRest; + +% Class1 vs Class2 dataset +MAV_Data_Class1vsClass2 = [MAV_class1 MAV_class2]; +MAV_Labels_Class1vsClass2 = [ones(1,length(MAV_class1)) 2*ones(1,length(MAV_class2))]; + +VAR_Data_Class1vsClass2 = [VAR_class1 VAR_class2]; +VAR_Labels_Class1vsClass2 = MAV_Labels_Class1vsClass2; + +%% +% Both feature datasets +MAVVAR_Data_Class1vsRest = [MAV_Data_Class1vsRest; VAR_Data_Class1vsRest]; +MAVVAR_Labels_Class1vsRest = MAV_Labels_Class1vsRest; + +MAVVAR_Data_Class2vsRest = [MAV_Data_Class2vsRest; VAR_Data_Class2vsRest]; +MAVVAR_Labels_Class2vsRest = MAV_Labels_Class2vsRest; + +MAVVAR_Data_Class1vsClass2 = [MAV_Data_Class1vsClass2; VAR_Data_Class1vsClass2]; +MAVVAR_Labels_Class1vsClass2 = MAV_Labels_Class1vsClass2; + +%% +% Classify all combinations (training set) +k = 10; % for k-fold cross validation +c1 = cvpartition(length(MAV_Labels_Class1vsRest),'KFold',k); +c2 = cvpartition(length(VAR_Labels_Class1vsRest),'KFold',k); +c3 = cvpartition(length(MAVVAR_Labels_Class1vsRest),'KFold',k); +c4 = cvpartition(length(MAV_Labels_Class2vsRest),'KFold',k); +c5 = cvpartition(length(VAR_Labels_Class2vsRest),'KFold',k); +c6 = cvpartition(length(MAVVAR_Labels_Class2vsRest),'KFold',k); +c7 = cvpartition(length(MAV_Labels_Class1vsClass2),'KFold',k); +c8 = cvpartition(length(VAR_Labels_Class1vsClass2),'KFold',k); +c9 = cvpartition(length(MAVVAR_Labels_Class1vsClass2),'KFold',k); + +% Repeat the following for i=1:k, and average performance metrics across all iterations +i=1; +% loop over all k-folds and avergae the performance +% for i=1:k + [TstMAVFC1Rest TstMAVErrC1Rest] = classify(MAV_Data_Class1vsRest(c1.test(i))',MAV_Data_Class1vsRest(c1.training(i))',MAV_Labels_Class1vsRest(c1.training(i))); + [TstCM_MAV_C1rest dum1 TstAcc_MAV_C1rest dum2] = confusion(MAV_Labels_Class1vsRest(c1.test(i)), TstMAVFC1Rest); + + [TstVARFC1Rest TstVARErrC1Rest] = classify(VAR_Data_Class1vsRest(c2.test(i))',VAR_Data_Class1vsRest(c2.training(i))',VAR_Labels_Class1vsRest(c2.training(i))); + [TstCM_VAR_C1rest dum1 TstAcc_VAR_C1rest dum2] = confusion(VAR_Labels_Class1vsRest(c2.test(i)), TstVARFC1Rest); + + [TstMAVVARFC1Rest TstMAVVARErrC1Rest] = classify(MAVVAR_Data_Class1vsRest(:,c3.test(i))',MAVVAR_Data_Class1vsRest(:,c3.training(i))',MAVVAR_Labels_Class1vsRest(c3.training(i))); + [TstCM_MAVVAR_C1rest dum1 TstAcc_MAVVAR_C1rest dum2] = confusion(MAVVAR_Labels_Class1vsRest(c3.test(i)), TstMAVVARFC1Rest); + + % Class2 vs Rest + [TstMAVFC2Rest TstMAVErrC2Rest] = classify(MAV_Data_Class2vsRest(c4.test(i))',MAV_Data_Class2vsRest(c4.training(i))',MAV_Labels_Class2vsRest(c4.training(i))); + [TstCM_MAV_C2rest dum1 TstAcc_MAV_C2rest dum2] = confusion(MAV_Labels_Class2vsRest(c4.test(i)), TstMAVFC2Rest); + + [TstVARFC2Rest TstVARErrC2Rest] = classify(VAR_Data_Class2vsRest(c5.test(i))',VAR_Data_Class2vsRest(c5.training(i))',VAR_Labels_Class2vsRest(c5.training(i))); + [TstCM_VAR_C2rest dum1 TstAcc_VAR_C2rest dum2] = confusion(VAR_Labels_Class2vsRest(c5.test(i)), TstVARFC2Rest); + + [TstMAVVARFC2Rest TstMAVVARErrC2Rest] = classify(MAVVAR_Data_Class2vsRest(:,c6.test(i))',MAVVAR_Data_Class2vsRest(:,c6.training(i))',MAVVAR_Labels_Class2vsRest(c6.training(i))); + [TstCM_MAVVAR_C2rest dum1 TstAcc_MAVVAR_C2rest dum2] = confusion(MAVVAR_Labels_Class2vsRest(c6.test(i)), TstMAVVARFC2Rest); + + % Class1 vs Class2 + [TstMAVFC1C2 TstMAVErrC1C2] = classify(MAV_Data_Class1vsClass2(c7.test(i))',MAV_Data_Class1vsClass2(c7.training(i))',MAV_Labels_Class1vsClass2(c7.training(i))); + [TstCM_MAV_C1C2 dum1 TstAcc_MAV_C1C2 dum2] = confusion(MAV_Labels_Class1vsClass2(c7.test(i)), TstMAVFC1C2); + + [TstVARFC1C2 TstVARErrC1C2] = classify(VAR_Data_Class1vsClass2(c8.test(i))',VAR_Data_Class1vsClass2(c8.training(i))',VAR_Labels_Class1vsClass2(c8.training(i))); + [TstCM_VAR_C1C2 dum1 TstAcc_VAR_C1C2 dum2] = confusion(VAR_Labels_Class1vsClass2(c8.test(i)), TstVARFC1C2); + + [TstMAVVARFC1C2 TstMAVVARErrC1C2] = classify(MAVVAR_Data_Class1vsClass2(:,c9.test(i))',MAVVAR_Data_Class1vsClass2(:,c9.training(i))',MAVVAR_Labels_Class1vsClass2(c9.training(i))); + [TstCM_MAVVAR_C1C2 dum1 TstAcc_MAVVAR_C1C2 dum2] = confusion(MAVVAR_Labels_Class1vsClass2(c9.test(i)), TstMAVVARFC1C2); +% end +%% diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/confusion.m b/HW_1_PNS_EE379K-385V_Neural_Eng/confusion.m new file mode 100644 index 0000000..c08bdc6 --- /dev/null +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/confusion.m @@ -0,0 +1,78 @@ +function [ConfMatClass ConfMatAll Accuracy Error] = confusion(TrueLabels, FoundLabels) + +% function [ConfMatClass ConfMatAll Accuracy Error] = +% eegc3_confusion_matrix(TrueLabels, FoundLabels) +% +% Function to recover a confusion matric of classification problem given +% the groiund truth classes and the classes found by a classifier +% +% Inputs: +% +% TrueLabels: Ground truth data labels +% +% FoundLabels: Labels predicted by a classifier +% +% Outputs: +% +% ConfMatClass: Confusion Matrx ClassNum x ClassNum in percentages per +% class +% +% ConfMatAll: Confusion Matrx ClassNum x ClassNum in overall percentages +% +% Accuracy: Total accuracy % (diagonal of Confusion Matrix) +% +% Error: Total Error % (100 - Accuracy) + +if(~isvector(TrueLabels) || ~isvector(FoundLabels)) + disp('[eegc3_confusion_matrix] Labels must be vectors'); + ConfMat = []; + Accuracy = []; + Error = []; + return; +end + +N1 = length(TrueLabels); +N2 = length(FoundLabels); + +if(N1 ~= N2) + disp('[confusion] True and predicted labels must be the same size'); + ConfMat = []; + Accuracy = []; + Error = []; + return; +else + NSample = N1; +end + +% Find number of classes +Classes = sort(unique(TrueLabels)); +NClass = length(Classes); + +if(~isequal(Classes,[1:NClass]) && ~isequal(Classes,[1:NClass]')) + disp('[confusion] Classes must be a vector [1:N] or [1:N]'); +end + +% Compute Confusion Matrix + +ConfMat = zeros(NClass); +ConfMatAll = zeros(NClass); +for i=1:NSample + ConfMat(TrueLabels(i),FoundLabels(i)) = ... + ConfMat(TrueLabels(i),FoundLabels(i)) + 1; +end + +ConfMatAll = 100*ConfMat/sum(ConfMat(:)); + +ConfMatClass = zeros(NClass); +for i=1:NClass + ConfMatClass(i,:) = 100*ConfMat(i,:)/sum(ConfMat(i,:)); +end + +Accuracy = trace(ConfMatAll); + +Error = 100 - Accuracy; + + + + + diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/data.mat b/HW_1_PNS_EE379K-385V_Neural_Eng/data.mat new file mode 100644 index 0000000..84009ad Binary files /dev/null and b/HW_1_PNS_EE379K-385V_Neural_Eng/data.mat differ diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/gettrigger.m b/HW_1_PNS_EE379K-385V_Neural_Eng/gettrigger.m new file mode 100644 index 0000000..441c70e --- /dev/null +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/gettrigger.m @@ -0,0 +1,60 @@ +function TRIG = gettrigger(s,TH,rfp); +% GETTRIGGER identifies trigger points +% +% TRIG = gettrigger( s [,TH [,rfp]]); % trigger at ascending edge +% TRIG = gettrigger(-s [,TH [,rfp]]); % trigger at decending edge +% +% input : s signal +% TH Threshold; default: (max+min)/2 +% rfp refractory period (default=0) +% output: TRIG TRIGGER time points +% +% see also: TRIGG + +% $Id: gettrigger.m,v 1.4 2008/04/15 14:58:50 schloegl Exp $ +% Copyright (C) 2002-2003,2008 by Alois Schloegl + +% This library is free software; you can redistribute it and/or +% modify it under the terms of the GNU Library General Public +% License as published by the Free Software Foundation; either +% Version 2 of the License, or (at your option) any later version. +% +% This library is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +% Library General Public License for more details. +% +% You should have received a copy of the GNU Library General Public +% License along with this library; if not, write to the +% Free Software Foundation, Inc., 59 Temple Place - Suite 330, +% Boston, MA 02111-1307, USA. + + +if nargin<2, + TH = (max(s)+min(s))/2; +end; +if nargin<3, + rfp = 0; +end; + +TRIG = find(diff(sign(s-TH))>0)+1; +% perform check of trigger points + +if (rfp<=0), return; end; + +% refractory period +k0=1; +k1=2; +while k1