diff --git a/.gitignore b/.gitignore index 1a48405..e7fb459 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,13 @@ *~ +libsvm/ +plots/ +log/ data/ +messor/ +src/ +matlabData/ *mex* Block* BCI2000* *.mat +*.zip diff --git a/oldMcode/readKin.m b/oldMcode/readKin.m new file mode 100644 index 0000000..12525cb --- /dev/null +++ b/oldMcode/readKin.m @@ -0,0 +1,22 @@ +maxFile=5; + +% addpath('/home/hohlochj/masterarbeit/usedMcode') +% pathToFile='/nfs/wsi/ti/messor/hohlochj/origData/'; +pathToFile='/home/jph/Uni/masterarbeit/Block1_ReachingMovements/'; + +[subjects,numbers]=namesAndNumbers(strcat(pathToFile,'dirs.txt')); + +tmp=cell([maxFile,1]); + +for i=1:size(subjects,2) + subject=subjects{i}; + for number=numbers{i} + parfor j=1:maxFile + tmp{j}=csvread(strcat(pathToFile,sprintf('kin/%s_B1S00%iR0%i.txt',subject,number,j)),1,0); + tmp{j}=tmp{j}(:,1:4); + end + eval(sprintf('kin%s%i = tmp;',subject,number)); + end +end + +save(strcat(pathToFile,'kinematics.mat')) \ No newline at end of file diff --git a/oldMcode/upsample_kin.m b/oldMcode/upsample_kin.m new file mode 100644 index 0000000..488de56 --- /dev/null +++ b/oldMcode/upsample_kin.m @@ -0,0 +1,59 @@ +function [kin_data, startingPoint]=upsample_kin(kin,bci_signal,bci_sf,bci_total_samples,channelNames) +% function for kin data upsampling to EMG data frequency (2500Hz) +% inputs: +% kin: kin data +% bci_signal: signal recorded in bci (emg+eeg+eog+synchro) +% bci_sf: sampling Frequency of bci signal +% bci_total_samples: length of emg signal (2500Hz) +% channelNames: names of channels recorded in bci (1:32 --> eeg, 33:42--> emg, 43: EOG, 44: synchro, 45: null) +% +% +% outputs: +% kin_data: kin data matrix (bci_total_samplesx13 kin variables) of same length as emg data, upsampled +% startingPoint: starting point is the position point in the EMG data (2500 Hz) where +% %kin data starts being recorded (synchronization step). the +% %staringPoint position in EMG data corresponds to the first data point +% %of kinematic data + + %upsampling --> get it equal to bci2000 data frequency + kin_total_samples=length(kin_filt); % kin number of data points + + %1. find the starting point in the bci2000 data within the first 4 seconds (step in synchronization trigger) + + %look for a step in the beginning of kinematic signal + [~,startingPoint]=max(abs(diff(bci_signal(1:4*bci_sf,(find(ismember(channelNames,'Synchro'))))))); + + %2. create a new vector for kinematic data as long as the vector for EMG + %data + + kin_data=NaN(bci_total_samples,13);%creates a new matrix of same lenght as EMG data + + + %for kin variables 2-13 (not time), from 1 to startingPoint, we asign a + %constant value, the same value that each variable has in the first + %data point + for j=2:13 + kin_data(1:startingPoint,j)=kin(1,j); % from 1 to startingPoint, the kinematic data + end + + %we take the time array + time=kin(:,1);%in ms + + for i=1:length(time) %length in kinematic data + + %assign the value of kin data in each time position in 17 Hz to the time + %position that would have in 2500Hz (the rest of the positions stay with a NaN value) + kin_data((startingPoint+floor(time(i)*bci_sf/1000)),2:13)=kin_filt(i,2:13); + + end + + %interpolate NaN values + + X=[1:bci_total_samples]'; + for i=2:13 + Y=(kin_data(1:bci_total_samples,i)); + ix=~isnan(Y); + kin_data(1:bci_total_samples,i)=interp1(X(ix),Y(ix),X,'CUBIC'); + end + +end diff --git a/text/TODO.txt b/text/TODO.txt index f9d55da..0fce01c 100644 --- a/text/TODO.txt +++ b/text/TODO.txt @@ -1,10 +1,17 @@ -filtfilt +filtfilt (done) -Alle VP auf eciton +Alle VP auf eciton (done) + +pBurgOrder bestimmen + +filter -> filtfilt EMG Matching Bewegung EEG - - durch interpolation (ridge) - + - Passt nicht zusammen - insgesamt ~40s zu wenig kin + + synchro Trigger (ch40) + + upsample_kin (in code) + - durch ridge-Regression (EMG dann EEG als feature) + - als Maß correlation ohne gridsearch vergleichen mit gridsearch diff --git a/text/collection.txt b/text/collection.txt index 77dfff3..cd309d7 100644 --- a/text/collection.txt +++ b/text/collection.txt @@ -12,6 +12,6 @@ Move/Rest EMG : 99.9535% Bewegungen EMG: 77.7352% bzw. 55.8440% wenn 0.5s vorher einbezogen EMG bei 1s Pause vor Bewegungsbeginn 70.4908% - - -Vergleich mit/ohne grid Search +allSubjectsAllDays: +EMG: 56.11%, EEG: 42.2% (pburgOrder 50, EMGwindow 0.2, 1s Pause, c=1 fest, minFreq=0, maxFreq=200) +EMG: 55.71%, EEG: 33.32% (see above, minFrequ=8, maxFreq=30) diff --git a/usedMcode/callAll.m b/usedMcode/callAll.m index b7333d3..3d40917 100644 --- a/usedMcode/callAll.m +++ b/usedMcode/callAll.m @@ -1,20 +1,23 @@ addpath('/home/hohlochj/masterarbeit/usedMcode') -pathToFile='/nfs/wsi/ti/messor/hohlochj/origData/dirs.txt'; -% pathToFile='/home/jph/Uni/masterarbeit/Block1_ReachingMovements/dirs.txt'; +pathToFile='/nfs/wsi/ti/messor/hohlochj/origData/'; +%pathToFile='/home/jph/Uni/masterarbeit/Block1_ReachingMovements/'; -windowEMG=0.2; +maxFile=5; +threshold=7500; %EMG is classified as movement + +windowEMG=0.2; %try 1s windowEEG=1; shift=0.2; pburgOrder=50; -minEEGFreq=8; -maxEEGFreq=30; +minEEGFreq=0; %try 8-30 +maxEEGFreq=200; k=10; maxExpC=0; maxPerClass=250; EEG=true; poolObj=parpool(32); -[subjects,numbers]=structForAllDays(pathToFile); +[subjects,numbers]=namesAndNumbers(pathToFile); j=0; for i=1:size(subjects,2) subject=subjects{i}; @@ -25,17 +28,24 @@ meanAccurancysEMG=zeros([j,1]); meanAccurancysEEG=zeros([j,1]); +maxCEMG=zeros([j,k,1]); +maxCEEG=zeros([j,k,1]); +cmScaledEMG=zeros([j,5,5]); +cmScaledEEG=zeros([j,5,5]); j=1; for i=1:size(subjects,2) subject=subjects{i}; for number=numbers{i} % fprintf('%s%i\n',subject,number); - readEEGEciton(subject,number,windowEMG,windowEEG,shift,5,7500,pburgOrder,minEEGFreq,maxEEGFreq); - meanAccurancysEMG(j)=svmEciton(subject,number,false,k,maxExpC,maxPerClass); - meanAccurancysEEG(j)=svmEciton(subject,number,EEG,k,maxExpC,maxPerClass); + readEEG(pathToFile,subject,number,windowEMG,windowEEG,shift,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq); + [meanAccurancysEMG(j),maxCEMG(j,:),cmScaledEMG(j,:,:)]=svmEciton(pathToFile,subject,number,false,k,maxExpC,maxPerClass); + [meanAccurancysEEG(j),maxCEEG(j,:),cmScaledEEG(j,:,:)]=svmEciton(pathToFile,subject,number,EEG,k,maxExpC,maxPerClass); j=j+1; end end -save(sprintf('/nfs/wsi/ti/messor/hohlochj/matlabData/%s_callAll.mat',datestr(datetime('now')))); + +save(strcat(pathToFile,sprintf('../matlabData/%s_callAll.mat',datestr(datetime('now'))))); delete(poolObj) + + diff --git a/usedMcode/generateTrainingData.m b/usedMcode/generateTrainingData.m index 6cc8601..5d94a48 100644 --- a/usedMcode/generateTrainingData.m +++ b/usedMcode/generateTrainingData.m @@ -1,23 +1,34 @@ -function [trainingDataEEG,trainingDataEMG] = generateTrainingData(signal,windowEMG,windowEEG,shift,params,pburgOrder,minEEGFreq,maxEEGFreq) - frequency=params.SamplingRate.NumericValue; - signalWindowEMG=frequency*windowEMG; - signalWindowEEG=frequency*windowEEG; +function [trainingDataEEG,trainingDataEMG,kinPerSec] = generateTrainingData(signal,kin,windowEMG,windowEEG,shift,params,pburgOrder,minEEGFreq,maxEEGFreq) + bci_sf=params.SamplingRate.NumericValue; + signalWindowEMG=bci_sf*windowEMG; + signalWindowEEG=bci_sf*windowEEG; tempEEG=zeros([32,fix(floor(size(signal,1)/signalWindowEEG-1)*windowEEG/shift+1),maxEEGFreq-minEEGFreq+1]); tempEMG=zeros(fix([size(signal,2)-3-32,fix(size(signal,1)/signalWindowEMG-1)*windowEMG/shift+1])); %Filter around 50Hz and below 2 Hz - [A,B]= butter(2,[48 52]/(frequency/2),'stop'); - [C,D]= butter(2,1/(frequency/2),'high'); - [E,F]= butter(2,[148 152]/(frequency/2),'stop'); + [A,B]= butter(2,[48 52]/(bci_sf/2),'stop'); + [C,D]= butter(2,1/(bci_sf/2),'high'); + [E,F]= butter(2,[148 152]/(bci_sf/2),'stop'); - parfor i=1:32 %filter single channel, w/o HEOG, Synchro and 0s - tempEEG(i,:,:)=shiftingPburg(filtfilt(double(E),double(F),filtfilt(double(C),double(D),filtfilt(double(A),double(B),double(signal(:,i))))),frequency,windowEEG,shift,pburgOrder,minEEGFreq,maxEEGFreq); + parfor i=1:32 %filter single channel, w/o EMG, HEOG, Synchro and 0s + tempEEG(i,:,:)=shiftingPburg(filtfilt(double(E),double(F),filtfilt(double(C),double(D),filtfilt(double(A),double(B),double(signal(:,i))))),bci_sf,windowEEG,shift,pburgOrder,minEEGFreq,maxEEGFreq); %filtfilt stat filter end % TODO: 8-30Hz, pburg bei kurzen Fenstern (Ordnung: Je nach SamplingRate, bei 2500: 32-50, Alpha, Beta peaks sollten sichtbar sein) trainingDataEEG=permute(tempEEG,[2 1 3]); parfor i=33:size(signal,2)-3 - tempEMG(i-32,:,:)=waveformLength(filter(C,D,filter(A,B,signal(:,i))),frequency,windowEMG,shift); + tempEMG(i-32,:,:)=waveformLength(filtfilt(double(E),double(F),filtfilt(double(C),double(D),filtfilt(double(A),double(B),double(signal(:,i))))),bci_sf,windowEMG,shift); end trainingDataEMG=permute(tempEMG,[2 1 3]); + + %shift kin according to synch channel in first minute + [~,startingPoint]=max(abs(diff(signal(1:60*bci_sf,size(signal,2)-1)))); + % add time offset in ms to kin-time + offset=startingPoint/bci_sf*1000; + kin(:,1)=kin(:,1)+offset; +% disp(offset) + %add dummy timestamp at the end to fill kinematics + kin=[kin;kin(end,:)]; + kin(end,1)=size(trainingDataEEG,1)*1000*shift+windowEEG*1000; + kinPerSec=shiftingKin(kin,windowEEG,shift); end diff --git a/usedMcode/namesAndNumbers.m b/usedMcode/namesAndNumbers.m new file mode 100644 index 0000000..cedd10a --- /dev/null +++ b/usedMcode/namesAndNumbers.m @@ -0,0 +1,21 @@ +function [names,numbers]=namesAndNumbers(pathToFile) + fileID=fopen(strcat(pathToFile,'dirs.txt')); + text=textscan(fileID,'%s','Delimiter','\n'); + text=text{1}; + j=0; + k=0; + names=cell(1); + for i=1:size(text) + tmp=char(text(i)); + name=strcat(tmp(1:2)); + if (j~=0 && all(squeeze(names{j}==name))) + k=k+1; + numbers{j}(k)=str2double(tmp(end)); + else + j=j+1; + names{j}=name; + k=1; + numbers{j}(k)=str2double(tmp(end)); + end + end +end \ No newline at end of file diff --git a/usedMcode/readEEG.m b/usedMcode/readEEG.m index f7d1925..7a7eace 100644 --- a/usedMcode/readEEG.m +++ b/usedMcode/readEEG.m @@ -1,16 +1,18 @@ -function readEEG(subject,number,windowEMG,windowEEG,shift,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq) +function readEEG(pathToFile, subject,number,windowEMG,windowEEG,shift,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq) - fprintf('start read %s%i %s\n',subject,number,datestr(datetime('now'))); + %fprintf('start read %s%i %s\n',subject,number,datestr(datetime('now'))); trainingDataEEGcell=cell(maxFile,1); trainingDataEMGcell=cell(maxFile,1); classesCell=cell(maxFile,1); + kin=cell([maxFile,1]); for i=1:maxFile - [sig, stat, params] = load_bcidat(sprintf('/home/jph/Uni/masterarbeit/Block1_ReachingMovements/%s/%s_B100%i/AO_B1S001R0%i.dat',subject,subject,number,i)); - [trainingDataEEGcell{i},trainingDataEMGcell{i}]=generateTrainingData(sig,windowEMG,windowEEG,shift,params,pburgOrder,minEEGFreq,maxEEGFreq); + [sig, stat, params] = load_bcidat(strcat(pathToFile,sprintf('%s/%s_B100%i/%s_B1S00%iR0%i.dat',subject,subject,number,subject,number,i))); + tmpKin=csvread(strcat(pathToFile,sprintf('kin/%s_B1S00%iR0%i.txt',subject,number,i)),1,0); + [trainingDataEEGcell{i},trainingDataEMGcell{i},kin{i}]=generateTrainingData(sig,tmpKin(:,1:4),windowEMG,windowEEG,shift,params,pburgOrder,minEEGFreq,maxEEGFreq); classesCell{i}=stat.StimulusCode; - fprintf('%ith file processed\n',i) + %fprintf('%ith file processed\n',i) end clear sig @@ -18,7 +20,8 @@ trainingDataEEG=cell2mat(trainingDataEEGcell); trainingDataEMG=cell2mat(trainingDataEMGcell); classesMat=cell2mat(classesCell); - clear trainingDataEEGcell trainingDataEMGcell classesCell + kinMat=cell2mat(kin); + %clear trainingDataEEGcell trainingDataEMGcell classesCell kin classificationWithPause=classifyAccordingToEMG(trainingDataEEG, trainingDataEMG,classesMat,shift,params.SamplingRate.NumericValue,threshold); clear classesMat @@ -33,10 +36,11 @@ trainingDataEEG=trainingDataEEG(smoothClassification~=-1,:,:); trainingDataEMG=trainingDataEMG(smoothClassification~=-1,:); classification=smoothClassification(smoothClassification~=-1); + kinematics=kinMat(smoothClassification~=-1,:); clear smoothClassification i - save(sprintf('/nfs/wsi/ti/messor/hohlochj/matlabData/%s%i200msWindowEMG1sWindowEEG200msShift1sPauseFreq0to200.mat',subject,number),trainingDataEEG,trainingDataEMG,classification); + save(strcat(pathToFile,sprintf('../matlabData/%s%i200msWindowEMG1sWindowEEG200msShift1sPauseFreq0to200.mat',subject,number)),'trainingDataEEG','trainingDataEMG','classification','kinematics','-v7.3'); - fprintf('end read %s%i %s\n',subject,number,datestr(datetime('now'))); -end \ No newline at end of file + %fprintf('finished reading %s%i %s\n',subject,number,datestr(datetime('now'))); +end diff --git a/usedMcode/readEEGEciton.m b/usedMcode/readEEGEciton.m deleted file mode 100644 index 69661ee..0000000 --- a/usedMcode/readEEGEciton.m +++ /dev/null @@ -1,42 +0,0 @@ -function readEEGEciton(subject,number,windowEMG,windowEEG,shift,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq) - - %fprintf('start read %s%i %s\n',subject,number,datestr(datetime('now'))); - - trainingDataEEGcell=cell(maxFile,1); - trainingDataEMGcell=cell(maxFile,1); - classesCell=cell(maxFile,1); - - for i=1:maxFile - [sig, stat, params] = load_bcidat(sprintf('/nfs/wsi/ti/messor/hohlochj/origData/%s/%s_B100%i/%s_B1S00%iR0%i.dat',subject,subject,number,subject,number,i)); - [trainingDataEEGcell{i},trainingDataEMGcell{i}]=generateTrainingData(sig,windowEMG,windowEEG,shift,params,pburgOrder,minEEGFreq,maxEEGFreq); - classesCell{i}=stat.StimulusCode; - %fprintf('%ith file processed\n',i) - end - - clear sig - - trainingDataEEG=cell2mat(trainingDataEEGcell); - trainingDataEMG=cell2mat(trainingDataEMGcell); - classesMat=cell2mat(classesCell); - clear trainingDataEEGcell trainingDataEMGcell classesCell - - classificationWithPause=classifyAccordingToEMG(trainingDataEEG, trainingDataEMG,classesMat,shift,params.SamplingRate.NumericValue,threshold); - clear classesMat - - smoothClassification=zeros(size(classificationWithPause)); - for i=1:size(classificationWithPause,1) - smoothClassification(i)=round(mode(classificationWithPause(max(i-2,1):min(i+2,end)))); - end - - clear classificationWithPause - - trainingDataEEG=trainingDataEEG(smoothClassification~=-1,:,:); - trainingDataEMG=trainingDataEMG(smoothClassification~=-1,:); - classification=smoothClassification(smoothClassification~=-1); - - clear smoothClassification i - save(sprintf('/nfs/wsi/ti/messor/hohlochj/matlabData/%s%i200msWindowEMG1sWindowEEG200msShift1sPauseFreq0to200.mat',subject,number),'trainingDataEEG','trainingDataEMG','classification','-v7.3'); - - - %fprintf('end read %s%i %s\n',subject,number,datestr(datetime('now'))); -end diff --git a/usedMcode/shiftingKin.m b/usedMcode/shiftingKin.m new file mode 100644 index 0000000..2b6475e --- /dev/null +++ b/usedMcode/shiftingKin.m @@ -0,0 +1,6 @@ +function [kinPerSec]=shiftingKin(kin, windowEEG, shift) + kinPerSec=zeros(fix((max(kin(:,1))-windowEEG*1000)/(shift*1000)),3); + for j=1:size(kinPerSec,1) + kinPerSec(j,:)=sum(diff(kin(kin(:,1)>(j-1)*shift*1000 & kin(:,1)<=(j-1)*shift*1000+windowEEG*1000,2:4))); + end +end \ No newline at end of file diff --git a/usedMcode/structForAllDays.m b/usedMcode/structForAllDays.m deleted file mode 100644 index d8c5f3a..0000000 --- a/usedMcode/structForAllDays.m +++ /dev/null @@ -1,21 +0,0 @@ -function [names,numbers]=structForAllDays(dirs) - fileID=fopen(dirs); - text=textscan(fileID,'%s','Delimiter','\n'); - text=text{1}; - j=0; - k=0; - names=cell(1); - for i=1:size(text) - tmp=char(text(i)); - name=strcat(tmp(1:2)); - if (j~=0 && all(squeeze(names{j}==name))) - k=k+1; - numbers{j}(k)=str2double(tmp(end)); - else - j=j+1; - names{j}=name; - k=1; - numbers{j}(k)=str2double(tmp(end)); - end - end -end \ No newline at end of file diff --git a/usedMcode/svmEciton.m b/usedMcode/svmEciton.m index c91d139..86c0611 100644 --- a/usedMcode/svmEciton.m +++ b/usedMcode/svmEciton.m @@ -1,5 +1,5 @@ -function [meanAccurancy]= svmEciton(subject,number,EEG,k,maxExpC,maxPerClass) - load(sprintf('/nfs/wsi/ti/messor/hohlochj/matlabData/%s%i200msWindowEMG1sWindowEEG200msShift1sPauseFreq0to200.mat',subject,number)); +function [meanAccurancy,maxC, cmScaled]= svmEciton(pathToFile,subject,number,EEG,k,maxExpC,maxPerClass) + load(strcat(pathToFile,sprintf('../matlabData/%s%i200msWindowEMG1sWindowEEG200msShift1sPauseFreq0to200.mat',subject,number))); % fprintf('%i,%i,%i',size(trainingDataEMG,1),size(trainingDataEEG,1),size(classification,1)) addpath('/nfs/wsi/ti/messor/hohlochj/libsvm/matlab'); @@ -32,7 +32,7 @@ [predictions,accurancy(i,:),~]=svmpredict(leaveClasses,leaveOut(:,:),model); cm=cm+confusionmat(leaveClasses,predictions); %confusion matrix end - meanAccurancy=mean(accurancy(:,1)) + meanAccurancy=mean(accurancy(:,1)); fig=figure; cmScaled=zeros(size(cm)); @@ -41,7 +41,7 @@ end imagesc(cmScaled) colorbar(); - saveas(fig,sprintf('/nfs/wsi/ti/messor/hohlochj/plots/%s%i%icm200ms1sPause.fig',subject,number,EEG),'fig'); - save(sprintf('/nfs/wsi/ti/messor/hohlochj/matlabData/%s%i%i200ms1sPause.mat',subject,number,EEG),'meanAccurancy','maxC','cmScaled','-v7.3'); + saveas(fig,strcat(pathToFile,sprintf('../plots/%s%i%icm200ms1sPause.fig',subject,number,EEG)),'fig'); + save(strcat(pathToFile,sprintf('../matlabData/%s%i%i200ms1sPause.mat',subject,number,EEG)),'meanAccurancy','maxC','cmScaled','-v7.3'); fprintf('%s%i%i finished %s\n',subject,number,EEG,datestr(datetime('now'))) end