diff --git a/oldMcode/meanCorrelations.m b/oldMcode/meanCorrelations.m new file mode 100644 index 0000000..2cbf202 --- /dev/null +++ b/oldMcode/meanCorrelations.m @@ -0,0 +1,23 @@ + + +% for different eegOffsets +corrEEGkin(eegOffset+1)=mean(correlationEEGkin); +corrEEGautoenc(eegOffset+1)=mean(correlationEEGautoenc); +corrEEGnnmf(eegOffset+1)=mean(correlationEEGnnmf); +corrEEGpca(eegOffset+1)=mean(correlationEEGpca); + +corrLFkin(eegOffset+1)=mean(correlationLFkin); +corrLFautoenc(eegOffset+1)=mean(correlationLFautoenc); +corrLFnnmf(eegOffset+1)=mean(correlationLFnnmf); +corrLFpca(eegOffset+1)=mean(correlationLFpca); +%end + +plot([corrEEGkin;corrEEGautoenc;corrEEGnnmf;corrEEGpca]',':o') +legend({'kin','autoenc','nnmf','pca'}) +xlabel('offset [x*0.2s]') +ylabel('correlation') + +plot([corrLFkin;corrLFautoenc;corrLFnnmf;corrLFpca]',':o') +legend({'kin','autoenc','nnmf','pca'}) +xlabel('offset [x*0.2s]') +ylabel('correlation') \ No newline at end of file diff --git a/plots/noSyn.fig b/plots/noSyn.fig new file mode 100644 index 0000000..9b8be0d --- /dev/null +++ b/plots/noSyn.fig Binary files differ diff --git a/readme.m b/readme.m index babc2a6..b18cd01 100644 --- a/readme.m +++ b/readme.m @@ -44,7 +44,10 @@ %Tipps % find, diff, xlim([0 0.1]) -%imagesc +%%imagesc +% fig=figure; +% imagesc(cmScaled) +% colorbar(); %libsvm % - svmTrain diff --git a/text/TODO.txt b/text/TODO.txt index e9cd377..6f78c2a 100644 --- a/text/TODO.txt +++ b/text/TODO.txt @@ -24,7 +24,7 @@ correlation zwischen Frequenzband und Koordinaten (r^2) (done) -topgraphisch plotten? (mail, elog-files, suchen) +topgraphisch plotten? (mail, elog-files, suchen) (done) low-frequency features (done) @@ -32,7 +32,7 @@ Titel: Movement Prediction from EEG based on Synergies Evaluation and Comparison to other methods -Prüfer: Rosenstiehl, Benda? (angefragt) +Prüfer: Rosenstiehl, Benda? (Zusage) topoplot_Wrapper (done) @@ -40,4 +40,13 @@ Synergien aus WL (50ms step 200 ms window) (done) -Synergien lernen +Synergien lernen (implementiert) + +shift echt variabel, bisher EEG und EMG gleich (done) + +offset für EEG (done) + +debug Output, Fortschritt + +Refactoring, Comments +(incl. save -append) diff --git a/text/collection.txt b/text/collection.txt index c993245..4218019 100644 --- a/text/collection.txt +++ b/text/collection.txt @@ -37,3 +37,6 @@ CV: maxRidgeParam: 1 (chosen from 0.001, 0.01, 0.1, 1) maxCexp: instabil: mean: EEG 0.1, EMG -0.7; std: EEG 1.729, EMG 1.567 + +eegOffset (0,1,2): mean(EEGautoenc) +EEG: 0.8277, 0.8245, 0.8188 diff --git a/usedMcode/callAll.m b/usedMcode/callAll.m index aa54e84..8982846 100644 --- a/usedMcode/callAll.m +++ b/usedMcode/callAll.m @@ -1,26 +1,28 @@ -addpath('/home/hohlochj/masterarbeit/usedMcode') +addpath('/home/hohlochj/masterarbeitNew/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=4; allSubjects=true; %run all subjects and days or only one random day for one random subject - +name='default'; %suffix for the output, has to be valid name for file (no space, /, ...) windowEMG=0.2; windowEEG=1; -shift=0.2; +shiftEMG=0.05; +shiftEEG=0.2; +eegOffset=0; %predict actions x*shiftEEG after EEG measurement pburgOrder=50; minEEGFreq=0; maxEEGFreq=200; pause=false; noLFsamples=5; -ridgeParams=[0.1]; +ridgeParams=[10]; k=10; maxExpC=0; maxPerClass=250; -EEG=true; -poolObj=parpool(32); +poolObj=parpool(18); [subjects,numbers]=namesAndNumbers(pathToFile); numbersMat=cell2mat(numbers); @@ -28,7 +30,7 @@ j=1; for i=1:size(subjects,2) subject=subjects{i}; - for number=numbers{i} + for number=numbers{i} subjectsForNumbers{j}=subject; j=j+1; end @@ -47,43 +49,77 @@ cmScaledEMG=zeros([j,5,5]); cmScaledEEG=zeros([j,5,5]); cmScaledLF=zeros([j,5,5]); - maxRidgeParamIndexEEG=zeros([j,3,k]); - maxRidgeParamIndexEMG=zeros([j,3,k]); - maxRidgeParamIndexLF=zeros([j,3,k]); - correlationEMG=zeros([j,3]); %x,y,angle - correlationEEG=zeros([j,3]); - correlationLF=zeros([j,3]); - + 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]); + parfor j=1:size(numbersMat,2) number=numbersMat(j); subject=subjectsForNumbers{j}; - savePath=readAll(pathToFile,subject,number,windowEMG,windowEEG,shift,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples); - [meanAccurancysEMG(j),maxCEMG(j,:),cmScaledEMG(j,:,:)]=svmEciton(savePath,'EMG',k,maxExpC,maxPerClass); - [meanAccurancysEEG(j),maxCEEG(j,:),cmScaledEEG(j,:,:)]=svmEciton(savePath,'EEG',k,maxExpC,maxPerClass); - [meanAccurancysLF(j),maxCLF(j,:),cmScaledLF(j,:,:)]=svmEciton(savePath,'LF',k,maxExpC,maxPerClass); - [correlationEMG(j,:),maxRidgeParamIndexEMG(j,:,:)]=ridgeCV(savePath,'EMG',k,ridgeParams); - [correlationEEG(j,:),maxRidgeParamIndexEEG(j,:,:)]=ridgeCV(savePath,'EEG',k,ridgeParams); - [correlationLF(j,:),maxRidgeParamIndexLF(j,:,:)]=ridgeCV(savePath,'LF',k,ridgeParams); + savePath=readAll(pathToFile,subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples); + [meanAccurancysEMG(j),maxCEMG(j,:),cmScaledEMG(j,:,:)]=svmEciton(savePath,'EMG',k,maxExpC,maxPerClass,eegOffset); + [meanAccurancysEEG(j),maxCEEG(j,:),cmScaledEEG(j,:,:)]=svmEciton(savePath,'EEG',k,maxExpC,maxPerClass,eegOffset); + [meanAccurancysLF(j),maxCLF(j,:),cmScaledLF(j,:,:)]=svmEciton(savePath,'LF',k,maxExpC,maxPerClass,eegOffset); + [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); fprintf('%s%i finished %s\n',subject,number,datestr(datetime('now'))) end - - save(strcat(pathToFile,sprintf('../matlabData/%s_callAll.mat',datestr(datetime('now'))))); + + save(strcat(pathToFile,sprintf('../matlabData/%s_callAll-%s.mat',datestr(datetime('now')),name))); else j=fix(rand()*size(numbersMat,2)+1); number=numbersMat(j); subject=subjectsForNumbers{j}; - savePath=readEEG(pathToFile,subject,number,windowEMG,windowEEG,shift,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels); - [meanAccurancysEMG,maxCEMG,cmScaledEMG(:,:)]=svmEciton(savePath,'EMG',k,maxExpC,maxPerClass); - [meanAccurancysEEG,maxCEEG,cmScaledEEG(:,:)]=svmEciton(savePath,'EEG',k,maxExpC,maxPerClass); - [meanAccurancysLF,maxCLF,cmScaledLF(:,:)]=svmEciton(savePath,'LF',k,maxExpC,maxPerClass); - [correlationEMG,maxRidgeParamIndexEMG(:,:)]=ridgeCV(savePath,'EMG',k,ridgeParams); - [correlationEEG,maxRidgeParamIndexEEG(:,:)]=ridgeCV(savePath,'EEG',k,ridgeParams); - [correlationLF,maxRidgeParamIndexLF(:,:)]=ridgeCV(savePath,'LF',k,ridgeParams); + savePath=readAll(pathToFile,subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels,noSynergies); + [meanAccurancysEMG,maxCEMG,cmScaledEMG(:,:)]=svmEciton(savePath,'EMG',k,maxExpC,maxPerClass,eegOffset); + [meanAccurancysEEG,maxCEEG,cmScaledEEG(:,:)]=svmEciton(savePath,'EEG',k,maxExpC,maxPerClass,eegOffset); + [meanAccurancysLF,maxCLF,cmScaledLF(:,:)]=svmEciton(savePath,'LF',k,maxExpC,maxPerClass,eegOffset); + [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); fprintf('%s%i finished %s\n',subject,number,datestr(datetime('now'))) - save(strcat(pathToFile,sprintf('../matlabData/%s_call%s%i.mat',datestr(datetime('now')),subject,number))); + save(strcat(pathToFile,sprintf('../matlabData/%s_call%s%i-%s.mat',datestr(datetime('now')),subject,number,name))); end delete(poolObj) - - diff --git a/usedMcode/classifyAccordingToEMG.m b/usedMcode/classifyAccordingToEMG.m index 048f540..92b261d 100644 --- a/usedMcode/classifyAccordingToEMG.m +++ b/usedMcode/classifyAccordingToEMG.m @@ -1,5 +1,6 @@ -function [classification]=classifyAccordingToEMG(sizeTrainingData, trainingDataEMG, classes,shift,frequency,threshold,pause) +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 diff --git a/usedMcode/downsample.m b/usedMcode/downsample.m deleted file mode 100644 index 931c84f..0000000 --- a/usedMcode/downsample.m +++ /dev/null @@ -1,4 +0,0 @@ -function [out]=downsample(in,noSamples) - stepSize=max(size(in))/noSamples; - out=in(fix((0:noSamples-1)*stepSize)+1); -end \ No newline at end of file diff --git a/usedMcode/generateTrainingData.m b/usedMcode/generateTrainingData.m index 0d5ca7a..747c807 100644 --- a/usedMcode/generateTrainingData.m +++ b/usedMcode/generateTrainingData.m @@ -1,10 +1,8 @@ -function [trainingDataEEG,trainingDataEEGlf,trainingDataEMG,kinPerSec] = generateTrainingData(signal,kin,windowEMG,windowEEG,shift,params,pburgOrder,minEEGFreq,maxEEGFreq,noLFsamples,EMGChannels) +function [trainingDataEEG,trainingDataEEGlf,trainingDataEMG,kinPerSec] = generateTrainingData(signal,kin,windowEMG,windowEEG,shiftEMG,shiftEEG,params,pburgOrder,minEEGFreq,maxEEGFreq,noLFsamples,EMGChannels) 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]); - tempEEGlf=zeros([32,fix(floor(size(signal,1)/signalWindowEEG-1)*windowEEG/shift+1),noLFsamples]); - tempEMG=zeros(fix([6,fix(size(signal,1)/signalWindowEMG-1)*windowEMG/shift+1])); + tempEEG=zeros([32,fix(floor(size(signal,1)/signalWindowEEG-1)*windowEEG/shiftEEG+1),maxEEGFreq-minEEGFreq+1]); + 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'); @@ -16,22 +14,21 @@ [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(E),double(F),filtfilt(double(C),double(D),filtfilt(double(A),double(B),double(signal(:,i))))),bci_sf,windowEEG,shift,pburgOrder,minEEGFreq,maxEEGFreq); - tempEEGlf(i,:,:)=shiftingDownsample(filtfilt(double(X),double(Y),filtfilt(double(V),double(W),double(signal(:,i)))),bci_sf,windowEEG,shift,noLFsamples); + tempEEG(i,:,:)=shiftingPburg(filtfilt(double(E),double(F),filtfilt(double(C),double(D),filtfilt(double(A),double(B),double(signal(:,i))))),bci_sf,windowEEG,shiftEEG,pburgOrder,minEEGFreq,maxEEGFreq); + 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,shift); + 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(trainingDataEEG,1)*1000*shift+windowEEG*1000; - kinPerSec=shiftingKin(kin,windowEEG,shift); + kin(end,1)=size(trainingDataEMG,1)*1000*shiftEMG+windowEMG*1000; + kinPerSec=shiftingKin(kin,windowEMG,shiftEMG); end diff --git a/usedMcode/kfoldCV.m b/usedMcode/kfoldCV.m index af80d84..c1ea267 100644 --- a/usedMcode/kfoldCV.m +++ b/usedMcode/kfoldCV.m @@ -1,32 +1,34 @@ function [model,maxC] = kfoldCV(classification,trainingData,k,cExpMax,maxDataPerClass) noClasses=size(unique(classification),1); -if cExpMax>0 %CV only if necessary - cvAccurancy=zeros(2*cExpMax+1,1); - for cExp=1:2*cExpMax+1 - c=2^(cExp-cExpMax-1); + if cExpMax>0 %CV only if necessary + cvAccurancy=zeros(2*cExpMax+1,1); + for cExp=1:2*cExpMax+1 + c=2^(cExp-cExpMax-1); - randomMapping=transpose(randperm(size(trainingData,1))); - accurancy=zeros(k,3); + randomMapping=transpose(randperm(size(trainingData,1))); + accurancy=zeros(k,3); - parfor i=1:k - 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)); - [~, accurancy(i,:), ~]=svmpredict(testClasses, testData(:,:), model,'-q'); + parfor i=1:k + 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)); + [~, accurancy(i,:), ~]=svmpredict(testClasses, testData(:,:), model,'-q'); + end + + cvAccurancy(cExp)=mean(accurancy(:,1)); end - - cvAccurancy(cExp)=mean(accurancy(:,1)); + [~,maxC]=max(cvAccurancy); + else + maxC=1; %no gridsearch since only one C possible end - [~,maxC]=max(cvAccurancy); -else - maxC=1; %no gridsearch since only one C possible -end bestC=2^(maxC-cExpMax-1); - [balancedClasses, balanceData]=balanceClasses(classification,trainingData,maxDataPerClass,noClasses); + [balancedClasses, balanceData]=balanceClasses(classification,trainingData,... + maxDataPerClass,noClasses); model=svmtrain(balancedClasses, balanceData(:,:),sprintf('-t 0 -c %f -q',bestC)); end diff --git a/usedMcode/myDownsample.m b/usedMcode/myDownsample.m new file mode 100644 index 0000000..b626167 --- /dev/null +++ b/usedMcode/myDownsample.m @@ -0,0 +1,4 @@ +function [out]=myDownsample(in,noSamples) + stepSize=max(size(in))/noSamples; + out=in(fix((0:noSamples-1)*stepSize)+1); +end \ No newline at end of file diff --git a/usedMcode/noOfSynergies.m b/usedMcode/noOfSynergies.m index 16cd0fa..933f4b5 100644 --- a/usedMcode/noOfSynergies.m +++ b/usedMcode/noOfSynergies.m @@ -41,9 +41,10 @@ r2Autoenc=zeros([maxSize,size(trainingData,2)]); +ae=cell([maxSize,1]); parfor i=1:maxSize - ae=trainAutoencoder(trainingData',i,'ShowProgressWindow',false); - predicted=predict(ae,trainingData'); + ae{i}=trainAutoencoder(trainingData',i,'ShowProgressWindow',false); + predicted=predict(ae{i},trainingData'); r2Autoenc(i,:)=correlation2(trainingData,predicted'); end diff --git a/usedMcode/readAll.m b/usedMcode/readAll.m index c382509..058d9f0 100644 --- a/usedMcode/readAll.m +++ b/usedMcode/readAll.m @@ -1,8 +1,11 @@ -function [savePath]=readAll(pathToFile, subject,number,windowEMG,windowEEG,shift,maxFile,threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels) +function [savePath]=readAll(pathToFile, subject,number,windowEMG,windowEEG,shiftEMG,shiftEEG,maxFile,... + threshold,pburgOrder,minEEGFreq,maxEEGFreq,pause,noLFsamples,EMGChannels,noSynergies) %fprintf('start read %s%i %s\n',subject,number,datestr(datetime('now'))); - savePath=strcat(pathToFile,sprintf('../matlabData/%s%i%imsWindowEMG%isWindowEEG%imsShiftFreq%ito%iPause%ipBurg%i.mat',subject,number,windowEMG*1000,windowEEG,shift*1000,minEEGFreq,maxEEGFreq,pause,pburgOrder)); - + savePath=strcat(pathToFile,... + sprintf('../matlabData/%s%i%imsWindowEMG%isWindowEEG%imsShiftEMG%imsShiftEEGFreq%ito%iPause%ipBurg%i.mat',... + subject,number,windowEMG*1000,windowEEG,shiftEMG*1000,shiftEEG*1000,minEEGFreq,maxEEGFreq,pause,pburgOrder)); + %only create file if it doesn't exist yet if exist(savePath, 'file') ~= 2 fprintf(strcat(savePath,' not existing; creating\n')); @@ -13,9 +16,15 @@ kin=cell([maxFile,1]); 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,shift,params,pburgOrder,minEEGFreq,maxEEGFreq,noLFsamples,EMGChannels); + [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 @@ -30,26 +39,49 @@ kinMat=cell2mat(kin); clear trainingDataEEGcell trainingDataEMGcell classesCell kin - classificationWithPause=classifyAccordingToEMG(size(trainingDataEEG,1), trainingDataEMG,classesMat,shift,bci_sf,threshold,pause); + classificationWithPause=classifyAccordingToEMG(trainingDataEMG,... + classesMat,shiftEMG,bci_sf,threshold,pause); clear classesMat - smoothClassification=zeros(size(classificationWithPause)); + smoothClassificationEMG=zeros(size(classificationWithPause)); for i=1:size(classificationWithPause,1) - smoothClassification(i)=round(mode(classificationWithPause(max(i-fix(bci_sf/1000),1):min(i+fix(bci_sf/1000),end)))); + smoothClassificationEMG(i)=round(mode(classificationWithPause(... + max(i-fix(bci_sf/1000),1):min(i+fix(bci_sf/1000),end)))); end clear classificationWithPause - trainingDataEEG=trainingDataEEG(smoothClassification~=-1,:,:); - trainingDataEEGlf=trainingDataEEGlf(smoothClassification~=-1,:,:); - trainingDataEMG=trainingDataEMG(smoothClassification~=-1,:); - classification=smoothClassification(smoothClassification~=-1); - kinematics=kinMat(smoothClassification~=-1,:); + 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 + + 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:4)'; + [~,H]=nnmf(trainingDataEMG,noSynergies); + synergiesNNMF=trainingDataEMG*H'; clear smoothClassification i - save(savePath,'trainingDataEEG','trainingDataEEGlf','trainingDataEMG','classification','kinematics','-v7.3'); + 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/usedMcode/ridgeCV.m b/usedMcode/ridgeCV.m index d1fccca..a4a1ff4 100644 --- a/usedMcode/ridgeCV.m +++ b/usedMcode/ridgeCV.m @@ -1,38 +1,61 @@ -function [correlation,maxRidgeParamIndex]=ridgeCV(savePath,EEG,k,ridgeParams) +function [correlation,maxRidgeParamIndex]=ridgeCV(savePath,data,goal,k,ridgeParams,eegOffset) load(savePath); clear classification; - - if strcmp(EEG,'EEG') - trainingData=trainingDataEEG; - elseif strcmp(EEG,'EMG') + eeg=false; + if strcmp(data,'EEG') + eeg=true; + trainingData=trainingDataEEG(eegOffset+1:end,:); + elseif strcmp(data,'EMG') trainingData=trainingDataEMG; - elseif strcmp(EEG,'LF') - trainingData=trainingDataEEGlf; + elseif strcmp(data,'LF') + eeg=true; + trainingData=trainingDataEEGlf(eegOffset+1:end,:); else - error('only EEG,EMG and LF are valid inputs'); + error('only EEG, EMG and LF are valid inputs for data'); end + factor=size(trainingDataEMG,1)/(size(trainingData,1)+eegOffset); + clear trainingDataEEG; + clear trainingDataEEGlf; clear trainingDataEMG; - - correlation=zeros(size(kinematics,2),1); - maxRidgeParamIndex=zeros(size(kinematics,2),k); - - for j=1:size(kinematics,2) + + 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); + else + error('only kin, Autoenc, PCA and NNMF are valid inputs for goal'); + end + + clear kinematics; + clear synergies*; + + if eeg + predicted=predicted(1:end-eegOffset,:); + end + + correlation=zeros(size(predicted,2),1); + maxRidgeParamIndex=zeros(size(predicted,2),k); + + for j=1:size(predicted,2) randMap=randperm(size(trainingData,1)); - kin=kinematics(:,j); + pred=predicted(:,j); correlations=zeros([k,1]); for i=1:k leaveData=trainingData(mod(randMap,k)==i-1,:); - leaveKin=kin(mod(randMap,k)==i-1); + leavePerd=pred(mod(randMap,k)==i-1); remainingData=trainingData(mod(randMap,k)~=i-1,:); - remainingKin=kin(mod(randMap,k)~=i-1); - %fprintf('%s create %ith model\n',datestr(datetime('now')),i) + remainingPred=pred(mod(randMap,k)~=i-1); - [coeffs,maxRidgeParamIndex(j,i)]=kFoldRidge(remainingData,remainingKin,k,ridgeParams); + [coeffs,maxRidgeParamIndex(j,i)]=kFoldRidge(remainingData,remainingPred,k,ridgeParams); - correlations(i)=ridgeCorrelation(leaveData,leaveKin,coeffs); + correlations(i)=ridgeCorrelation(leaveData,leavePerd,coeffs); end correlation(j)=mean(correlations); diff --git a/usedMcode/shiftingDownsample.m b/usedMcode/shiftingDownsample.m index 79850ad..22da433 100644 --- a/usedMcode/shiftingDownsample.m +++ b/usedMcode/shiftingDownsample.m @@ -5,6 +5,6 @@ noOfWindows=fix(floor(size(signal,1)/signalWindow-1)*shiftProp+1); out=zeros([noOfWindows,noSamples]); for j=1:noOfWindows - out(j,:)=downsample(signal((j-1)*signalShift+1:(j-1)*signalShift+signalWindow),noSamples); + out(j,:)=myDownsample(signal((j-1)*signalShift+1:(j-1)*signalShift+signalWindow),noSamples); end end \ No newline at end of file diff --git a/usedMcode/shiftingKin.m b/usedMcode/shiftingKin.m index f549076..d997a8d 100644 --- a/usedMcode/shiftingKin.m +++ b/usedMcode/shiftingKin.m @@ -1,10 +1,10 @@ -function [kinPerSec]=shiftingKin(kin, windowEEG, shift) - kinPerSec=zeros(fix((max(kin(:,1))-windowEEG*1000)/(shift*1000)),3); +function [kinPerSec]=shiftingKin(kin, windowEMG, shiftEMG) + 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)*shift*1000 & kin(:,1)<=(j-1)*shift*1000+windowEEG*1000,2:4))); + 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)*shift*1000 & kin(:,1)<=(i-1)*shift*1000+windowEEG*1000,2:4))); + 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; diff --git a/usedMcode/shiftingMean.m b/usedMcode/shiftingMean.m new file mode 100644 index 0000000..7c417ae --- /dev/null +++ b/usedMcode/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(ceil(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/usedMcode/svmEciton.m b/usedMcode/svmEciton.m index 857d546..4b98740 100644 --- a/usedMcode/svmEciton.m +++ b/usedMcode/svmEciton.m @@ -1,22 +1,27 @@ -function [meanAccurancy,maxC, cmScaled]= svmEciton(savePath,EEG,k,maxExpC,maxPerClass) +function [meanAccurancy,maxC, cmScaled]= svmEciton(savePath,EEG,k,maxExpC,maxPerClass,eegOffset) 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; + trainingData=trainingDataEEG(1:end-eegOffset,:); + classification=classificationEEG(eegOffset+1:end); elseif strcmp(EEG,'EMG') trainingData=trainingDataEMG; + classification=classificationEMG; elseif strcmp(EEG,'LF') - trainingData=trainingDataEEGlf; + trainingData=trainingDataEEGlf(1:end-eegOffset,:); + classification=classificationEEG(eegOffset+1:end); else - error('only EEG,EMG and LF are valid inputs'); + error('only EEG, EMG and LF are valid inputs'); end - clear trainingDataEEG; - clear trainingDataEEGlf; - clear trainingDataEMG; + clear trainingDataE*; + clear classifiactionE*; + clear kinematics; + clear synergies*; + accurancy=zeros(k,3); maxC=zeros(k,1); @@ -24,13 +29,11 @@ cm=zeros(noClasses); randMap=randperm(size(trainingData,1)); %disp('startCV') - classes=classification; %necessary since otherwise classes are not passed to workers - clear classification; parfor i=1:k leaveOut=trainingData(mod(randMap,k)==i-1,:,:); - leaveClasses=classes(mod(randMap,k)==i-1); + leaveClasses=classification(mod(randMap,k)==i-1); remaining=trainingData(mod(randMap,k)~=i-1,:,:); - remainingClasses=classes(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); @@ -44,9 +47,4 @@ for i=1:size(cm,1) cmScaled(i,:)=cm(i,:)/sum(cm(i,:)); end -% fig=figure; -% imagesc(cmScaled) -% colorbar(); - %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'); end diff --git a/usedMcode/waveformLengthAll.m b/usedMcode/waveformLengthAll.m index 3697b40..34eadcc 100644 --- a/usedMcode/waveformLengthAll.m +++ b/usedMcode/waveformLengthAll.m @@ -1,4 +1,4 @@ -function [EMG]=waveformLengthAll(sig,bci_sf,windowEMG,windowShift) +function [EMG]=waveformLengthAll(sig,bci_sf,windowEMG,shiftEMG) %Filter around 50Hz and below 2 Hz [A,B]= butter(2,[48 52]/(bci_sf/2),'stop'); @@ -6,9 +6,9 @@ [E,F]= butter(2,[148 152]/(bci_sf/2),'stop'); signalWindow=bci_sf*windowEMG; - shiftProp=windowEMG/windowShift; + 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,windowShift); + 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