diff --git a/TODOthesis.pdf b/TODOthesis.pdf new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/TODOthesis.pdf diff --git a/bashscripts/noOfSyn.bash b/bashscripts/noOfSyn.bash new file mode 100755 index 0000000..462a105 --- /dev/null +++ b/bashscripts/noOfSyn.bash @@ -0,0 +1,4 @@ +#!/bin/bash +cd /nfs/wsi/ti/messor/hohlochj/matlabData + +matlab < /home/hohlochj/masterarbeit/usedMcode/noOfSynergies.m > /nfs/wsi/ti/messor/hohlochj/log/$(date +%y-%m-%d--%T)-synergiesAll.log diff --git a/bashscripts/runOnEciton.bash b/bashscripts/runOnEciton.bash new file mode 100755 index 0000000..bf22f51 --- /dev/null +++ b/bashscripts/runOnEciton.bash @@ -0,0 +1,10 @@ +#!/bin/bash + +cd /nfs/wsi/ti/messor/hohlochj/origData + +#find -mindepth 2 -maxdepth 2 -type d > dirs.txt + +#sed 's,\./,,g' -i dirs.txt +#sed '/^[[:blank:]]*$/d' -i dirs.txt + +matlab < /home/hohlochj/masterarbeit/usedMcode/callAll.m > /nfs/wsi/ti/messor/hohlochj/log/$(date +%y-%m-%d--%T)-MoveRest.log diff --git a/bashscripts/runOnEcitonPos.bash b/bashscripts/runOnEcitonPos.bash new file mode 100755 index 0000000..bbbcc95 --- /dev/null +++ b/bashscripts/runOnEcitonPos.bash @@ -0,0 +1,10 @@ +#!/bin/bash + +cd /nfs/wsi/ti/messor/hohlochj/origData + +#find -mindepth 2 -maxdepth 2 -type d > dirs.txt + +#sed 's,\./,,g' -i dirs.txt +#sed '/^[[:blank:]]*$/d' -i dirs.txt + +matlab < /home/hohlochj/masterarbeit/usedMcode/callAllPos.m > /nfs/wsi/ti/messor/hohlochj/log/$(date +%y-%m-%d--%T)-ViaToPos.log diff --git a/matlabCode/balanceClasses.m b/matlabCode/balanceClasses.m new file mode 100644 index 0000000..59bdcd6 --- /dev/null +++ b/matlabCode/balanceClasses.m @@ -0,0 +1,21 @@ +function [trainClasses,trainData]=balanceClasses(trainClasses,trainData,maxPerClass,noClasses) + if maxPerClass > min(histcounts(trainClasses,noClasses)) + maxPerClass=min(histcounts(trainClasses)); + warning('not enough Data per Class; set maxPerClass to %i', maxPerClass); + end + i=1; + % number of elements per class that has to be deleted + delPerClass=histcounts(trainClasses,noClasses)-maxPerClass; + delIndex=zeros(sum(delPerClass),1); % collects the indices that will be deleted + while sum(delPerClass)>0 + r=fix(size(trainClasses,1)*rand(1)+1); + if delPerClass(trainClasses(r)+1)>0 && all(delIndex~=r) + delIndex(i)=r; + i=i+1; + delPerClass(trainClasses(r)+1)=delPerClass(trainClasses(r)+1)-1; + end + end + trainClasses(delIndex)=[]; + trainData(delIndex,:)=[]; + assert(size(trainClasses,1)>0) +end \ No newline at end of file diff --git a/matlabCode/callAll.m b/matlabCode/callAll.m new file mode 100644 index 0000000..5f35c2d --- /dev/null +++ b/matlabCode/callAll.m @@ -0,0 +1,172 @@ +addpath('/home/hohlochj/masterarbeit/usedMcode') +pathToFile='/nfs/wsi/ti/messor/hohlochj/origData/'; +%pathToFile='/home/jph/Uni/masterarbeit/origData/'; +maxFile=5; +threshold=10000; %EMG is classified as movement +EMGChannels={'AbdPolLo','Biceps','Triceps','FrontDelt','MidDelt','BackDelt'}; +noSynergies=3; + +allSubjects=true; %run all subjects and days or only one random day for one random subject +name='MoveRest'; %suffix for the output, has to be valid name for file (no space, /, ...) + +% window and shift given in seconds +windowEMG=0.2; +windowEEG=1; +shiftEMG=0.05; +shiftEEG=0.2; + +eegOffset=0; %predict actions x*shiftEEG after EEG measurement +pburgOrder=250; +minEEGFreq=2; %[Hz] +maxEEGFreq=49; +pause=false; % is all data 1s befor movement dropped? +noLFsamples=5; + +% Parameters tried for ridge regression +ridgeParams=100; +k=10; % k-fold CV (usually 10) +% Parameters tried for SVM: 2^-x to 2^x in steps of 1 in x +maxExpC=0; +maxPerClass=250; +poolObj=parpool(32); + +% structure data of runs +[subjects,numbers]=namesAndNumbers(pathToFile); +numbersMat=cell2mat(numbers); +subjectsForNumbers=cell(size(numbersMat,2),1); +j=1; +for i=1:size(subjects,2) + subject=subjects{i}; + for number=numbers{i} + subjectsForNumbers{j}=subject; + j=j+1; + end +end +j=j-1; %number of trial-days for all subjects + +% run for all 51 runs +if allSubjects + +% meanAccuracysEMG=zeros([j,1]); +% meanAccuracysEEG=zeros([j,1]); +% meanAccuracysLF=zeros([j,1]); +% maxCEMG=zeros([j,k,1]); +% maxCEEG=zeros([j,k,1]); +% maxCLF=zeros([j,k,1]); +% meanAccuracysEMGmovements=zeros([j,1]); +% meanAccuracysEEGmovements=zeros([j,1]); +% meanAccuracysLFmovements=zeros([j,1]); +% maxCEMGmovements=zeros([j,k,1]); +% maxCEEGmovements=zeros([j,k,1]); +% maxCLFmovements=zeros([j,k,1]); +% cmScaledEMG=zeros([j,5,5]); +% cmScaledEEG=zeros([j,5,5]); +% cmScaledLF=zeros([j,5,5]); +% maxRidgeParamIndexEEGkin=zeros([j,3,k]); +% maxRidgeParamIndexEMGkin=zeros([j,3,k]); +% maxRidgeParamIndexLFkin=zeros([j,3,k]); +% correlationEMGkin=zeros([j,3]); %x,y,angle +% correlationEEGkin=zeros([j,3]); +% correlationLFkin=zeros([j,3]); +% maxRidgeParamIndexEEGautoenc=zeros([j,noSynergies,k]); +% maxRidgeParamIndexEMGautoenc=zeros([j,noSynergies,k]); +% maxRidgeParamIndexLFautoenc=zeros([j,noSynergies,k]); +% correlationEMGautoenc=zeros([j,noSynergies]); +% correlationEEGautoenc=zeros([j,noSynergies]); +% correlationLFautoenc=zeros([j,noSynergies]); +% maxRidgeParamIndexEEGpca=zeros([j,noSynergies,k]); +% maxRidgeParamIndexEMGpca=zeros([j,noSynergies,k]); +% maxRidgeParamIndexLFpca=zeros([j,noSynergies,k]); +% correlationEMGpca=zeros([j,noSynergies]); +% correlationEEGpca=zeros([j,noSynergies]); +% correlationLFpca=zeros([j,noSynergies]); +% maxRidgeParamIndexEEGnnmf=zeros([j,noSynergies,k]); +% maxRidgeParamIndexEMGnnmf=zeros([j,noSynergies,k]); +% maxRidgeParamIndexLFnnmf=zeros([j,noSynergies,k]); +% correlationEMGnnmf=zeros([j,noSynergies]); +% correlationEEGnnmf=zeros([j,noSynergies]); +% correlationLFnnmf=zeros([j,noSynergies]); +% maxRidgeParamIndexEEGemg=zeros([j,size(EMGChannels,2),k]); +% correlationEEGemg=zeros([j,size(EMGChannels,2)]); +% maxRidgeParamIndexAutoencKin=zeros([j,noSynergies,k]); + maxRidgeParamIndexPCAKin=zeros([j,noSynergies,k]); + maxRidgeParamIndexNNMFKin=zeros([j,noSynergies,k]); +% correlationAutoencKin=zeros([j,3]); + correlationPCAKin=zeros([j,3]); + correlationNNMFKin=zeros([j,3]); +% correlationViaAutoenc=zeros([j,3]); +% correlationViaPCA=zeros([j,3]); +% correlationViaNNMF=zeros([j,3]); +% viaCorrelationAutoenc=zeros(j,noSynergies); +% viaCorrelationPCA=zeros(j,noSynergies); +% viaCorrelationNNMF=zeros(j,noSynergies); +% correlationEMGViaAutoenc=zeros([j,size(EMGChannels,2)]); + + parfor j=1:size(numbersMat,2) + number=numbersMat(j); + subject=subjectsForNumbers{j}; + savePath=readAll(pathToFile,subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels,noSynergies); +% [meanAccuracysEMG(j),maxCEMG(j,:),~]=svmEciton(savePath,'EMG',k,maxExpC,maxPerClass,eegOffset,true); +% [meanAccuracysEEG(j),maxCEEG(j,:),~]=svmEciton(savePath,'EEG',k,maxExpC,maxPerClass,eegOffset,true); +% [meanAccuracysLF(j),maxCLF(j,:),~]=svmEciton(savePath,'LF',k,maxExpC,maxPerClass,eegOffset,true); +% [meanAccuracysEMGmovements(j),maxCEMGmovements(j,:),cmScaledEMG(j,:,:)]=svmEciton(savePath,'EMG',k,maxExpC,maxPerClass,eegOffset,false); +% [meanAccuracysEEGmovements(j),maxCEEGmovements(j,:),cmScaledEEG(j,:,:)]=svmEciton(savePath,'EEG',k,maxExpC,maxPerClass,eegOffset,false); +% [meanAccuracysLFmovements(j),maxCLFmovements(j,:),cmScaledLF(j,:,:)]=svmEciton(savePath,'LF',k,maxExpC,maxPerClass,eegOffset,false); +% [correlationEMGkin(j,:),maxRidgeParamIndexEMGkin(j,:,:)]=ridgeCV(savePath,'EMG','kin',k,ridgeParams,eegOffset); +% [correlationEEGkin(j,:),maxRidgeParamIndexEEGkin(j,:,:)]=ridgeCV(savePath,'EEG','kin',k,ridgeParams,eegOffset); +% [correlationLFkin(j,:),maxRidgeParamIndexLFkin(j,:,:)]=ridgeCV(savePath,'LF','kin',k,ridgeParams,eegOffset); +% [correlationEMGautoenc(j,:),maxRidgeParamIndexEMGautoenc(j,:,:)]=ridgeCV(savePath,'EMG','Autoenc',k,ridgeParams,eegOffset); +% [correlationEEGautoenc(j,:),maxRidgeParamIndexEEGautoenc(j,:,:)]=ridgeCV(savePath,'EEG','Autoenc',k,ridgeParams,eegOffset); +% [correlationLFautoenc(j,:),maxRidgeParamIndexLFautoenc(j,:,:)]=ridgeCV(savePath,'LF','Autoenc',k,ridgeParams,eegOffset); +% [correlationEMGpca(j,:),maxRidgeParamIndexEMGpca(j,:,:)]=ridgeCV(savePath,'EMG','PCA',k,ridgeParams,eegOffset); +% [correlationEEGpca(j,:),maxRidgeParamIndexEEGpca(j,:,:)]=ridgeCV(savePath,'EEG','PCA',k,ridgeParams,eegOffset); +% [correlationLFpca(j,:),maxRidgeParamIndexLFpca(j,:,:)]=ridgeCV(savePath,'LF','PCA',k,ridgeParams,eegOffset); +% [correlationEMGnnmf(j,:),maxRidgeParamIndexEMGnnmf(j,:,:)]=ridgeCV(savePath,'EMG','NNMF',k,ridgeParams,eegOffset); +% [correlationEEGnnmf(j,:),maxRidgeParamIndexEEGnnmf(j,:,:)]=ridgeCV(savePath,'EEG','NNMF',k,ridgeParams,eegOffset); +% [correlationLFnnmf(j,:),maxRidgeParamIndexLFnnmf(j,:,:)]=ridgeCV(savePath,'LF','NNMF',k,ridgeParams,eegOffset); +% [correlationEEGemg(j,:),maxRidgeParamIndexEEGemg(j,:,:)]=ridgeCV(savePath,'EEG','EMG',k,ridgeParams,eegOffset); +% [correlationAutoencKin(j,:),maxRidgeParamIndexAutoencKin(j,:,:)]=ridgeCV(savePath,'Autoenc','kin',k,ridgeParams,eegOffset); + [correlationPCAKin(j,:),maxRidgeParamIndexPCAKin(j,:,:)]=ridgeCV(savePath,'PCA','kin',k,ridgeParams,eegOffset); + [correlationNNMFKin(j,:),maxRidgeParamIndexNNMFKin(j,:,:)]=ridgeCV(savePath,'NNMF','kin',k,ridgeParams,eegOffset); +% [correlationViaAutoenc(j,:),viaCorrelationAutoenc(j,:)]=ridgeCVvia(savePath,'EEG','Autoenc','kin',k,ridgeParams,eegOffset); +% [correlationViaPCA(j,:),viaCorrelationPCA(j,:)]=ridgeCVvia(savePath,'EEG','PCA','kin',k,ridgeParams,eegOffset); +% [correlationViaNNMF(j,:),viaCorrelationNNMF(j,:)]=ridgeCVvia(savePath,'EEG','NNMF','kin',k,ridgeParams,eegOffset); +% [correlationEMGViaAutoenc(j,:),~]=ridgeCVvia(savePath,'EEG','Autoenc','EMG',k,ridgeParams,eegOffset); + + fprintf('%s%i finished %s\n',subject,number,datestr(datetime('now'))) + end + + save(strcat(pathToFile,sprintf('../matlabData/%s_callAll-%s.mat',datestr(datetime('now')),name))); +else + % run only for one random run + j=fix(rand()*size(numbersMat,2)+1); + + number=numbersMat(j); + subject=subjectsForNumbers{j}; + savePath=readAll(pathToFile,subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels,noSynergies); +% [meanAccuracysEMG,maxCEMG,cmScaledEMG(:,:)]=svmEciton(savePath,'EMG',k,maxExpC,maxPerClass,eegOffset,false); +% [meanAccuracysEEG,maxCEEG,cmScaledEEG(:,:)]=svmEciton(savePath,'EEG',k,maxExpC,maxPerClass,eegOffset,false); +% [meanAccuracysLF,maxCLF,cmScaledLF(:,:)]=svmEciton(savePath,'LF',k,maxExpC,maxPerClass,eegOffset,false); +% [correlationEMGkin,maxRidgeParamIndexEMGkin(:,:)]=ridgeCV(savePath,'EMG','kin',k,ridgeParams,eegOffset); +% [correlationEEGkin,maxRidgeParamIndexEEGkin(:,:)]=ridgeCV(savePath,'EEG','kin',k,ridgeParams,eegOffset); +% [correlationLFkin,maxRidgeParamIndexLFkin(:,:)]=ridgeCV(savePath,'LF','kin',k,ridgeParams,eegOffset); +% [correlationEMGautoenc,maxRidgeParamIndexEMGautoenc(:,:)]=ridgeCV(savePath,'EMG','Autoenc',k,ridgeParams,eegOffset); +% [correlationEEGautoenc,maxRidgeParamIndexEEGautoenc(:,:)]=ridgeCV(savePath,'EEG','Autoenc',k,ridgeParams,eegOffset); +% [correlationLFautoenc,maxRidgeParamIndexLFautoenc(:,:)]=ridgeCV(savePath,'LF','Autoenc',k,ridgeParams,eegOffset); +% [correlationEMGpca,maxRidgeParamIndexEMGpca(:,:)]=ridgeCV(savePath,'EMG','PCA',k,ridgeParams,eegOffset); +% [correlationEEGpca,maxRidgeParamIndexEEGpca(:,:)]=ridgeCV(savePath,'EEG','PCA',k,ridgeParams,eegOffset); +% [correlationLFpca,maxRidgeParamIndexLFpca(:,:)]=ridgeCV(savePath,'LF','PCA',k,ridgeParams,eegOffset); +% [correlationEMGnnmf,maxRidgeParamIndexEMGnnmf(:,:)]=ridgeCV(savePath,'EMG','NNMF',k,ridgeParams,eegOffset); +% [correlationEEGnnmf,maxRidgeParamIndexEEGnnmf(:,:)]=ridgeCV(savePath,'EEG','NNMF',k,ridgeParams,eegOffset); +% [correlationLFnnmf,maxRidgeParamIndexLFnnmf(:,:)]=ridgeCV(savePath,'LF','NNMF',k,ridgeParams,eegOffset); + [correlationEEGemg,maxRidgeParamIndexEEGemg(:,:)]=ridgeCV(savePath,'EEG','EMG',k,ridgeParams,eegOffset); + [correlationAutoencKin,maxRidgeParamIndexAutoencKin(:,:)]=ridgeCV(savePath,'Autoenc','kin',k,ridgeParams,eegOffset); + [correlationViaAutoenc,viaCorrelationAutoenc]=ridgeCVvia(savePath,'EEG','Autoenc','kin',k,[100],eegOffset); + [correlationViaPCA,viaCorrelationPCA]=ridgeCVvia(savePath,'EEG','PCA','kin',k,[100],eegOffset); + [correlationViaNNMF,viaCorrelationNNMF]=ridgeCVvia(savePath,'EEG','NNMF','kin',k,[100],eegOffset); + [correlationEMGViaAutoenc,~]=ridgeCVvia(savePath,'EEG','Autoenc','EMG',k,ridgeParams,eegOffset); + fprintf('%s%i finished %s\n',subject,number,datestr(datetime('now'))) + + save(strcat(pathToFile,sprintf('../matlabData/%s_call%s%i-%s.mat',datestr(datetime('now')),subject,number,name))); +end +delete(poolObj) diff --git a/matlabCode/callAllPos.m b/matlabCode/callAllPos.m new file mode 100644 index 0000000..4654eed --- /dev/null +++ b/matlabCode/callAllPos.m @@ -0,0 +1,105 @@ +addpath('/home/hohlochj/masterarbeit/usedMcode') +pathToFile='/nfs/wsi/ti/messor/hohlochj/origData/'; +%pathToFile='/home/jph/Uni/masterarbeit/origData/'; +maxFile=1; +threshold=10000; %EMG is classified as movement +EMGChannels={'AbdPolLo','Biceps','Triceps','FrontDelt','MidDelt','BackDelt'}; +noSynergies=2; + +allSubjects=true; %run all subjects and days or only one random day for one random subject +name='Via2Syn'; %suffix for the output, has to be valid name for file (no space, /, ...) + +% window and shift given in seconds +windowEMG=0.2; +windowEEG=1; +shiftEMG=0.05; +shiftEEG=0.2; + +eegOffset=0; %predict actions x*shiftEEG after EEG measurement +pburgOrder=250; +minEEGFreq=2; +maxEEGFreq=49; +pause=false; +noLFsamples=5; +% Parameters tried for ridge regression +ridgeParams=100; +k=10; % k-fold CV (usually 10) +% Parameters tried for SVM: 2^-x to 2^x in steps of 1 in x +maxExpC=0; +maxPerClass=250; +poolObj=parpool(32); + +% structure data of runs +[subjects,numbers]=namesAndNumbers(pathToFile); +numbersMat=cell2mat(numbers); +subjectsForNumbers=cell(size(numbersMat,2),1); +j=1; +for i=1:size(subjects,2) + subject=subjects{i}; + for number=numbers{i} + subjectsForNumbers{j}=subject; + j=j+1; + end +end +j=j-1; %number of trial-days for all subjects + +% run for all 51 runs +if allSubjects +% maxRidgeParamIndexEEGpos=zeros([j,3,k]); +% maxRidgeParamIndexEMGpos=zeros([j,3,k]); +% maxRidgeParamIndexLFpos=zeros([j,3,k]); +% correlationEMGpos=zeros([j,3]); %x,y,angle +% correlationEEGpos=zeros([j,3]); +% correlationLFpos=zeros([j,3]); + +% maxRidgeParamIndexAutoencKin=zeros([j,noSynergies,k]); +% correlationAutoencKin=zeros([j,3]); + correlationViaAutoenc=zeros([j,3]); + correlationViaPCA=zeros([j,3]); + correlationViaNNMF=zeros([j,3]); + viaCorrelationAutoenc=zeros(j,noSynergies); + viaCorrelationPCA=zeros(j,noSynergies); + viaCorrelationNNMF=zeros(j,noSynergies); + + parfor j=1:size(numbersMat,2) + number=numbersMat(j); + subject=subjectsForNumbers{j}; + savePath=readAllPos(pathToFile, subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,... + threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels,noSynergies); +% [correlationEMGpos(j,:),maxRidgeParamIndex(j,:,:)]=ridgeCV(savePath,'EMG','kin',k,ridgeParams,eegOffset); +% [correlationEEGpos(j,:),maxRidgeParamIndex(j,:,:)]=ridgeCV(savePath,'EEG','kin',k,ridgeParams,eegOffset); +% [correlationLFpos(j,:),maxRidgeParamIndex(j,:,:)]=ridgeCV(savePath,'LF','kin',k,ridgeParams,eegOffset); + +% [correlationAutoencKin(j,:),maxRidgeParamIndexAutoencKin(j,:,:)]=ridgeCV(savePath,'Autoenc','kin',k,ridgeParams,eegOffset); + [correlationViaAutoenc(j,:),viaCorrelationAutoenc(j,:)]=ridgeCVvia(savePath,'EEG','Autoenc','kin',k,ridgeParams,eegOffset); + [correlationViaPCA(j,:),viaCorrelationPCA(j,:)]=ridgeCVvia(savePath,'EEG','PCA','kin',k,ridgeParams,eegOffset); + [correlationViaNNMF(j,:),viaCorrelationNNMF(j,:)]=ridgeCVvia(savePath,'EEG','NNMF','kin',k,ridgeParams,eegOffset); + + fprintf('%s%i finished %s\n',subject,number,datestr(datetime('now'))) + end + + save(strcat(pathToFile,sprintf('../matlabData/%s_callAll-%s.mat',datestr(datetime('now')),name))); +else + % run only for one random run + j=fix(rand()*size(numbersMat,2)+1); + + number=numbersMat(j); + subject=subjectsForNumbers{j}; + savePath=readAllPos(pathToFile, subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,... + threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels,noSynergies); + [correlationEMGpos,maxRidgeParamIndex(:,:)]=ridgeCV(savePath,'EMG','kin',k,ridgeParams,eegOffset); + [correlationEEGpos,maxRidgeParamIndex(:,:)]=ridgeCV(savePath,'EEG','kin',k,ridgeParams,eegOffset); + [correlationLFpos,maxRidgeParamIndex(:,:)]=ridgeCV(savePath,'LF','kin',k,ridgeParams,eegOffset); + + [correlationEEGemg,maxRidgeParamIndexEEGemg(:,:)]=ridgeCV(savePath,'EEG','EMG',k,ridgeParams,eegOffset); + [correlationAutoencKin,maxRidgeParamIndexAutoencKin(:,:)]=ridgeCV(savePath,'Autoenc','kin',k,ridgeParams,eegOffset); + [correlationViaAutoenc,viaCorrelationAutoenc]=ridgeCVvia(savePath,'EEG','Autoenc','kin',k,[100],eegOffset); + [correlationViaPCA,viaCorrelationPCA]=ridgeCVvia(savePath,'EEG','PCA','kin',k,[100],eegOffset); + [correlationViaNNMF,viaCorrelationNNMF]=ridgeCVvia(savePath,'EEG','NNMF','kin',k,[100],eegOffset); + [correlationEMGViaAutoenc,~]=ridgeCVvia(savePath,'EEG','Autoenc','EMG',k,ridgeParams,eegOffset); + + save(strcat(pathToFile,sprintf('../matlabData/%s_call%s%i-%s.mat',datestr(datetime('now')),subject,number,name))); +end +delete(poolObj) + + diff --git a/matlabCode/classifyAccordingToEMG.m b/matlabCode/classifyAccordingToEMG.m new file mode 100644 index 0000000..92b261d --- /dev/null +++ b/matlabCode/classifyAccordingToEMG.m @@ -0,0 +1,31 @@ +function [classification]=classifyAccordingToEMG(trainingDataEMG, classes,shift,frequency,threshold,pause) +% classifyAccordingToEMG data classified as -1 should be taken out + sizeTrainingData=size(trainingDataEMG,1); + classification=zeros([sizeTrainingData,1]); + oldclass=0; + for i=1:sizeTrainingData + if mean(trainingDataEMG(i,:)) > threshold % is movement + class=mode(classes(max(round(shift*frequency*i-frequency),1):min(round(i*shift*frequency-1),end))); + if class <= 5 % task is to move + if oldclass==0 + for j=max(fix(i-1/shift),1):fix(i-0.5/shift) % Pause before beginning of movement + classification(j)=-1; %code for pause + end + for j=max(fix(i-0.5/shift),1):i % half a second before executed movement there should be activity (esp. EEG) + if pause + classification(j)=-1; + else + classification(j)=class; + end + end + oldclass=class; + end + classification(i)=class; + else % task is to rest but movement not finished + classification(i)=oldclass; + end + else + oldclass=0; + end + end +end \ No newline at end of file diff --git a/matlabCode/correlation2.m b/matlabCode/correlation2.m new file mode 100644 index 0000000..6405984 --- /dev/null +++ b/matlabCode/correlation2.m @@ -0,0 +1,11 @@ +function [R2]=correlation2(A,B) +% correlation2 calculates the elementwise squared correlation of two +% matrices + if size(A)~=size(B) + error('dimension mismatch') + end + R2=zeros(size(A,2),1); + for i=1:size(A,2) + R2(i)=corr(A(:,i),B(:,i)).^2; + end +end \ No newline at end of file diff --git a/matlabCode/evaluation.m b/matlabCode/evaluation.m new file mode 100644 index 0000000..1558984 --- /dev/null +++ b/matlabCode/evaluation.m @@ -0,0 +1,63 @@ +% Collection of calls to create .mat for evaluation +load('/home/jph/Uni/masterarbeit/evaluation.mat') + +% the data for the corresponding setting has to be loaded + +accuracys.EEG.(sprintf('%s',name))=meanAccuracysEEG; +accuracys.EMG.(sprintf('%s',name))=meanAccuracysEMG; +accuracys.LF.(sprintf('%s',name))=meanAccuracysLF; + +correlations.EEG.Autoenc.(sprintf('%s',name))=correlationEEGautoenc; +correlations.EEG.PCA.(sprintf('%s',name))=correlationEEGpca; +correlations.EEG.NNMF.(sprintf('%s',name))=correlationEEGnnmf; +correlations.EEG.kin.(sprintf('%s',name))=correlationEEGkin; + +correlations.EMG.Autoenc.(sprintf('%s',name))=correlationEMGautoenc; +correlations.EMG.PCA.(sprintf('%s',name))=correlationEMGpca; +correlations.EMG.NNMF.(sprintf('%s',name))=correlationEMGnnmf; +correlations.EMG.kin.(sprintf('%s',name))=correlationEMGkin; + +correlations.LF.Autoenc.(sprintf('%s',name))=correlationLFautoenc; +correlations.LF.PCA.(sprintf('%s',name))=correlationLFpca; +correlations.LF.NNMF.(sprintf('%s',name))=correlationLFnnmf; +correlations.LF.kin.(sprintf('%s',name))=correlationLFkin; + +correlations.EEG.pos.(sprintf('%s',name))=correlationEEGpos; +correlations.EMG.pos.(sprintf('%s',name))=correlationEMGpos; +correlations.LF.pos.(sprintf('%s',name))=correlationLFpos; + +correlations.EEG.emg.Default=correlationEEGemg; +correlations.EEG.kin.ViaAutoenc=correlationViaAutoenc; +correlations.EEG.kin.ViaNNMF=correlationViaNNMF; +correlations.EEG.kin.ViaPCA=correlationViaPCA; + +correlations.Autoenc.kin.Default=correlationAutoencKin; + + +correlations.EEG.pos.ViaAutoenc=correlationViaAutoenc; +correlations.EEG.pos.ViaNNMF=correlationViaNNMF; +correlations.EEG.pos.ViaPCA=correlationViaPCA; +correlations.EEG.emg.ViaAutoenc=correlationEMGViaAutoenc; + +correlations.Autoenc.pos.Default=correlationAutoencKin; + + +synergies.Autoenc.posVia4=correlationViaAutoenc; +synergies.PCA.posVia4=correlationViaPCA; +synergies.NNMF.posVia4=correlationViaNNMF; + +synergies.Autoenc.Default4=viaCorrelationAutoenc; +synergies.PCA.Default4=viaCorrelationPCA; +synergies.NNMF.Default4=viaCorrelationNNMF; + + +synergies.Autoenc.posVia2=correlationViaAutoenc; +synergies.PCA.posVia2=correlationViaPCA; +synergies.NNMF.posVia2=correlationViaNNMF; + +synergies.Autoenc.Default2=viaCorrelationAutoenc; +synergies.PCA.Default2=viaCorrelationPCA; +synergies.NNMF.Default2=viaCorrelationNNMF; + +save('/home/jph/Uni/masterarbeit/evaluation.mat','accuracys','correlations','synergies') +clear diff --git a/matlabCode/evaluationAccuracys.m b/matlabCode/evaluationAccuracys.m new file mode 100644 index 0000000..f460e86 --- /dev/null +++ b/matlabCode/evaluationAccuracys.m @@ -0,0 +1,48 @@ +% Collection of calls to evaluate the accuracys +load('/home/jph/Uni/masterarbeit/evaluation.mat') + +figureSavePath='/home/jph/Uni/masterarbeit/text/thesis/pictures/results/'; +% mySaveFigure(gcf,strcat(figureSavePath,'plot')) +%% compare forms of recording +eegAcc=struct2array(accuracys.EEG); +emgAcc=struct2array(accuracys.EMG); +LFAcc=struct2array(accuracys.LF); +EEG=cell(size(eegAcc(:,1:end-1),2),1); +EEG(:)={'EEG'}; +EMG=cell(size(emgAcc(:,1:end-1),2),1); +EMG(:)={'EMG'}; +LF=cell(size(LFAcc(:,1:end-1),2),1); +LF(:)={'LF'}; + +anova1(cat(2,eegAcc(:,1:end-1),emgAcc(:,1:end-1),LFAcc(:,1:end-1)),cat(1,EEG,EMG,LF)) +ylabel('% classified correctly') +title('ANOVA for EEG, EMG and LF') + +%% evaluate single form of recording +input=accuracys.EMG; +sizeY=2; +limits_y=[100,0]; + +names=fieldnames(input); +noOfPlots=size(names,1); + +for i=1:noOfPlots + limits_y=[min(min(input.(sprintf('%s',names{i}))),limits_y(1)),... + max(max(input.(sprintf('%s',names{i}))),limits_y(2))]; +end +limits_y=[limits_y(1)-0.1*diff(limits_y),limits_y(2)+0.1*diff(limits_y)] + +figure(); +sizeX=ceil(noOfPlots/sizeY); +for i=1:noOfPlots + subplot(sizeY,sizeX,i) + boxplot(input.(sprintf('%s',names{i}))) + title(names{i}) + ylim(limits_y) + if ceil(i/sizeX) > ceil((i-1)/sizeX) + ylabel('% classified correctly') + end +end + +anova1(cat(2,input.default3Syn,input.offset1Syn3,input.offset2Syn3,input.pause1Syn3,input.pause1Off1Syn3),[0,0,0,1,1]) +anova1(cat(2,input.default3Syn,input.offset1Syn3,input.pause1Syn3,input.pause1Off1Syn3),[0,1,0,1]) diff --git a/matlabCode/evaluationCorrelations.m b/matlabCode/evaluationCorrelations.m new file mode 100644 index 0000000..0d0bb81 --- /dev/null +++ b/matlabCode/evaluationCorrelations.m @@ -0,0 +1,180 @@ +% Collection of calls to evaluate the correlations +load('/home/jph/Uni/masterarbeit/evaluation.mat') + +figureSavePath='/home/jph/Uni/masterarbeit/text/thesis/pictures/results/'; +% mySaveFigure(gcf,strcat(figureSavePath,'plot')) +%% inspect methods of recoding +input=correlations.EEG.pos.Default; +anova1(input,{'x','y','theta'}) +anova1(input(:,1:2),{'x','y'}) +anova1([input(:,1),input(:,3)],{'x','theta'}) + +%% compare methods of recording +% velocities +eegCorrKin=pickFromStruct(correlations.EEG.kin,1:3); +emgCorrKin=correlations.EMG.kin.Default; +LFCorrKin=pickFromStruct(correlations.LF.kin,1:3); +EEG=cell(3*size(eegCorrKin,1),1); +EEG(:)={'EEG'}; +EMG=cell(3,1); +EMG(:)={'EMG'}; +LF=cell(3*size(LFCorrKin,1),1); +LF(:)={'LF'}; + +anova1(cell2mat(cat(1,eegCorrKin,emgCorrKin,LFCorrKin)'),cat(1,EEG,EMG,LF)) + +ylabel('correlation recorded - predicted') +title('prediction of velocities for EEG, EMG and LF') +anova1(cell2mat(cat(1,eegCorrKin,emgCorrKin)'),cat(1,EEG,EMG)) +anova1(cell2mat(cat(1,eegCorrKin,LFCorrKin)'),cat(1,EEG,LF)) +anova1(cell2mat(cat(1,emgCorrKin,LFCorrKin)'),cat(1,EMG,LF)) + +% positions +eegCorrPos=pickFromStruct(correlations.EEG.pos,1:3); +emgCorrPos=correlations.EMG.pos.Default; +LFCorrPos=pickFromStruct(correlations.LF.pos,1:3); +EEG=cell(3*size(eegCorrPos,1),1); +EEG(:)={'EEG'}; +EMG=cell(3,1); +EMG(:)={'EMG'}; +LF=cell(3*size(LFCorrPos,1),1); +LF(:)={'LF'}; + +anova1(cell2mat(cat(1,eegCorrPos,emgCorrPos,LFCorrPos)'),cat(1,EEG,EMG,LF)) +ylabel('correlation recorded - predicted') +title('prediction of positions for EEG, EMG and LF') +anova1(cell2mat(cat(1,eegCorrPos,emgCorrPos)'),cat(1,EEG,EMG)) +anova1(cell2mat(cat(1,eegCorrPos,LFCorrPos)'),cat(1,EEG,LF)) +anova1(cell2mat(cat(1,emgCorrPos,LFCorrPos)'),cat(1,EMG,LF)) + +%% compare direct and via +%kin +direct=correlations.EEG.kin.Default; +viaAutoenc=correlations.EEG.kin.ViaAutoenc; +viaPCA=correlations.EEG.kin.ViaPCA; +viaNNMF=correlations.EEG.kin.ViaNNMF; + +anova1([direct(:),viaAutoenc(:),viaPCA(:),viaNNMF(:)],{'direct','Autoenc','PCA','NMF'}) +ylabel('correlation recorded - predicted') +title('Compare velocity prediction direct and via') + +anova1([viaAutoenc(:),viaPCA(:),viaNNMF(:)],{'Autoenc','PCA','NMF'}) + +%pos +direct=correlations.EEG.pos.Default; +viaAutoenc=correlations.EEG.pos.ViaAutoenc; +viaPCA=correlations.EEG.pos.ViaPCA; +viaNNMF=correlations.EEG.pos.ViaNNMF; + +anova1([direct(:),viaAutoenc(:),viaPCA(:),viaNNMF(:)],{'direct','Autoenc','PCA','NMF'}) +ylabel('correlation recorded - predicted') +title('Compare position prediction direct and via') + +anova1([viaAutoenc(:),viaPCA(:),viaNNMF(:)],{'Autoenc','PCA','NMF'}) + +% EMG +direct=correlations.EEG.emg.Default; +viaAutoenc=correlations.EEG.emg.ViaAutoenc; + +anova1([direct(:),viaAutoenc(:)],{'direct','Autoenc'}) +ylabel('correlation recorded - predicted') +title('Compare EMG prediction direct and viaAutoenc') +mySaveFigure(gcf,'predictEMGfromEEG') + +%% evaluate single combination +input=correlations.LF.Autoenc; +sizeY=2; +limits_y=[1,-1]; + +names=fieldnames(input); +noOfPlots=size(names,1); + +for i=1:noOfPlots + limits_y=[min(min(min(input.(sprintf('%s',names{i})))),limits_y(1)),... + max(max(max(input.(sprintf('%s',names{i})))),limits_y(2))]; +end +limits_y=[limits_y(1)-0.1*diff(limits_y),limits_y(2)+0.1*diff(limits_y)] + +figure(); + +sizeX=ceil(noOfPlots/sizeY); +for i=1:noOfPlots + subplot(sizeY,sizeX,i) + boxplot(input.(sprintf('%s',names{i}))) + title(names{i}) + ylim(limits_y) + if ceil(i/sizeX) > ceil((i-1)/sizeX) + ylabel('correlation recorded - predicted') + end +end + +anova1(cat(2,input.Default(:),input.Offset1(:),input.Offset2(:)),[0,1,2]) +ylabel('Offset') +anova1(cat(2,input.Default(:),input.Offset1(:),input.Pause1(:),input.Pause1Offset1(:)),{'0.5s Pause','0.5s Pause', '1s Pause', '1s Pause'}) + +%% EEG to EMG +EMGChannels={'AbdPolLo','Biceps','Triceps','FrontDelt','MidDelt','BackDelt'}; + +anova1(correlations.EEG.emg.Default,EMGChannels) + +title('EMG from EEG') +mySaveFigure(gcf,strcat(figureSavePath,'EEGemg')) + +%% compare EMG with synergies +anova1([correlations.EMG.pos.Default,correlations.Autoenc.pos.Default],{'EMG','EMG','EMG','Autoenc','Autoenc','Autoenc'}) + +title('Predict Positions from EMG or Synergies') +mySaveFigure(gcf,strcat(figureSavePath,'EMGautoencPos')) + + +anova1([correlations.EMG.kin.Default,correlations.Autoenc.kin.Default],{'EMG','EMG','EMG','Autoenc','Autoenc','Autoenc'}) + +title('Predict Velocities from EMG or Synergies') +mySaveFigure(gcf,strcat(figureSavePath,'EMGautoenc')) + +% prediction from EEG +input=correlations.EEG; + +emgPrediction=input.emg.Default; +autoencPrediction=input.Autoenc.Default; +pcaPrediction=input.PCA.Default; +nnmfPrediction=input.NNMF.Default; + +EMG=cell(size(emgPrediction,2),1); +EMG(:)={'EMG'}; +Autoenc=cell(size(autoencPrediction,2),1); +Autoenc(:)={'Autoenc'}; +PCA=cell(size(pcaPrediction,2),1); +PCA(:)={'PCA'}; +NNMF=cell(size(nnmfPrediction,2),1); +NNMF(:)={'NMF'}; + +anova1([emgPrediction,autoencPrediction,pcaPrediction,nnmfPrediction],cat(1,EMG,Autoenc,PCA,NNMF)) + +title('Predict Synergies or EMG from EEG') +ylabel('correlation') +mySaveFigure(gcf,strcat(figureSavePath,'predictEMGSyn')) + + +anova1([emgPrediction,autoencPrediction],cat(1,EMG,Autoenc)) + +title('Predict Autoencoder or EMG from EEG') +ylabel('correlation') +mySaveFigure(gcf,strcat(figureSavePath,'predictEMGAutoenc')) + +%% compare predictions via different synergies +input=correlations.EEG.kin; + +viaAutoenc=input.ViaAutoenc; +viaPCA=input.ViaPCA; +viaNNMF=input.ViaNNMF; + + +Autoenc=cell(size(viaAutoenc,2),1); +Autoenc(:)={'Autoenc'}; +PCA=cell(size(viaPCA,2),1); +PCA(:)={'PCA'}; +NNMF=cell(size(viaNNMF,2),1); +NNMF(:)={'NMF'}; + +anova1([viaAutoenc,viaPCA,viaNNMF],cat(1,Autoenc,PCA,NNMF)) \ No newline at end of file diff --git a/matlabCode/evaluationSynergies.m b/matlabCode/evaluationSynergies.m new file mode 100644 index 0000000..9102e2d --- /dev/null +++ b/matlabCode/evaluationSynergies.m @@ -0,0 +1,72 @@ +% Collection of calls to evaluate the synergies +load('/home/jph/Uni/masterarbeit/evaluation.mat') + +synergies.Autoenc.Default3=correlations.EEG.Autoenc.Default; +synergies.PCA.Default3=correlations.EEG.PCA.Default; +synergies.NNMF.Default3=correlations.EEG.NNMF.Default; + +synergies.Autoenc.posVia3=correlations.EEG.pos.ViaAutoenc; +synergies.PCA.posVia3=correlations.EEG.pos.ViaPCA; +synergies.NNMF.posVia3=correlations.EEG.pos.ViaNNMF; + +clear correlations accuracys + +%% predict Synergies from EEG +autoencData=[synergies.Autoenc.Default3,synergies.Autoenc.Default4]; +pcaData=[synergies.PCA.Default3,synergies.PCA.Default4]; +nnmfData=[synergies.NNMF.Default3,synergies.NNMF.Default4]; + +% autoencData=autoencData(4:end); +% pcaData=pcaData(4:end); +% nnmfData=nnmfData(4:end); + +autoencGroup=cell(size(autoencData,2),1); +autoencGroup(:)={'Autoenc'}; +pcaGroup=cell(size(pcaData,2),1); +pcaGroup(:)={'PCA'}; +nnmfGroup=cell(size(nnmfData,2),1); +nnmfGroup(:)={'NMF'}; +groups=cat(1,autoencGroup,pcaGroup,nnmfGroup); + +% noSyn=[2,2,3,3,3,4,4,4,4,2,2,3,3,3,4,4,4,4,2,2,3,3,3,4,4,4,4]; + +anova1([autoencData,pcaData,nnmfData],groups) + +%% Predict from predicted Synergies +autoencData=[synergies.Autoenc.posVia2,synergies.Autoenc.posVia3,synergies.Autoenc.posVia4]; +pcaData=[synergies.PCA.posVia2,synergies.PCA.posVia3,synergies.PCA.posVia4]; +nnmfData=[synergies.NNMF.posVia2,synergies.NNMF.posVia3,synergies.NNMF.posVia4]; + +autoencGroup=cell(size(autoencData,2),1); +autoencGroup(:)={'Autoenc'}; +pcaGroup=cell(size(pcaData,2),1); +pcaGroup(:)={'PCA'}; +nnmfGroup=cell(size(nnmfData,2),1); +nnmfGroup(:)={'NMF'}; +groups=cat(1,autoencGroup,pcaGroup,nnmfGroup); + +noSyn=[2,2,2,3,3,3,4,4,4,2,2,2,3,3,3,4,4,4,2,2,2,3,3,3,4,4,4]; + +anova1([autoencData,pcaData,nnmfData],groups) +anova1([autoencData,pcaData,nnmfData],noSyn) + + +%compare 3 and 4 synergies + +autoencData=[synergies.Autoenc.posVia3,synergies.Autoenc.posVia4]; +pcaData=[synergies.PCA.posVia3,synergies.PCA.posVia4]; +nnmfData=[synergies.NNMF.posVia3,synergies.NNMF.posVia4]; + +% autoencData=autoencData(4:end); +% pcaData=pcaData(4:end); +% nnmfData=nnmfData(4:end); + +autoencGroup=cell(size(autoencData,2),1); +autoencGroup(:)={'Autoenc'}; +pcaGroup=cell(size(pcaData,2),1); +pcaGroup(:)={'PCA'}; +nnmfGroup=cell(size(nnmfData,2),1); +nnmfGroup(:)={'NMF'}; +groups=cat(1,autoencGroup,pcaGroup,nnmfGroup); + +noSyn=[3,3,3,4,4,4,3,3,3,4,4,4,3,3,3,4,4,4]; \ No newline at end of file diff --git a/matlabCode/generateTrainingData.m b/matlabCode/generateTrainingData.m new file mode 100644 index 0000000..d959c3a --- /dev/null +++ b/matlabCode/generateTrainingData.m @@ -0,0 +1,38 @@ +function [trainingDataEEG,trainingDataEEGlf,trainingDataEMG,kinPerSec] = generateTrainingData(signal,kin,windowEMG,windowEEG,shiftEMG,shiftEEG,params,pburgOrder,minEEGFreq,maxEEGFreq,noLFsamples,EMGChannels) +% generateTrainingData filters and transforms the signal and creates +% according velocities + bci_sf=params.SamplingRate.NumericValue; + signalWindowEEG=bci_sf*windowEEG; + tempEEG=zeros([32,fix(floor(size(signal,1)/signalWindowEEG-1)*windowEEG/shiftEEG+1),201]); + tempEEGlf=zeros([32,fix(floor(size(signal,1)/signalWindowEEG-1)*windowEEG/shiftEEG+1),noLFsamples]); + + %Filter around 50 Hz and 150 Hz + % filter for the range minEEGFreq to maxEEGFreq + [A,B]= butter(2,[48 52]/(bci_sf/2),'stop'); + [C,D]= butter(2,minEEGFreq/(bci_sf/2),'high'); + [E,F]= butter(2,[148 152]/(bci_sf/2),'stop'); + [G,H]= butter(2,maxEEGFreq/(bci_sf/2),'low'); + + % filter for low Frequencies + [V,W]= butter(2,0.01/(bci_sf/2),'high'); + [X,Y]= butter(2,1/(bci_sf/2),'low'); + + parfor i=1:32 %filter single channel, w/o EMG, HEOG, Synchro and 0s + tempEEG(i,:,:)=shiftingPburg(filtfilt(double(G),double(H),filtfilt(double(E),double(F),filtfilt(double(C),double(D),filtfilt(double(A),double(B),double(signal(:,i)))))),bci_sf,windowEEG,shiftEEG,pburgOrder,0,200); + tempEEGlf(i,:,:)=shiftingDownsample(filtfilt(double(X),double(Y),filtfilt(double(V),double(W),double(signal(:,i)))),bci_sf,windowEEG,shiftEEG,noLFsamples); + end + trainingDataEEG=permute(tempEEG,[2 1 3]); + trainingDataEEGlf=permute(tempEEGlf,[2 1 3]); + + trainingDataEMG=waveformLengthAll(signal(:,ismember(params.ChannelNames.Value,EMGChannels)),bci_sf,windowEMG,shiftEMG); + + %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; + %add dummy timestamp at the end to fill kinematics + kin=[kin;kin(end,:)]; + kin(end,1)=size(trainingDataEMG,1)*1000*shiftEMG+windowEMG*1000; + kinPerSec=shiftingKin(kin,windowEMG,shiftEMG); +end diff --git a/matlabCode/generateTrainingDataPos.m b/matlabCode/generateTrainingDataPos.m new file mode 100644 index 0000000..1c8249b --- /dev/null +++ b/matlabCode/generateTrainingDataPos.m @@ -0,0 +1,38 @@ +function [trainingDataEEG,trainingDataEEGlf,trainingDataEMG,kinPerSec] = generateTrainingDataPos(signal,kin,windowEMG,windowEEG,shiftEMG,shiftEEG,params,pburgOrder,minEEGFreq,maxEEGFreq,noLFsamples,EMGChannels) +% generateTrainingData filters and transforms the signal and creates +% according positions + bci_sf=params.SamplingRate.NumericValue; + signalWindowEEG=bci_sf*windowEEG; + tempEEG=zeros([32,fix(floor(size(signal,1)/signalWindowEEG-1)*windowEEG/shiftEEG+1),201]); + tempEEGlf=zeros([32,fix(floor(size(signal,1)/signalWindowEEG-1)*windowEEG/shiftEEG+1),noLFsamples]); + + %Filter around 50Hz and below 2 Hz + [A,B]= butter(2,[48 52]/(bci_sf/2),'stop'); + [C,D]= butter(2,minEEGFreq/(bci_sf/2),'high'); + [E,F]= butter(2,[148 152]/(bci_sf/2),'stop'); + [G,H]= butter(2,maxEEGFreq/(bci_sf/2),'low'); + + % filter for low Frequencies + [V,W]= butter(2,0.01/(bci_sf/2),'high'); + [X,Y]= butter(2,1/(bci_sf/2),'low'); + + parfor i=1:32 %filter single channel, w/o EMG, HEOG, Synchro and 0s + tempEEG(i,:,:)=shiftingPburg(filtfilt(double(G),double(H),filtfilt(double(E),double(F),filtfilt(double(C),double(D),filtfilt(double(A),double(B),double(signal(:,i)))))),bci_sf,windowEEG,shiftEEG,pburgOrder,0,200); + tempEEGlf(i,:,:)=shiftingDownsample(filtfilt(double(X),double(Y),filtfilt(double(V),double(W),double(signal(:,i)))),bci_sf,windowEEG,shiftEEG,noLFsamples); + end + trainingDataEEG=permute(tempEEG,[2 1 3]); + trainingDataEEGlf=permute(tempEEGlf,[2 1 3]); + + trainingDataEMG=waveformLengthAll(signal(:,ismember(params.ChannelNames.Value,EMGChannels)),bci_sf,windowEMG,shiftEMG); + + %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(trainingDataEMG,1)*1000*shiftEMG+windowEMG*1000; + kinPerSec=shiftingPos(kin,windowEMG,shiftEMG); +end diff --git a/matlabCode/kFoldRidge.m b/matlabCode/kFoldRidge.m new file mode 100644 index 0000000..8a15767 --- /dev/null +++ b/matlabCode/kFoldRidge.m @@ -0,0 +1,31 @@ +function [coeffs,maxRidgeParamIndex]=kFoldRidge(data,kin,k,ridgeParams) +% kFoldRidge does a k-fold Crossvalidation of parameter lambda for Ridge +% Regression + if any(size(ridgeParams)>1) + + correlationPerParam=zeros(size(ridgeParams)); + + for param=ridgeParams + + randomMapping=transpose(randperm(size(data,1))); + correlations=zeros(k,1); + + parfor i=1:k + %split into training and test set + trainData=data(mod(randomMapping,k)+1~=i,:); + testData=data(mod(randomMapping,k)+1==i,:); + trainKin=kin(mod(randomMapping,k)+1~=i); + testKin=kin(mod(randomMapping,k)+1==i); + + coeffs=ridge(trainKin,trainData,param,0); + correlations(i)=ridgeCorrelation(testData,testKin,coeffs); + end + % use mean correlation as correlation for the current parameter + correlationPerParam(ridgeParams==param)=mean(correlations); + end + [~,maxRidgeParamIndex]=max(correlationPerParam); + else + maxRidgeParamIndex=1; + end + coeffs=ridge(kin,data,ridgeParams(maxRidgeParamIndex),0); +end \ No newline at end of file diff --git a/matlabCode/kfoldCV.m b/matlabCode/kfoldCV.m new file mode 100644 index 0000000..5fc6ea7 --- /dev/null +++ b/matlabCode/kfoldCV.m @@ -0,0 +1,36 @@ +function [model,maxC] = kfoldCV(classification,trainingData,k,cExpMax,maxDataPerClass) +% kfoldCV does a k-fold Crossvalidation of parameter c for SVM + noClasses=size(unique(classification),1); + if cExpMax>0 %CV only if necessary + cvAccuracy=zeros(2*cExpMax+1,1); + for cExp=1:2*cExpMax+1 % try different values for c + c=2^(cExp-cExpMax-1); + + randomMapping=transpose(randperm(size(trainingData,1))); + accuracy=zeros(k,3); + + parfor i=1:k + % split in training and test set + trainData=trainingData(mod(randomMapping,k)+1~=i,:,:); + testData=trainingData(mod(randomMapping,k)+1==i,:,:); + trainClasses=classification(mod(randomMapping,k)+1~=i); + testClasses=classification(mod(randomMapping,k)+1==i); + %fprintf('i=%i, k=%i, c=%i\n',i,k,c) + [trainClasses,trainData]=balanceClasses(trainClasses,trainData,... + maxDataPerClass,noClasses); + model=svmtrain(trainClasses,trainData(:,:),sprintf('-t 0 -c %f -q',c)); + [~, accuracy(i,:), ~]=svmpredict(testClasses, testData(:,:), model,'-q'); + end + %for each value for c take mean accuracy as accuracy + cvAccuracy(cExp)=mean(accuracy(:,1)); + end + [~,maxC]=max(cvAccuracy); + else + maxC=1; %no gridsearch since only one C possible + end + bestC=2^(maxC-cExpMax-1); + [balancedClasses, balanceData]=balanceClasses(classification,trainingData,... + maxDataPerClass,noClasses); + model=svmtrain(balancedClasses, balanceData(:,:),sprintf('-t 0 -c %f -q',bestC)); + +end diff --git a/matlabCode/load_bcidat.mexa64 b/matlabCode/load_bcidat.mexa64 new file mode 100644 index 0000000..dd30372 --- /dev/null +++ b/matlabCode/load_bcidat.mexa64 Binary files differ diff --git a/matlabCode/load_bcidat.mexw32 b/matlabCode/load_bcidat.mexw32 new file mode 100644 index 0000000..76ec378 --- /dev/null +++ b/matlabCode/load_bcidat.mexw32 Binary files differ diff --git a/matlabCode/mem.mexa64 b/matlabCode/mem.mexa64 new file mode 100644 index 0000000..12e8af0 --- /dev/null +++ b/matlabCode/mem.mexa64 Binary files differ diff --git a/matlabCode/mem.mexglx b/matlabCode/mem.mexglx new file mode 100644 index 0000000..f314b22 --- /dev/null +++ b/matlabCode/mem.mexglx Binary files differ diff --git a/matlabCode/mem.mexmaci b/matlabCode/mem.mexmaci new file mode 100644 index 0000000..dd58428 --- /dev/null +++ b/matlabCode/mem.mexmaci Binary files differ diff --git a/matlabCode/mem.mexw64 b/matlabCode/mem.mexw64 new file mode 100644 index 0000000..bba28b3 --- /dev/null +++ b/matlabCode/mem.mexw64 Binary files differ diff --git a/matlabCode/myDownsample.m b/matlabCode/myDownsample.m new file mode 100644 index 0000000..14141e7 --- /dev/null +++ b/matlabCode/myDownsample.m @@ -0,0 +1,5 @@ +function [out]=myDownsample(in,noSamples) +% myDownsample samples as noSamples equidistant values + stepSize=max(size(in))/noSamples; + out=in(fix((0:noSamples-1)*stepSize)+1); +end \ No newline at end of file diff --git a/matlabCode/mySaveFigure.m b/matlabCode/mySaveFigure.m new file mode 100644 index 0000000..c52c2aa --- /dev/null +++ b/matlabCode/mySaveFigure.m @@ -0,0 +1,9 @@ +function []=mySaveFigure(fig,filename) +% mySaveFigure saves fig (e.g. gcf) as filename + set(fig,'Position',[251 49 1105 636]); + axes = findobj(fig,'type','axes'); + set(axes,'FontSize',12); + set(fig,'PaperPositionMode','auto'); + print(fig,'-dpng','-zbuffer','-r300',filename) + disp(strcat('file saved at ',filename)) +end \ No newline at end of file diff --git a/matlabCode/namesAndNumbers.m b/matlabCode/namesAndNumbers.m new file mode 100644 index 0000000..79cf196 --- /dev/null +++ b/matlabCode/namesAndNumbers.m @@ -0,0 +1,23 @@ +function [names,numbers]=namesAndNumbers(pathToFile) +% namesAndNubers returns names and numbers of subjects and runs according +% to the file given (created by bashscript) + 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/matlabCode/noOfSynergies.m b/matlabCode/noOfSynergies.m new file mode 100644 index 0000000..8d76d0c --- /dev/null +++ b/matlabCode/noOfSynergies.m @@ -0,0 +1,98 @@ +% plots performance (squared elementwise correlation) for 1 to 6 synergies for Autoencoder, PCA and NNMF +addpath('/home/hohlochj/masterarbeit/usedMcode') +pathToFile='/nfs/wsi/ti/messor/hohlochj/origData/'; +% pathToFile='/home/jph/Uni/masterarbeit/origData/'; +showFigure=true; +oneSubject=true; + +EMGChannels={'AbdPolLo','Biceps','Triceps','FrontDelt','MidDelt','BackDelt'}; +maxFile=5; +windowEMG=0.2; +windowShift=0.05; + +maxSize=6; + + +[subjects,numbers]=namesAndNumbers(pathToFile); +numbersMat=cell2mat(numbers); +subjectsForNumbers=cell(size(numbersMat,2),1); +j=1; +for i=1:size(subjects,2) + subject=subjects{i}; + for number=numbers{i} + subjectsForNumbers{j}=subject; + j=j+1; + end +end + +if oneSubject + numbersMat=numbersMat(1); + number=numbersMat(1); + subject=subjectsForNumbers{1}; +end + +EMG=cell(size(numbersMat,2),1); + +parfor j=1:size(numbersMat,2) + EMG{j}=readEMGSig(pathToFile,subjectsForNumbers{j},numbersMat(j),... + maxFile,windowEMG,windowShift,EMGChannels); +end + +trainingData=cell2mat(EMG); + + +r2Autoenc=zeros([maxSize,size(trainingData,2)]); +ae=cell([maxSize,1]); +parfor i=1:maxSize + ae{i}=trainAutoencoder(trainingData',i,'ShowProgressWindow',false); + predicted=predict(ae{i},trainingData'); + r2Autoenc(i,:)=correlation2(trainingData,predicted'); +end + +%PCA +[COEFF,SCORE,~] = pca(trainingData,'Centered',false); +r2PCA=zeros([maxSize,size(trainingData,2)]); +parfor i=1:maxSize + predicted=SCORE(:,1:i)*COEFF(:,1:i)'; + r2PCA(i,:)=correlation2(trainingData,predicted); +end + +%NNMF +r2NNMF=zeros([maxSize,size(trainingData,2)]); +parfor i=1:maxSize + [W,H]=nnmf(trainingData,i); + predicted=W*H; + r2NNMF(i,:)=correlation2(trainingData,predicted); +end + +if showFigure + legendNames=EMGChannels; + legendNames{size(EMGChannels,2)+1}='mean'; + fig=figure(); + subplot(1,3,1); + plot([r2Autoenc,mean(r2Autoenc,2)]) + legend(legendNames','Location','SouthEast') + title('Autoencoder') + xlabel('number of syergies') + ylabel('R^2') + + subplot(1,3,2); + plot([r2PCA,mean(r2PCA,2)]) + legend(legendNames','Location','SouthEast') + title('PCA') + xlabel('number of syergies') + ylabel('R^2') + + subplot(1,3,3); + plot([r2NNMF,mean(r2NNMF,2)]) + legend(legendNames','Location','SouthEast') + title('NNMF') + xlabel('number of syergies') + ylabel('R^2') +else + if oneSubject + save(strcat(pathToFile,sprintf('../matlabData/%s_noSyn%s%i.mat',datestr(datetime('now')),subject,number))); + else + save(strcat(pathToFile,sprintf('../matlabData/%s_noSynAll.mat',datestr(datetime('now'))))); + end +end \ No newline at end of file diff --git a/matlabCode/overview.m b/matlabCode/overview.m new file mode 100644 index 0000000..85cd814 --- /dev/null +++ b/matlabCode/overview.m @@ -0,0 +1,7 @@ +% outputs mean, std, min and max of declared input +input=correlations.EEG.emg.Default; + +m=mean(input) +s=std(input) +ma=max(input) +mi=min(input) \ No newline at end of file diff --git a/matlabCode/pickFromStruct.m b/matlabCode/pickFromStruct.m new file mode 100644 index 0000000..e268217 --- /dev/null +++ b/matlabCode/pickFromStruct.m @@ -0,0 +1,10 @@ +function [ output ] = pickFromStruct( struct,indices ) +%pickFromStruct picks ith entry of a struct + names=fieldnames(struct); + pick=names(indices); + output=cell([size(pick,1),1]); + for i=1:size(pick,1) + output{i}=struct.(sprintf('%s',pick{i})); + end +end + diff --git a/matlabCode/plotCM.m b/matlabCode/plotCM.m new file mode 100644 index 0000000..fd9d69b --- /dev/null +++ b/matlabCode/plotCM.m @@ -0,0 +1,23 @@ +function [fig] = plotCM(cm,figureTitle) +% plotCM plots a confusion matrix with the given title + imagesc(squeeze(mean(cm,1))) + colorbar(); + fig=gcf; + axes =gca; + if size(cm,2)==5 + axes.YAxis.TickLabels={'Rest','1','2','3','4'}; + axes.YAxis.TickValues=[1,2,3,4,5]; + axes.XAxis.TickLabels={'Rest','1','2','3','4'}; + axes.XAxis.TickValues=[1,2,3,4,5]; + elseif size(cm,2)==4 + axes.YAxis.TickLabels={'1','2','3','4'}; + axes.YAxis.TickValues=[1,2,3,4]; + axes.XAxis.TickLabels={'1','2','3','4'}; + axes.XAxis.TickValues=[1,2,3,4]; + else + error('CM has neither 5 nor 4 entries - change plotCM to match your needs') + end + xlabel('classified as') + ylabel('actual class') + title(figureTitle) +end \ No newline at end of file diff --git a/matlabCode/psdPlot.m b/matlabCode/psdPlot.m new file mode 100644 index 0000000..af7b7a7 --- /dev/null +++ b/matlabCode/psdPlot.m @@ -0,0 +1,14 @@ +% compares pwelch and pburg +[sig, state, params] = load_bcidat(strcat(pathToFile,... + sprintf('%s/%s_B100%i/%s_B1S00%iR02.dat',subject,subject,... + number,subject,number))); + +fig=figure(); + +subplot(2,1,1); +pwelch(sig(:,14),[],[],2048,2500); +xlim([0.01,0.5]) + +subplot(2,1,2); +pburg(double(sig(:,14)),pburgOrder,2048,2500); +xlim([0.01,0.5]) \ No newline at end of file diff --git a/matlabCode/readAll.m b/matlabCode/readAll.m new file mode 100644 index 0000000..1d18cab --- /dev/null +++ b/matlabCode/readAll.m @@ -0,0 +1,95 @@ +function [savePath]=readAll(pathToFile, subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,... + threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels,noSynergies) +% readAll reads the data for given subject and run and saves it with the +% according parameters; savePath is returned +% if the file exists, nothing is done + %fprintf('start read %s%i %s\n',subject,number,datestr(datetime('now'))); + savePath=strcat(pathToFile,... + sprintf('../matlabData/%s%i%imsWindowEMG%isWindowEEG%imsShiftEMG%imsShiftEEGFreq%ito%iPause%ipBurg%inoSyn%i.mat',... + subject,number,windowEMG*1000,windowEEG,shiftEMG*1000,shiftEEG*1000,minEEGFreq,maxEEGFreq,pause,pburgOrder,noSynergies)); + + %only create file if it doesn't exist yet + if exist(savePath, 'file') ~= 2 + fprintf(strcat(savePath,' not existing; creating\n')); + trainingDataEEGcell=cell(maxFile,1); + trainingDataEEGlfCell=cell(maxFile,1); + trainingDataEMGcell=cell(maxFile,1); + classesCell=cell(maxFile,1); + kin=cell([maxFile,1]); + + % accumulate preprocessed data in cells + for i=1:maxFile + [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},trainingDataEEGlfCell{i},... + trainingDataEMGcell{i},kin{i}]=generateTrainingData(sig,... + tmpKin(:,1:4),windowEMG,windowEEG,shiftEMG,shiftEEG,params,pburgOrder,... + minEEGFreq,maxEEGFreq,noLFsamples,EMGChannels); + classesCell{i}=stat.StimulusCode; + %fprintf('%ith file processed\n',i) + end + bci_sf=params.SamplingRate.NumericValue; + + clear sig + + % reformat to matrices + trainingDataEEG=cell2mat(trainingDataEEGcell); + trainingDataEEGlf=cell2mat(trainingDataEEGlfCell); + trainingDataEMG=cell2mat(trainingDataEMGcell); + classesMat=cell2mat(classesCell); + kinMat=cell2mat(kin); + clear trainingDataEEGcell trainingDataEMGcell classesCell kin + + % reclassification according to EMG + classificationWithPause=classifyAccordingToEMG(trainingDataEMG,... + classesMat,shiftEMG,bci_sf,threshold,pause); + clear classesMat + + % Smoother classifcation + smoothClassificationEMG=zeros(size(classificationWithPause)); + for i=1:size(classificationWithPause,1) + smoothClassificationEMG(i)=round(mode(classificationWithPause(... + max(i-fix(bci_sf/1000),1):min(i+fix(bci_sf/1000),end)))); + end + + clear classificationWithPause + + %adjust to length of EEG data + smoothClassificationEEG=zeros([size(trainingDataEEG,1),1]); + for i=1:size(smoothClassificationEEG,1) + smoothClassificationEEG(i)=mode(smoothClassificationEMG(fix(... + round(shiftEEG/shiftEMG*(i-1))+1):min(fix(round(... + shiftEEG/shiftEMG*i)),max(size(smoothClassificationEMG))))); + end + + % drop data according to classification + trainingDataEEG=trainingDataEEG(smoothClassificationEEG~=-1,:,:); + trainingDataEEGlf=trainingDataEEGlf(smoothClassificationEEG~=-1,:,:); + trainingDataEMG=trainingDataEMG(smoothClassificationEMG~=-1,:); + classificationEMG=smoothClassificationEMG(smoothClassificationEMG~=-1); + classificationEEG=smoothClassificationEEG(smoothClassificationEEG~=-1); + + % all created matching the size of EMGdata, if necessary downsampling + % has to be done when used + kinematics=kinMat(smoothClassificationEMG~=-1,:); + ae=trainAutoencoder(trainingDataEMG',noSynergies,'ShowProgressWindow',false); + synergiesAutoenc=encode(ae,trainingDataEMG')'; + COEFF=pca(trainingDataEMG,'Centered',false); + synergiesPCA=trainingDataEMG/COEFF(:,1:noSynergies)'; + [~,H]=nnmf(trainingDataEMG,noSynergies); + synergiesNNMF=trainingDataEMG*H'; + + clear smoothClassification i + save(savePath,'trainingDataEEG','trainingDataEEGlf','trainingDataEMG',... + 'classificationEMG','classificationEEG','kinematics','synergiesAutoenc',... + 'synergiesPCA','synergiesNNMF','-v7.3'); + + + %fprintf('finished reading %s%i %s\n',subject,number,datestr(datetime('now'))); + else + fprintf(strcat(savePath,' yet existing; nothing to do\n')); + end +end diff --git a/matlabCode/readAllPos.m b/matlabCode/readAllPos.m new file mode 100644 index 0000000..92d214f --- /dev/null +++ b/matlabCode/readAllPos.m @@ -0,0 +1,95 @@ +function [savePath]=readAllPos(pathToFile, subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,... + threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels,noSynergies) +% readAllPos reads the data for given subject and run and saves it with the +% according parameters; savePath is returned +% if the file exists, nothing is done + %fprintf('start read %s%i %s\n',subject,number,datestr(datetime('now'))); + savePath=strcat(pathToFile,... + sprintf('../matlabData/%s%i%imsWindowEMG%isWindowEEG%imsShiftEMG%imsShiftEEGFreq%ito%iPause%ipBurg%inoSyn%iPos.mat',... + subject,number,windowEMG*1000,windowEEG,shiftEMG*1000,shiftEEG*1000,minEEGFreq,maxEEGFreq,pause,pburgOrder,noSynergies)); + + %only create file if it doesn't exist yet + if exist(savePath, 'file') ~= 2 + fprintf(strcat(savePath,' not existing; creating\n')); + trainingDataEEGcell=cell(maxFile,1); + trainingDataEEGlfCell=cell(maxFile,1); + trainingDataEMGcell=cell(maxFile,1); + classesCell=cell(maxFile,1); + kin=cell([maxFile,1]); + + % accumulate preprocessed data in cells + for i=1:maxFile + [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},trainingDataEEGlfCell{i},... + trainingDataEMGcell{i},kin{i}]=generateTrainingDataPos(sig,... + tmpKin(:,1:4),windowEMG,windowEEG,shiftEMG,shiftEEG,params,pburgOrder,... + minEEGFreq,maxEEGFreq,noLFsamples,EMGChannels); + classesCell{i}=stat.StimulusCode; + %fprintf('%ith file processed\n',i) + end + bci_sf=params.SamplingRate.NumericValue; + + clear sig + + % reformat to matrices + trainingDataEEG=cell2mat(trainingDataEEGcell); + trainingDataEEGlf=cell2mat(trainingDataEEGlfCell); + trainingDataEMG=cell2mat(trainingDataEMGcell); + classesMat=cell2mat(classesCell); + kinMat=cell2mat(kin); + clear trainingDataEEGcell trainingDataEMGcell classesCell kin + + % reclassification according to EMG + classificationWithPause=classifyAccordingToEMG(trainingDataEMG,... + classesMat,shiftEMG,bci_sf,threshold,pause); + clear classesMat + + % Smoother classifcation + smoothClassificationEMG=zeros(size(classificationWithPause)); + for i=1:size(classificationWithPause,1) + smoothClassificationEMG(i)=round(mode(classificationWithPause(... + max(i-fix(bci_sf/1000),1):min(i+fix(bci_sf/1000),end)))); + end + + clear classificationWithPause + + %adjust to length of EEG data + smoothClassificationEEG=zeros([size(trainingDataEEG,1),1]); + for i=1:size(smoothClassificationEEG,1) + smoothClassificationEEG(i)=mode(smoothClassificationEMG(fix(... + round(shiftEEG/shiftEMG*(i-1))+1):min(fix(round(... + shiftEEG/shiftEMG*i)),max(size(smoothClassificationEMG))))); + end + + % drop data according to classification + trainingDataEEG=trainingDataEEG(smoothClassificationEEG~=-1,:,:); + trainingDataEEGlf=trainingDataEEGlf(smoothClassificationEEG~=-1,:,:); + trainingDataEMG=trainingDataEMG(smoothClassificationEMG~=-1,:); + classificationEMG=smoothClassificationEMG(smoothClassificationEMG~=-1); + classificationEEG=smoothClassificationEEG(smoothClassificationEEG~=-1); + + % all created matching the size of EMGdata, if necessary downsampling + % has to be done when used + kinematics=kinMat(smoothClassificationEMG~=-1,:); + ae=trainAutoencoder(trainingDataEMG',noSynergies,'ShowProgressWindow',false); + synergiesAutoenc=encode(ae,trainingDataEMG')'; + COEFF=pca(trainingDataEMG,'Centered',false); + synergiesPCA=trainingDataEMG/COEFF(:,1:noSynergies)'; + [~,H]=nnmf(trainingDataEMG,noSynergies); + synergiesNNMF=trainingDataEMG*H'; + + clear smoothClassification i + save(savePath,'trainingDataEEG','trainingDataEEGlf','trainingDataEMG',... + 'classificationEMG','classificationEEG','kinematics','synergiesAutoenc',... + 'synergiesPCA','synergiesNNMF','-v7.3'); + + + %fprintf('finished reading %s%i %s\n',subject,number,datestr(datetime('now'))); + else + fprintf(strcat(savePath,' yet existing; nothing to do\n')); + end +end diff --git a/matlabCode/readEMGSig.m b/matlabCode/readEMGSig.m new file mode 100644 index 0000000..e58b117 --- /dev/null +++ b/matlabCode/readEMGSig.m @@ -0,0 +1,13 @@ +function [EMG]=readEMGSig(pathToFile,subject,number,maxFile,windowEMG,windowShift,EMGChannels) +% readEMGSig reads the EMG only + EMGcell=cell(maxFile,1); + + for i=1:maxFile + [sig, ~, params] = load_bcidat(strcat(pathToFile,... + sprintf('%s/%s_B100%i/%s_B1S00%iR0%i.dat',subject,subject,number,subject,number,i))); + EMGcell{i}=waveformLengthAll(sig(:,ismember(params.ChannelNames.Value,EMGChannels)),... + params.SamplingRate.NumericValue,windowEMG,windowShift); + end + + EMG=cell2mat(EMGcell); +end \ No newline at end of file diff --git a/matlabCode/ridgeCV.m b/matlabCode/ridgeCV.m new file mode 100644 index 0000000..7b06653 --- /dev/null +++ b/matlabCode/ridgeCV.m @@ -0,0 +1,88 @@ +function [correlation,maxRidgeParamIndex]=ridgeCV(savePath,data,goal,k,ridgeParams,eegOffset) +% ridgeCV runs a ridge regression to predict 'goal' from 'data' + load(savePath); + clear classification; %classification is not need + eeg=false; + if strcmp(data,'EEG') + eeg=true; + trainingData=trainingDataEEG(eegOffset+1:end,:); + elseif strcmp(data,'EMG') + trainingData=trainingDataEMG; + elseif strcmp(data,'LF') + eeg=true; + trainingData=trainingDataEEGlf(eegOffset+1:end,:); + elseif strcmp(data,'Autoenc') + trainingData=synergiesAutoenc; + elseif strcmp(data,'PCA') + trainingData=synergiesPCA; + elseif strcmp(data,'NNMF') + trainingData=synergiesNNMF; + else + error('only EEG, EMG, LF, Autoenc, PCA and NNMF are valid inputs for data'); + end + % if data is smaller then goal, goal is scaled down + factor=size(trainingDataEMG,1)/(size(trainingData,1)+eegOffset); + + + if strcmp(goal,'Autoenc') + predicted=shiftingMean(synergiesAutoenc,factor); + elseif strcmp(goal,'PCA') + predicted=shiftingMean(synergiesPCA,factor); + elseif strcmp(goal,'NNMF') + predicted=shiftingMean(synergiesNNMF,factor); + elseif strcmp(goal,'kin') + predicted=shiftingMean(kinematics,factor); + elseif strcmp(goal,'EMG') + predicted=shiftingSum(trainingDataEMG,factor); + else + error('only kin, Autoenc, PCA, NNMF and EMG are valid inputs for goal'); + end + + clear trainingDataEEG; + clear trainingDataEEGlf; + clear trainingDataEMG; + clear synergies*; + clear kinematics; + + % is eeg Offset applied? + if eeg + predicted=predicted(1:end-eegOffset,:); + end + + %disp([size(trainingData),size(predicted)]) + %adjust small errors when unsung down scaling + minSize=min([size(trainingData,1),size(predicted,1)]); + if minSize+5 via + viaDataPredictedLeave=viaDataPredicted(mod(randMap,k)==i-1,:); + viaDataPredictedRemaining=viaDataPredicted(mod(randMap,k)~=i-1,:); + + for j=1:size(predicted,2) + + coeffs=ridge(remainingPred(:,j),viaDataPredictedRemaining,ridgeParams(1),0); + + finalCorr(i,j)=ridgeCorrelation(viaDataPredictedLeave,leavePred(:,j),coeffs); + + end + end + viaCorrelation=mean(viaCorr,1); + correlation=mean(finalCorr,1); +end diff --git a/matlabCode/ridgeCorrelation.m b/matlabCode/ridgeCorrelation.m new file mode 100644 index 0000000..b4ac947 --- /dev/null +++ b/matlabCode/ridgeCorrelation.m @@ -0,0 +1,7 @@ +function [correlation]=ridgeCorrelation(data,kin,coeffs) +% ridgeCorrelation calculates the correlation with given coefficients for +% ridge regression + prediction=coeffs(1)+data*coeffs(2:end); + + correlation=corr(kin,prediction); +end \ No newline at end of file diff --git a/matlabCode/shiftingDownsample.m b/matlabCode/shiftingDownsample.m new file mode 100644 index 0000000..e923b60 --- /dev/null +++ b/matlabCode/shiftingDownsample.m @@ -0,0 +1,11 @@ +function [out]=shiftingDownsample(signal,frequency,windowEEG,shift,noSamples) +% shiftingDownsample downsamples in shifting windows + signalShift=frequency*shift; + signalWindow=frequency*windowEEG; + shiftProp=windowEEG/shift; + noOfWindows=fix(floor(size(signal,1)/signalWindow-1)*shiftProp+1); + out=zeros([noOfWindows,noSamples]); + for j=1:noOfWindows + out(j,:)=myDownsample(signal((j-1)*signalShift+1:(j-1)*signalShift+signalWindow),noSamples); + end +end \ No newline at end of file diff --git a/matlabCode/shiftingKin.m b/matlabCode/shiftingKin.m new file mode 100644 index 0000000..b5a14fb --- /dev/null +++ b/matlabCode/shiftingKin.m @@ -0,0 +1,13 @@ +function [kinPerSec]=shiftingKin(kin, windowEMG, shiftEMG) +% calculates velocity from kinematic data in shifting window + kinPerSec=zeros(fix((max(kin(:,1))-windowEMG*1000)/(shiftEMG*1000)),3); + for j=1:size(kinPerSec,1) + tmp=sum(diff(kin(kin(:,1)>(j-1)*shiftEMG*1000 & kin(:,1)<=(j-1)*shiftEMG*1000+windowEMG*1000,2:4))); + i=j; + while isnan(tmp) %interval is empty + tmp=sum(diff(kin(kin(:,1)>(i-1)*shiftEMG*1000 & kin(:,1)<=(i-1)*shiftEMG*1000+windowEMG*1000,2:4))); + i=i+1; + end + kinPerSec(j,:)=tmp; + end +end \ No newline at end of file diff --git a/matlabCode/shiftingMean.m b/matlabCode/shiftingMean.m new file mode 100644 index 0000000..edc2517 --- /dev/null +++ b/matlabCode/shiftingMean.m @@ -0,0 +1,10 @@ +function [ output ] = shiftingMean( data,factor ) +%shiftingMean downsampling by taking the mean +% taking the mean of $factor data as new datum + dataSize=size(data); + output=zeros([fix(dataSize(1)/factor),dataSize(2:end)]); + for i=1:size(output,1) + output(i,:)=mean(data(fix(round((i-1)*factor))+1:min(fix(round(i*factor)),dataSize(1)),:),1); + end +end + diff --git a/matlabCode/shiftingPburg.m b/matlabCode/shiftingPburg.m new file mode 100644 index 0000000..c6c1daf --- /dev/null +++ b/matlabCode/shiftingPburg.m @@ -0,0 +1,12 @@ +function [ EEG ] = shiftingPburg( signal,frequency,windowEEG,shift,order,minEEGFreq,maxEEGFreq) +%shiftingPburg uses pburg to calculate the power spectrum + signalShift=frequency*shift; + signalWindow=frequency*windowEEG; + shiftProp=windowEEG/shift; + noOfWindows=fix(floor(size(signal,1)/signalWindow-1)*shiftProp+1); + EEG=zeros([noOfWindows,fix(maxEEGFreq-minEEGFreq+1)]); + for j=1:noOfWindows + [EEG(j,:),~]=pburg(signal((j-1)*signalShift+1:(j-1)*signalShift+signalWindow),order,minEEGFreq:maxEEGFreq,frequency); + end +end + diff --git a/matlabCode/shiftingPos.m b/matlabCode/shiftingPos.m new file mode 100644 index 0000000..c8701a1 --- /dev/null +++ b/matlabCode/shiftingPos.m @@ -0,0 +1,13 @@ +function [kinPerSec]=shiftingPos(kin, windowEMG, shiftEMG) +% shiftingPos calculates position in window by taking the mean + kinPerSec=zeros(fix((max(kin(:,1))-windowEMG*1000)/(shiftEMG*1000)),3); + for j=1:size(kinPerSec,1) + tmp=mean(kin(kin(:,1)>(j-1)*shiftEMG*1000 & kin(:,1)<=(j-1)*shiftEMG*1000+windowEMG*1000,2:4)); + i=j; + while isnan(tmp) %interval is empty, next existen data will be used + tmp=mean(kin(kin(:,1)>i*shiftEMG*1000 & kin(:,1)<=i*shiftEMG*1000+windowEMG*1000,2:4)); + i=i+1; + end + kinPerSec(j,:)=tmp; + end +end \ No newline at end of file diff --git a/matlabCode/shiftingSum.m b/matlabCode/shiftingSum.m new file mode 100644 index 0000000..9495598 --- /dev/null +++ b/matlabCode/shiftingSum.m @@ -0,0 +1,10 @@ +function [ output ] = shiftingSum( data,factor ) +%shiftingSum downsampling by taking the sum +% taking the sum of $factor data as new datum + dataSize=size(data); + output=zeros([fix(dataSize(1)/factor),dataSize(2:end)]); + for i=1:size(output,1) + output(i,:)=sum(data(fix(round((i-1)*factor))+1:min(fix(round(i*factor)),dataSize(1)),:),1); + end +end + diff --git a/matlabCode/svmEciton.m b/matlabCode/svmEciton.m new file mode 100644 index 0000000..a55d904 --- /dev/null +++ b/matlabCode/svmEciton.m @@ -0,0 +1,56 @@ +function [meanAccurancy,maxC, cmScaled]= svmEciton(savePath,EEG,k,maxExpC,maxPerClass,eegOffset,moveRest) +% svmEciton runs SVM in parallel to decide class for data with different Cs + load(savePath) + % fprintf('%i,%i,%i',size(trainingDataEMG,1),size(trainingDataEEG,1),size(classification,1)) + addpath('/nfs/wsi/ti/messor/hohlochj/libsvm/matlab'); + + %choose to estimate based on EEG or EMG + if strcmp(EEG,'EEG') + trainingData=trainingDataEEG(1:end-eegOffset,:); + classification=classificationEEG(eegOffset+1:end); + elseif strcmp(EEG,'EMG') + trainingData=trainingDataEMG; + classification=classificationEMG; + elseif strcmp(EEG,'LF') + trainingData=trainingDataEEGlf(1:end-eegOffset,:); + classification=classificationEEG(eegOffset+1:end); + else + error('only EEG, EMG and LF are valid inputs'); + end + + clear trainingDataE*; + clear classifiactionE*; + clear kinematics; + clear synergies*; + + % only intersting whether moving or resting + if moveRest + classification=double(classification~=0); + end + + accuracy=zeros(k,3); + maxC=zeros(k,1); + noClasses=size(unique(classification),1); + cm=zeros(noClasses); + randMap=randperm(size(trainingData,1)); + %disp('startCV') + parfor i=1:k %k-fold CV + % split data into training and test set + leaveOut=trainingData(mod(randMap,k)==i-1,:,:); + leaveClasses=classification(mod(randMap,k)==i-1); + remaining=trainingData(mod(randMap,k)~=i-1,:,:); + remainingClasses=classification(mod(randMap,k)~=i-1); + %fprintf('%s create %ith model\n',datestr(datetime('now')),i) + + [model,maxC(i)]=kfoldCV(remainingClasses,remaining,k,maxExpC,maxPerClass); + + [predictions,accuracy(i,:),~]=svmpredict(leaveClasses,leaveOut(:,:),model); + cm=cm+confusionmat(leaveClasses,predictions); %confusion matrix + end + meanAccurancy=mean(accuracy(:,1)); % accuracy is the mean of accuracys + + cmScaled=zeros(size(cm)); % scale confusion matrix to 0-1 + for i=1:size(cm,1) + cmScaled(i,:)=cm(i,:)/sum(cm(i,:)); + end +end diff --git a/matlabCode/topoplot/electrocaps_all.locs b/matlabCode/topoplot/electrocaps_all.locs new file mode 100644 index 0000000..9b71dee --- /dev/null +++ b/matlabCode/topoplot/electrocaps_all.locs @@ -0,0 +1,32 @@ +1 -18 .352 Fp1. +2 18 .352 Fp2. +3 -46 0.50556 F7.. +4 -38.574506 0.315310 F3.. +5 -3.871021 0.249176 FZ.. +6 34.612754 0.286614 F4.. +7 47.267805 0.425391 F8.. +8 -66 0.38889 FC5. +9 -36 0.16667 FC1. +10 36 0.16667 FC2. +11 66 0.38889 FC6. +12 -90 0.51667 T7.. +13 -86.767948 0.218632 C3.. +14 27.092627 0.014595 Cz.. +15 85.892653 0.214625 C4.. +16 87.166427 0.424522 T8.. +17 -117 0.65 TP9. +18 -115 0.36667 CP5. +19 -141 0.15556 CP1. +20 141 0.15556 CP2. +21 115 0.36667 CP6. +22 117 0.65 TP10 +23 -128.419583 0.389577 P7.. +24 -140.004936 0.253792 P3.. +25 176.354756 0.194046 Pz.. +26 136.232553 0.259532 P4.. +27 125.426297 0.398099 P8.. +28 -142 0.55 PO9. +29 -165.045704 0.369739 O1.. +30 178.207328 0.368312 OZ.. +31 162.455749 0.376674 O2.. +32 142 0.55 PO10 diff --git a/matlabCode/topoplot/plotOneSubjectOneDay.m b/matlabCode/topoplot/plotOneSubjectOneDay.m new file mode 100644 index 0000000..d6b5f68 --- /dev/null +++ b/matlabCode/topoplot/plotOneSubjectOneDay.m @@ -0,0 +1,66 @@ +function plotOneSubjectOneDay() + pathToFile='/home/jph/Uni/masterarbeit/origData/'; + eloc_file='/home/jph/Uni/masterarbeit/topoplot/electrocaps_all.locs'; + maxFile=5; + stepSize=1; + +% [subjects,numbers]=namesAndNumbers(pathToFile); +% numbersMat=cell2mat(numbers); +% subjectsForNumbers=cell(size(numbersMat,2),1); +% j=1; +% for i=1:size(subjects,2) +% subject=subjects{i}; +% for number=numbers{i} +% subjectsForNumbers{j}=subject; +% j=j+1; +% end +% end +% j=fix(rand()*size(numbersMat,2)+1); +% number=numbersMat(j); +% subject=subjectsForNumbers{j}; +% fprintf('%s%i',subject,number) + + subject='FS'; + number=3; + + maxFile=5; + threshold=10000; %EMG is classified as movement + EMGChannels={'AbdPolLo','Biceps','Triceps','FrontDelt','MidDelt','BackDelt'}; + noSynergies=3; + noLFsamples=5; + pause=false; + pburgOrder=250; + windowEMG=0.2; + windowEEG=1; + shiftEMG=0.05; + shiftEEG=0.2; + minEEGFreq=2; + maxEEGFreq=49; + + + savePath=readAll(pathToFile,subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,... + threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels,noSynergies); + + load(savePath) + clear trainingDataEMG; + clear trainingDataEEGlf; + clear classifiactionEMG; + clear kinematics; + clear synergies*; + + meanRest=squeeze(mean(trainingDataEEG(classificationEEG==0,:,minEEGFreq:maxEEGFreq),1)); + meanMove=squeeze(mean(trainingDataEEG(classificationEEG>0,:,minEEGFreq:maxEEGFreq),1)); + + samplingFrequency=1; + + %alpha band + meanRest=mean(meanRest(:,7:13),2); + meanMove=mean(meanMove(:,7:13),2); + + %beta band + meanRest=mean(meanRest(:,13:20),2); + meanMove=mean(meanMove(:,13:20),2); + + + topoplot_Wrapper((meanMove./meanRest - 1)',minEEGFreq:maxEEGFreq,eloc_file,stepSize,samplingFrequency) +end \ No newline at end of file diff --git a/matlabCode/topoplot/topoplot_Wrapper.m b/matlabCode/topoplot/topoplot_Wrapper.m new file mode 100644 index 0000000..83e26a1 --- /dev/null +++ b/matlabCode/topoplot/topoplot_Wrapper.m @@ -0,0 +1,38 @@ +function topoplot_Wrapper(EEG,codes,eloc_file,stepSize,samplingFrequency) +% Wrapper for topoplot_martin() +% Allows for plottig all EEG-Data of one timeline and sliding through them +% +% EEG - consecutive EEG data (no of channels has to match the number in eloc_file) +% +% eloc_file - path to the eloc-file, type topoplot_martin('demo') for an example +% +% stepSize - number of Datapoints the slider skips in one step (1 for none) +% +% samplingFrequency - sampling frequency for the EEG data + + +slider=uicontrol('Style','slider','Min',1,'Max',fix(size(EEG,1)/stepSize),... + 'Value',1,'SliderStep',[stepSize/size(EEG,1),stepSize*10/size(EEG,1)],... + 'Callback',{@topoplotSlide,EEG,codes,eloc_file,stepSize,samplingFrequency},... + 'Position',[10 10 450 20]); + +timeField=uicontrol('Style','edit','Callback',... + {@topoplotText,EEG,codes,eloc_file,stepSize,samplingFrequency},... + 'Position',[460 10 50 20]); +text=uicontrol('Style','text','Position',[510 10 30 15]); +set(text,'String','Hz'); + + function topoplotSlide(source,~,EEG,codes,eloc_file,stepSize,samplingFrequency) + idx=fix(stepSize*source.Value); + set(timeField,'String',num2str(fix(idx/samplingFrequency))); + %disp(codes(idx)) + topoplot_martin(double(EEG(idx,:)),eloc_file) + end + + function topoplotText(source,~,EEG,codes,eloc_file,stepSize,samplingFrequency) + idx=fix(str2double(source.String)*samplingFrequency); + set(slider,'Value',fix(idx/stepSize)); + %disp(codes(idx)) + topoplot_martin(double(EEG(idx,:)),eloc_file) + end +end \ No newline at end of file diff --git a/matlabCode/topoplot/topoplot_martin.m b/matlabCode/topoplot/topoplot_martin.m new file mode 100644 index 0000000..68a3296 --- /dev/null +++ b/matlabCode/topoplot/topoplot_martin.m @@ -0,0 +1,340 @@ +% topoplot() - plot a topographic map of an EEG field as a 2-D +% circular view (looking down at the top of the head) +% using cointerpolation on a fine cartesian grid. +% Usage: +% >> topoplot(datavector,'eloc_file'); +% >> topoplot(datavector,'eloc_file', 'Param1','Value1', ...) +% Inputs: +% datavector = vector of values at the corresponding locations. +% 'eloc_file' = name of an EEG electrode position file {0 -> 'chan_file'} +% +% Optional Parameters & Values (in any order): +% Param Value +% 'colormap' - any sized colormap +% 'interplimits' - 'electrodes' to furthest electrode +% 'head' to edge of head +% {default 'head'} +% 'gridscale' - scaling grid size {default 67} +% 'maplimits' - 'absmax' +/- the absolute-max +% 'maxmin' scale to data range +% [clim1,clim2] user-definined lo/hi +% {default = 'absmax'} +% 'style' - 'straight' colormap only +% 'contour' contour lines only +% 'both' both colormap and contour lines +% 'fill' constant color between lines +% 'blank' just head and electrodes +% {default = 'both'} +% 'numcontour' - number of contour lines +% {default = 6} +% 'shading' - 'flat','interp' {default = 'flat'} +% 'headcolor' - Color of head cartoon {default black} +% 'electrodes' - 'on','off','labels','numbers' +% 'efontsize','electcolor','emarker','emarkersize' - details +% +% Note: topoplot() only works when map limits are >= the max and min +% interpolated data values. +% Eloc_file format: +% chan_number degrees radius reject_level amp_gain channel_name +% (Angle-0 = Cz-to-Fz; C3-angle =-90; Radius at edge of image = 0.5) +% +% For a sample eloc file: >> topoplot('example') +% +% Note: topoplot() will ignore any electrode with a position outside +% the head (radius > 0.5) + +% Topoplot Version 2.1 + +% Begun by Andy Spydell, NHRC, 7-23-96 +% 8-96 Revised by Colin Humphries, CNL / Salk Institute, La Jolla CA +% -changed surf command to imagesc (faster) +% -can now handle arbitrary scaling of electrode distances +% -can now handle non integer angles in eloc_file +% 4-4-97 Revised again by Colin Humphries, reformat by SM +% -added parameters +% -changed eloc_file format +% 2-26-98 Revised by Colin +% -changed image back to surface command +% -added fill and blank styles +% -removed extra background colormap entry (now use any colormap) +% -added parameters for electrode colors and labels +% -now each topoplot axes use the caxis command again. +% -removed OUTPUT parameter +% 3-11-98 changed default emarkersize, improve help msg -sm + +function handle = topoplot_martin(Vl,loc_file,p1,v1,p2,v2,p3,v3,p4,v4,p5,v5,p6,v6,p7,v7,p8,v8,p9,v9) + +% User Defined Defaults: +MAXCHANS = 256; +DEFAULT_ELOC = 'eloc64.txt'; +INTERPLIMITS = 'head'; % head, electrodes +MAPLIMITS = 'absmax'; % absmax, maxmin, [values] +GRID_SCALE = 67; +CONTOURNUM = 6; +STYLE = 'both'; % both,straight,fill,contour,blank +HCOLOR = [0 0 0]; +ECOLOR = [0 0 0]; +CONTCOLOR = [0 0 0]; +ELECTROD = 'labels'; % ON OFF LABEL +EMARKERSIZE = 6; +EFSIZE = get(0,'DefaultAxesFontSize'); +HLINEWIDTH = 2; +EMARKER = '.'; +SHADING = 'flat'; % flat or interp +LINES = []; + +%%%%%%%%%%%%%%%%%%%%%%% +nargs = nargin; +if nargs < 2 + loc_file = DEFAULT_ELOC; +end +if nargs == 1 + if isstr(Vl) + if any(strcmp(lower(Vl),{'example','demo'})) + fprintf(['This is an example of an electrode location file,\n',... + 'an ascii file consisting of the following four columns:\n',... + ' channel_number degrees arc_length channel_name\n\n',... + 'Example:\n',... + ' 1 -18 .352 Fp1.\n',... + ' 2 18 .352 Fp2.\n',... + ' 5 -90 .181 C3..\n',... + ' 6 90 .181 C4..\n',... + ' 7 -90 .500 A1..\n',... + ' 8 90 .500 A2..\n',... + ' 9 -142 .231 P3..\n',... + '10 142 .231 P4..\n',... + '11 0 .181 Fz..\n',... + '12 0 0 Cz..\n',... + '13 180 .181 Pz..\n\n',... + 'The model head sphere has a diameter of 1.\n',... + 'The vertex (Cz) has arc length 0. Channels with arc \n',... + 'lengths > 0.5 are not plotted nor used for interpolation.\n'... + 'Zero degrees is towards the nasion. Positive angles\n',... + 'point to the right hemisphere; negative to the left.\n',... + 'Channel names should each be four chars, padded with\n',... + 'periods (in place of spaces).\n']) + return + + end + end +end +if isempty(loc_file) + loc_file = 0; +end +if loc_file == 0 + loc_file = DEFAULT_ELOC; +end + +if nargs > 2 + if ~(round(nargs/2) == nargs/2) + error('topoplot(): Odd number of inputs?') + end + for i = 3:2:nargs + Param = eval(['p',int2str((i-3)/2 +1)]); + Value = eval(['v',int2str((i-3)/2 +1)]); + if ~isstr(Param) + error('topoplot(): Parameter must be a string') + end + Param = lower(Param); + switch lower(Param) + case 'colormap' + if size(Value,2)~=3 + error('topoplot(): Colormap must be a n x 3 matrix') + end + colormap(Value) + case {'interplimits','headlimits'} + if ~isstr(Value) + error('topoplot(): interplimits value must be a string') + end + Value = lower(Value); + if ~strcmp(Value,'electrodes') & ~strcmp(Value,'head') + error('topoplot(): Incorrect value for interplimits') + end + INTERPLIMITS = Value; + case 'maplimits' + MAPLIMITS = Value; + case 'gridscale' + GRID_SCALE = Value; + case 'style' + STYLE = lower(Value); + case 'numcontour' + CONTOURNUM = Value; + case 'electrodes' + ELECTROD = lower(Value); + case 'emarker' + EMARKER = Value; + case {'headcolor','hcolor'} + HCOLOR = Value; + case {'electcolor','ecolor'} + ECOLOR = Value; + case {'emarkersize','emsize'} + EMARKERSIZE = Value; + case {'efontsize','efsize'} + EFSIZE = Value; + case 'lines' + LINES = Value; + case 'shading' + SHADING = lower(Value); + if ~any(strcmp(SHADING,{'flat','interp'})) + error('Invalid Shading Parameter') + end + otherwise + error('Unknown parameter.') + end + end +end + +[r,c] = size(Vl); +if r>1 & c>1, + error('topoplot(): data should be a single vector\n'); +end +fid = fopen(loc_file); +if fid<1, + fprintf('topoplot(): cannot open eloc_file (%s).\n',loc_file); + return +end +A = fscanf(fid,'%d %f %f %s',[7 MAXCHANS]); +fclose(fid); + +A = A'; + +if length(Vl) ~= size(A,1), + fprintf(... + 'topoplot(): data vector must have the same rows (%d) as eloc_file (%d)\n',... + length(Vl),size(A,1)); + A + error(''); +end + +labels = setstr(A(:,4:7)); +idx = find(labels == '.'); % some labels have dots +labels(idx) = setstr(abs(' ')*ones(size(idx))); % replace them with spaces + +Th = pi/180*A(:,2); % convert degrees to radians +Rd = A(:,3); +ii = find(Rd <= 0.5); % interpolate on-head channels only +Th = Th(ii); +Rd = Rd(ii); +Vl = Vl(ii); +labels = labels(ii,:); + +[x,y] = pol2cart(Th,Rd); % transform from polar to cartesian coordinates +rmax = 0.5; + +ha = gca; +cla +hold on + +if ~strcmp(STYLE,'blank')&& ~strcmp(STYLE,'lines') && ~strcmp(STYLE,'lines2') + % find limits for interpolation + if strcmp(INTERPLIMITS,'head') + xmin = min(-.5,min(x)); xmax = max(0.5,max(x)); + ymin = min(-.5,min(y)); ymax = max(0.5,max(y)); + else + xmin = max(-.5,min(x)); xmax = min(0.5,max(x)); + ymin = max(-.5,min(y)); ymax = min(0.5,max(y)); + end + + xi = linspace(xmin,xmax,GRID_SCALE); % x-axis description (row vector) + yi = linspace(ymin,ymax,GRID_SCALE); % y-axis description (row vector) + + [Xi,Yi,Zi] = griddata(y,x,Vl,yi',xi,'v4'); % Interpolate data + + % Take data within head + mask = (sqrt(Xi.^2+Yi.^2) <= rmax); + ii = find(mask == 0); + Zi(ii) = NaN; + + % calculate colormap limits + m = size(colormap,1); + if isstr(MAPLIMITS) + if strcmp(MAPLIMITS,'absmax') + amin = -max(max(abs(Zi))); + amax = max(max(abs(Zi))); + elseif strcmp(MAPLIMITS,'maxmin') + amin = min(min(Zi)); + amax = max(max(Zi)); + end + else + amin = MAPLIMITS(1); + amax = MAPLIMITS(2); + end + delta = xi(2)-xi(1); % length of grid entry + + % Draw topoplot on head + if strcmp(STYLE,'contour') + contour(Xi,Yi,Zi,CONTOURNUM,'k'); + elseif strcmp(STYLE,'both') + surface(Xi-delta/2,Yi-delta/2,zeros(size(Zi)),Zi,'EdgeColor','none',... + 'FaceColor',SHADING); + contour(Xi,Yi,Zi,CONTOURNUM,'k'); + elseif strcmp(STYLE,'straight') + surface(Xi-delta/2,Yi-delta/2,zeros(size(Zi)),Zi,'EdgeColor','none',... + 'FaceColor',SHADING); + elseif strcmp(STYLE,'fill') + contourf(Xi,Yi,Zi,CONTOURNUM,'k'); + else + error('Invalid style') + end + caxis([amin amax]) % set coloraxis +end + +%% Plot lines +if strcmp(STYLE, 'lines') + plot(y(LINES),x(LINES)); +end + +if strcmp(STYLE, 'lines2') + for n=1:size(LINES,2) + if LINES(3,n)==0 + lcolor='k'; + elseif LINES(3,n)==1 + lcolor='b'; + elseif LINES(3,n)==2 + lcolor='r'; + end + + plot(y(LINES(1:2,n)),x(LINES(1:2,n)),lcolor,'LineWidth',LINES(4,n)); + end +end + +set(ha,'Xlim',[-rmax*1.3 rmax*1.3],'Ylim',[-rmax*1.3 rmax*1.3]) + +% %%% Draw Head %%%% +l = 0:2*pi/100:2*pi; +basex = .18*rmax; +tip = rmax*1.15; base = rmax-.004; +EarX = [.497 .510 .518 .5299 .5419 .54 .547 .532 .510 .489]; +EarY = [.0555 .0775 .0783 .0746 .0555 -.0055 -.0932 -.1313 -.1384 -.1199]; + +% Plot Electrodes +if strcmp(ELECTROD,'on') + hp2 = plot(y,x,EMARKER,'Color',ECOLOR,'markersize',EMARKERSIZE); +elseif strcmp(ELECTROD,'labels') + for i = 1:size(labels,1) + text(y(i),x(i),labels(i,:),'HorizontalAlignment','center',... + 'VerticalAlignment','middle','Color',ECOLOR,... + 'FontSize',EFSIZE) + end +elseif strcmp(ELECTROD,'numbers') +whos y x + for i = 1:size(labels,1) + text(y(i),x(i),int2str(i),'HorizontalAlignment','center',... + 'VerticalAlignment','middle','Color',ECOLOR,... + 'FontSize',EFSIZE) + end +end + +% Plot Head, Ears, Nose +plot(cos(l).*rmax,sin(l).*rmax,... + 'color',HCOLOR,'Linestyle','-','LineWidth',HLINEWIDTH); + +plot([.18*rmax;0;-.18*rmax],[base;tip;base],... + 'Color',HCOLOR,'LineWidth',HLINEWIDTH); + +plot(EarX,EarY,'color',HCOLOR,'LineWidth',HLINEWIDTH) +plot(-EarX,EarY,'color',HCOLOR,'LineWidth',HLINEWIDTH) + +hold off +axis off + diff --git a/matlabCode/waveformLength.m b/matlabCode/waveformLength.m new file mode 100644 index 0000000..3b927f3 --- /dev/null +++ b/matlabCode/waveformLength.m @@ -0,0 +1,12 @@ +function [WL]=waveformLength(signal,frequency,window,windowShift) +%waveformLength computes waveformLength of signal with given frequency, window and +%shift + signalShift=frequency*windowShift; + signalWindow=frequency*window; + shiftProp=window/windowShift; + WL=zeros((floor(size(signal,1)/signalWindow)-1)*shiftProp+1,1); + for j=1:fix(floor(size(signal,1)/signalWindow-1)*shiftProp+1) + WL(j)=sum(abs(diff(signal((j-1)*signalShift+1:(j-1)*signalShift+signalWindow)))); + end +end + diff --git a/matlabCode/waveformLengthAll.m b/matlabCode/waveformLengthAll.m new file mode 100644 index 0000000..04e19d4 --- /dev/null +++ b/matlabCode/waveformLengthAll.m @@ -0,0 +1,15 @@ +function [EMG]=waveformLengthAll(sig,bci_sf,windowEMG,shiftEMG) +% waveformLengthAll computes waveformlength after filtering sig + + %Filter around 50Hz and below 2 Hz + [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'); + + signalWindow=bci_sf*windowEMG; + shiftProp=windowEMG/shiftEMG; + EMG=zeros((floor(size(sig,1)/signalWindow)-1)*shiftProp+1,size(sig,2)); + parfor i=1:size(sig,2) + EMG(:,i)=waveformLength(filtfilt(double(E),double(F),filtfilt(double(C),double(D),filtfilt(double(A),double(B),double(sig(:,i))))),bci_sf,windowEMG,shiftEMG); + end +end \ No newline at end of file