first draft for first deadline.

toni/octave/functions.m

@@ -1,7 +1,6 @@
 display("functions")
 source("settings.m");
 
-
 function files = getDataFiles(clsName, trainsetPerClass)
 
   global setDataDir;
@@ -35,19 +34,21 @@ function samples = getSamplesForClass(clsName, trainsetPerClass, start, percent)
     raw = {Gyro(start:pEnd,:), Accel(start:pEnd,:), Magnet(start:pEnd,:)}; # create cell array of the wanted samples
     samples{i}.raw = raw; # write them into cell array for each Class
   end
-  display(strcat("training data for '", clsName, "': ", num2str(length(samples)), " samples"));
+  disp(strcat("training data for '", clsName, "': ", num2str(length(samples)), " samples"));
 end
 
 function data = getRawTrainData()
   
   #global trainsetPerClass;
+  global signalStart;
+  global signalEnd;
 
   data = {};
-  data{1}.samples = getSamplesForClass("forwardbend", 9, 10, 0.95);
+  data{1}.samples = getSamplesForClass("forwardbend", 9, signalStart, signalEnd);
-  data{2}.samples = getSamplesForClass("kneebend", 11, 10, 0.95);
+  data{2}.samples = getSamplesForClass("kneebend", 11, signalStart, signalEnd);
-  data{3}.samples = getSamplesForClass("pushups", 11, 10, 0.95);
+  data{3}.samples = getSamplesForClass("pushups", 11, signalStart, signalEnd);
-  data{4}.samples = getSamplesForClass("situps", 8, 10, 0.95);
+  data{4}.samples = getSamplesForClass("situps", 8, signalStart, signalEnd);
-  data{5}.samples = getSamplesForClass("jumpingjack", 13, 10, 0.95);
+  data{5}.samples = getSamplesForClass("jumpingjack", 13, signalStart, signalEnd);
 end
 
 function plotData(data,outPath)
@@ -107,6 +108,12 @@ function windowedData = windowData(data)
   global setWindowSize;
   global setWindowSliding;
   windowedData = {};
+  
+  global setMagnet;
+  global setAccel;
+  global setGyro;
+  global useSignalPCA;
+  global useSignalMG;
 
   for k = 1:numel(data)
     for j = 1:numel(data{k}.samples)
@@ -135,60 +142,68 @@ function windowedData = windowData(data)
       
       for i = winBuffer:setWindowSliding:winEnd-winBuffer #steps for sliding/overlapping windows 
         
-        #accel
+        if setAccel == true
-        winAccel = window(data{k}.samples{j}.raw{1}, i);
+          #accel
-        winsAccelX = [winsAccelX winAccel(:,1)];
+          winAccel = window(data{k}.samples{j}.raw{1}, i);
-        winsAccelY = [winsAccelY winAccel(:,2)];
+          winsAccelX = [winsAccelX winAccel(:,1)];
-        winsAccelZ = [winsAccelZ winAccel(:,3)];
+          winsAccelY = [winsAccelY winAccel(:,2)];
+          winsAccelZ = [winsAccelZ winAccel(:,3)];
+          
+          if useSignalPCA == true
+            [coeff1, score1] = princomp(winAccel');
+            winsAccelPCA = [winsAccelPCA score1(:,1)];  #choose the first axis (eigvec with the biggest eigvalue)
+          endif
+          
+          if useSignalMG == true
+            winAccelMG = getMagnitude(winAccel);
+            winsAccelMG = [winsAccelMG winAccelMG];
+          endif
+          
+        endif
         
-        [coeff1, score1] = princomp(winAccel);
+        if setGyro == true
-        winsAccelPCA = [winsAccelPCA score1(:,1)];  #choose the first axis (eigvec with the biggest eigvalue)
+          #gyro
+          winGyro = window(data{k}.samples{j}.raw{2}, i);
+          winsGyroX = [winsGyroX winGyro(:,1)];
+          winsGyroY = [winsGyroY winGyro(:,2)];
+          winsGyroZ = [winsGyroZ winGyro(:,3)];
+          
+          if useSignalPCA == true
+            [coeff2, score2] = princomp(winGyro');
+            winsGyroPCA = [winsGyroPCA score2(:,1)];
+          endif
+          
+          if useSignalMG == true
+            winGyroMG = getMagnitude(winGyro);
+            winsGyroMG = [winsGyroMG winGyroMG];
+          endif
+            
+        endif
         
-        winAccelMG = getMagnitude(winAccel);
+        if setMagnet == true  
-        winsAccelMG = [winsAccelMG winAccelMG];
+          #magnet
-        
+          winMagnet = window(data{k}.samples{j}.raw{3}, i);
-        #gyro
+          winsMagnetX = [winsMagnetX winMagnet(:,1)];
-        winGyro = window(data{k}.samples{j}.raw{2}, i);
+          winsMagnetY = [winsMagnetY winMagnet(:,2)];
-        winsGyroX = [winsGyroX winGyro(:,1)];
+          winsMagnetZ = [winsMagnetZ winMagnet(:,3)];        
-        winsGyroY = [winsGyroY winGyro(:,2)];
+          
-        winsGyroZ = [winsGyroZ winGyro(:,3)];
+          if useSignalPCA == true
-        
+            [coeff3, score3] = princomp(winMagnet');
-        [coeff2, score2] = princomp(winGyro);
+            winsMagnetPCA = [winsMagnetPCA score3(:,1)];
-        winsGyroPCA = [winsGyroPCA score2(:,1)];
+          endif
-        
+          
-        winGyroMG = getMagnitude(winGyro);
+          if useSignalMG == true
-        winsGyroMG = [winsGyroMG winGyroMG];
+            winMagnetMG = getMagnitude(winMagnet);
-        
+            winsMagnetMG = [winsMagnetMG winMagnetMG];
-        #magnet
+          endif
-        winMagnet = window(data{k}.samples{j}.raw{3}, i);
+          
-        winsMagnetX = [winsMagnetX winMagnet(:,1)];
+        endif
-        winsMagnetY = [winsMagnetY winMagnet(:,2)];
-        winsMagnetZ = [winsMagnetZ winMagnet(:,3)];        
-        
-        [coeff3, score3] = princomp(winMagnet);
-        winsMagnetPCA = [winsMagnetPCA score3(:,1)];
-        
-        winMagnetMG = getMagnitude(winMagnet);
-        winsMagnetMG = [winsMagnetMG winMagnetMG];
                 
       end
-      
+       
-      #write back data
+       #write back data      
-      windowedData{k}.samples{j}.raw{1}.wins = winsAccelX; #X
+      windowedData{k}.samples{j}.raw = [winsAccelX; winsAccelY; winsAccelZ; winsGyroX; winsGyroY; winsGyroZ; winsMagnetX; winsMagnetY; winsMagnetZ; winsAccelPCA; winsGyroPCA; winsMagnetPCA; winsAccelMG; winsGyroMG; winsMagnetMG];
-      windowedData{k}.samples{j}.raw{2}.wins = winsAccelY; #Y
+   
-      windowedData{k}.samples{j}.raw{3}.wins = winsAccelZ; #Z
-      windowedData{k}.samples{j}.raw{4}.wins = winsGyroX;
-      windowedData{k}.samples{j}.raw{5}.wins = winsGyroY;
-      windowedData{k}.samples{j}.raw{6}.wins = winsGyroZ;
-      windowedData{k}.samples{j}.raw{7}.wins = winsMagnetX;
-      windowedData{k}.samples{j}.raw{8}.wins = winsMagnetY;
-      windowedData{k}.samples{j}.raw{9}.wins = winsMagnetZ; 
-      windowedData{k}.samples{j}.raw{10}.wins = winsAccelPCA;
-      windowedData{k}.samples{j}.raw{11}.wins = winsGyroPCA;
-      windowedData{k}.samples{j}.raw{12}.wins = winsMagnetPCA;          
-      windowedData{k}.samples{j}.raw{13}.wins = winsAccelMG;
-      windowedData{k}.samples{j}.raw{14}.wins = winsGyroMG;
-      windowedData{k}.samples{j}.raw{15}.wins = winsMagnetMG;         
     end
   end  
 end
@@ -198,44 +213,66 @@ function features = featureCalculation(data)
   global setAutocorrelationBinSize;
   global setPSDBinSize;
   features = [];
+  
+  global setAutoCorr;
+  global setRMS;
+  global setPSD;
+  global setStatistic;
+  global setPCA;
 
   for k = 1:numel(data)
     for j = 1:numel(data{k}.samples)
-      for i = 1:numel(data{k}.samples{j}.raw)
+      for i = 1:size(data{k}.samples{j}.raw, 1)
-        for m = 1:numel(data{k}.samples{j}.raw{i}.wins)
+        for m = 1:size(data{k}.samples{j}.raw, 2)
-          currentWindow = data{k}.samples{j}.raw{i}.wins{m};
+          currentWindow = data{k}.samples{j}.raw{i, m};
           
-          #autocorrelation on window. split into 5 evenly spaced bins (frequencies are evenly spaced, not number of values ;) ) and calculate mean of bin. 
+          if setAutoCorr == true
-          [autoCorr] = xcorr(currentWindow);
+            #autocorrelation on window. split into 5 evenly spaced bins (amplitudes are evenly spaced, not number of values ;) ) and calculate mean of bin. 
-          [binNum, binCenter] = hist(autoCorr, setAutocorrelationBinSize); #define bins for the data.
+            [autoCorr] = xcorr(currentWindow);
-          binSize = abs(binCenter(end-1) - binCenter(end));
+            [binNum, binCenter] = hist(autoCorr, setAutocorrelationBinSize); #define bins for the data.
-          binEdges = linspace(binCenter(1)-(binSize/2), binCenter(end)+(binSize/2), setAutocorrelationBinSize+1);
+            binSize = abs(binCenter(end-1) - binCenter(end));
-          [binNumc, binIdx] = histc(autoCorr, binEdges);
+            binEdges = linspace(binCenter(1)-(binSize/2), binCenter(end)+(binSize/2), setAutocorrelationBinSize+1);
-          binMeans = getBinMean(autoCorr, binIdx, setAutocorrelationBinSize);
+            [binNumc, binIdx] = histc(autoCorr, binEdges);
+            binMeans = getBinMean(autoCorr, binIdx, setAutocorrelationBinSize);
+            
+          endif
           
-          #calculate the root-mean-square (RMS) of the signal
+          if setRMS == true           
-          rms = sqrt(mean(currentWindow.^2));
+            #calculate the root-mean-square (RMS) of the signal
+            rms = sqrt(mean(currentWindow.^2));
+          endif
+            
+          if setPSD == true
+            #power bands 0.5 to 25hz (useful if the windows are greater then 4s and window sizes to 256, 512..) 
+            [powerBand, w] = periodogram(currentWindow); #fills up fft with zeros            
+            powerEdges = logspace(log10(0.5), log10(25), setPSDBinSize + 2);  #logarithmic bin spaces for 10 bins
+            triFilter = getTriangularFunction(powerEdges, length(powerBand)*2 - 2);    
+            for l = 1:numel(triFilter)
+              filteredBand = triFilter{l} .* powerBand;
+              psd(l) = sum(filteredBand); #sum freq (no log and no dct)
+            end
+          endif
+           
+          if setStatistic == true
+            #statistical features
+            windowMean = mean(currentWindow);
+            windowSTD = std(currentWindow); #standard deviation
+            windowVariance = var(currentWindow);
+            windowKurtosis = kurtosis(currentWindow); #(ger. Wölbung)
+            windowIQR = iqr(currentWindow); #interquartile range
+          endif
           
-          #power bands 0.5 to 25hz (useful if the windows are greater then 4s and window sizes to 256, 512..) 
+          pcaFeature = [];
-          [powerBand, w] = periodogram(currentWindow); #fills up fft with zeros            
+          if setPCA == true 
-          powerEdges = logspace(log10(0.5), log10(25), setPSDBinSize + 2);  #logarithmic bin spaces for 10 bins
+            #pca on windows
-          triFilter = getTriangularFunction(powerEdges, length(powerBand)*2 - 2);    
+            [coeff, score] = princomp(currentWindow');
-          for l = 1:numel(triFilter)
+            pcaFeature = currentWindow' * coeff(:, 1:10); #10 pca feature           
-            filteredBand = triFilter{l} .* powerBand;
+          endif
-            psd(l) = sum(filteredBand); #sum freq (no log and no dct)
-          end
-          
-          #statistical features
-          windowMean = mean(currentWindow);
-          windowSTD = std(currentWindow); #standard deviation
-          windowVariance = var(currentWindow);
-          windowKurtosis = kurtosis(currentWindow); #(ger. Wölbung)
-          windowIQR = iqr(currentWindow); #interquartile range
           
           #put everything together
           classLabel = k; #what class?
           sampleLabel = j; #what sample?
-          features = [features; sampleLabel, classLabel, binMeans, rms, psd, windowMean, windowSTD, windowVariance, windowKurtosis, windowIQR];
+          features = [features; sampleLabel, classLabel, binMeans, rms, psd, windowMean, windowSTD, windowVariance, windowKurtosis, windowIQR, pcaFeature];
         end
       end 
     end
@@ -258,7 +295,7 @@ function value = getBinMean(data, idx, numBins)
     if length(binMembers) == 0
       value(i) = 0;
     else
-      value(i) = mean(binMembers);
+      value(i) = sum(binMembers);
     endif
   end  
 end

toni/octave/genfeatures.m

@@ -13,6 +13,8 @@ classes[1 to 5]
 access the cells using classes{u}.samples{v}.raw{w}
 #}
 source("functions.m");
+source("settings.m")
+global setWindowSize;
 
 classes = {};
 classes = getRawTrainData();
@@ -45,12 +47,12 @@ classes[1 to 5]
       3x WindowsMagnitude (Accel, Gyro, Magnet)
         Win, Win, Win, Win ... <--- single matrices  
         
-access the cells using classes{u}.samples{v}.raw{w}.wins{}
+access the cells using classes{u}.samples{v}.raw{signal, window)
 
 pca uses the eigenvector with the heighest eigenvalue as axis and projects the signals onto it for each sensor. 
 magnitude is calculated using sqrt(x^2 + y^2 + z^2) for each sensor. 
 #}
-windowedClasses = windowData(filteredClasses);
+windowedClasses = windowData(classes);
 
 #calculated features for the 5 signales (x, y, z, MG, PCA) of a sensor
 #{
@@ -60,7 +62,8 @@ data structure of features
 features = featureCalculation(windowedClasses);
 
 #save features
-save features.txt features;
+name = strcat("features_", num2str(setWindowSize),"_xyz_nomag_psd18.txt");
+save(name, 'features');
 display("saved features into features.txt");
 
 

toni/octave/settings.m

@@ -1,8 +1,22 @@
 #global trainsetPerClass = 6; #number of used trainsets for one class
 global setDataDir = "/home/toni/Documents/handygames/HandyGames/daten/";
-global setWindowSize = 256; #in samples per window. even integer! 
+global setWindowSize = 512; #in samples per window. even integer! 
-global setWindowSliding = 320; #in ms - (sampling rate is 10 ms.. so numSamples*10)
+global setWindowSliding = 160; #in ms - (sampling rate is 10 ms.. so numSamples*10)
 global setAutocorrelationBinSize = 5;
-global setPSDBinSize = 10; 
+global setPSDBinSize = 18; 
 global samplerateHZ = 100;
 
+global setMagnet = false;
+global setAccel = true;
+global setGyro = true;
+global useSignalPCA = false;
+global useSignalMG = false;
+
+global signalStart = 1510;
+global signalEnd = 0.8;
+
+global setAutoCorr = true;
+global setRMS = true;
+global setPSD = true;
+global setStatistic = true;
+global setPCA = false; #not implementend / buggy

toni/octave/todo.txt

@@ -23,7 +23,7 @@ autocorrelation bins. jedes signal window wird in 5 frequenz bins unterteilt und
 
 energy features:
 RMS
-power spectrum (|fft|^2). features über logarithmisch überlappende vierrecksfilter, da dreieicksfilter hier keinen sinn ergeben, oder? Wie bei Audio noch eine DCT auf die Ergebnisse des Dreiecksfilters um unwichtige rauszuwerfen?
+power spectrum (|fft|^2). features über logarithmisch überlappende dreiecksfilter Wie bei Audio noch eine DCT auf die Ergebnisse des Dreiecksfilters um unwichtige rauszuwerfen?
 
 statistical features:
 mean
@@ -47,10 +47,13 @@ Training:
 one leave out.
 
 
-
+Prasi:
-
+
-
+pure klassifikation. ohne zu schecken was macht das fenster vorher.
-
+
+svm aufwand?
+
+autokorrelation: ist ein signal eigentlich periodisch, wenn ich paar lags einbaue? peak in der mitte sagt ja. wir sind jetzt an der amplitude interessiert. da die meisten übungen wohl periodisch sind ;).
 
 
 

toni/octave/training.m

@@ -1,9 +1,12 @@
 #train features using svm
 display("Train Features")
+page_screen_output(0);
+page_output_immediately(1);
+more off;
 
 #load all features
 # features = [sampleLabel, classLabel, binMeans, rms, psd, windowMean, windowSTD, windowVariance, windowKurtosis, windowIQR];
-load "features.txt"; #matrix is also called features
+load "eval/512/features_512_xyz_nomag_psd18.txt"; #matrix is also called features
 
 # split features into training and test features using leave-one-out method
 # class idx:
@@ -13,37 +16,79 @@ load "features.txt"; #matrix is also called features
 # idx 4 -> situps
 # idx 5 -> jumpingjack
 
-# define which sampleSet is used as testset and not for training.
+#remove features for evaluation
-leaveOut = find(features(:,1) == 3 & features(:,2) == 2); #sampleset 3 class 2
+#features(:, 3:7) = []; #remove binMeans
-testFeatures = features(leaveOut, :); #set testSignal
+#features(:, 8:25) = []; #remove psd
-features(leaveOut,:) = []; #remove the testFeatures
+#features(:, 26) = []; #remove rms
-features(:,1) = []; #remove the sampleLabel
+#features(:, 27:31) = []; #remove statistical stuff
 
-# bring the feature matrix into libsvm format. 
+size(features, 2)
-# 1. the label vector: 
-trainLabel = features(:,1);
 
-# 2. sparse matrix with every feature in one column:
+#samples_class = features(:,1:2);
-features(:,1) = []; #remove the classLabel
+#for numClass = 1 : max(samples_class(:,2)) 
-trainFeatures = sparse(features);
+#  for numSample = 1 : max(samples_class(find(samples_class(:,2)==1),1))
+    
+    t=cputime;
 
-# before training we need to scale the feature values
+    # define which sampleSet is used as testset and not for training.    
-minimums = min(trainFeatures);
+    leaveOut = randperm(length(features), 1000);
-ranges = max(trainFeatures) - minimums;
+    
+    numSample = 2;
+    numClass = 5;
+    #leaveOut = find(features(:,1) == numSample & features(:,2) == numClass); #sampleset x class y
+    testFeatures = features(leaveOut, :); #set testSignal
+    features(leaveOut,:) = []; #remove the testFeatures
+    features(:,1) = []; #remove the sampleLabel 
+    
 
-trainFeatures = (trainFeatures - repmat(minimums, size(trainFeatures, 1), 1)) ./ repmat(ranges, size(trainFeatures, 1), 1); 
+    # bring the feature matrix into libsvm format. 
+    # 1. the label vector: 
+    trainLabel = features(:,1);
 
-# training: svm with default settings
+    # 2. sparse matrix with every feature in one column:
-model = svmtrain(trainLabel, trainFeatures);
+    features(:,1) = []; #remove the classLabel
+    trainFeatures = sparse(features);
 
-display("Classify Features")
+    #write out libsvm file
-# for testing we need to scale again
+    #libsvmwrite(strcat("trainFeatures_512_", num2str(numClass), "_", num2str(numSample), ".txt"), trainLabel, trainFeatures);
-testLabel = testFeatures(:,2);
-testFeatures(:,1:2) = []; #remove the labels
-testFeatures = (testFeatures - repmat(minimums, size(testFeatures, 1), 1)) ./ repmat(ranges, size(testFeatures, 1), 1);
-
-# classification
-[predict_label, accuracy, dec_values] = svmpredict(testLabel, testFeatures, model);
 
+    # before training we need to scale the feature values
+    minimums = min(trainFeatures);
+    ranges = max(trainFeatures) - minimums;
+
+    trainFeatures = (trainFeatures - repmat(minimums, size(trainFeatures, 1), 1)) ./ repmat(ranges, size(trainFeatures, 1), 1); 
+
+    # training: svm with default settings
+    model = svmtrain(trainLabel, trainFeatures, '-h 0 -c 32768 -g 8 -q');
+
+    disp("Classify Features");
+    # for testing we need to scale again
+    testLabel = testFeatures(:,2);
+    testFeatures(:,1:2) = []; #remove the labels
+    testFeatures = sparse(testFeatures);
+
+    #write out libsvm file
+    #libsvmwrite(strcat("testFeatures_512_", num2str(numClass), "_", num2str(numSample), ".txt"), testLabel, testFeatures);
+
+    testFeatures = (testFeatures - repmat(minimums, size(testFeatures, 1), 1)) ./ repmat(ranges, size(testFeatures, 1), 1);
+
+    # classification
+    [predict_label, accuracy, dec_values] = svmpredict(testLabel, testFeatures, model);
+    
+       
+    # get evaluation matrix
+    evaluation = [];
+    for i = 1:5
+      for j = 1:5          
+        evaluation(i,j) = sum(testLabel == i & predict_label == j);
+      end
+    end
+    disp(evaluation);
+    
+    save evaluation_512.txt evaluation accuracy;
+       
+    printf('Total cpu time: %f seconds\n', cputime-t);
+#  end
+#end