From 12c5cae2535a7d1f897a4c10485e2a2b942107d9 Mon Sep 17 00:00:00 2001
From: toni <toni@localhost.localdomain>
Date: Mon, 27 Nov 2017 21:17:30 +0100
Subject: [PATCH] ref #6 - autocorrelation implemented

---
 .../conductorswatch/BpmEstimator.java         |   4 +-
 .../conductorswatch/WorkerFragment.java       |   6 -
 .../conductorswatch/utilities/Utils.java      | 142 ++++++++++++++----
 android/ConductorsWatch/build.gradle          |   2 +-
 octave/AutoCorrMethodNew_Watch.m              |   6 +-
 5 files changed, 120 insertions(+), 40 deletions(-)

diff --git a/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/BpmEstimator.java b/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/BpmEstimator.java
index 4a9f6ec..658b0fe 100644
--- a/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/BpmEstimator.java
+++ b/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/BpmEstimator.java
@@ -6,7 +6,9 @@ import android.hardware.SensorEvent;
 import android.hardware.SensorEventListener;
 import android.hardware.SensorManager;
 import android.os.Handler;
+import android.util.Log;
 
+import java.util.Arrays;
 import java.util.List;
 import java.util.concurrent.CopyOnWriteArrayList;
 import java.util.concurrent.ThreadLocalRandom;
@@ -17,8 +19,6 @@ import de.tonifetzer.conductorswatch.utilities.Utils;
  * Created by toni on 13/11/17.
  */
 
-//TODO: diesen estimator testen. kommen alle messungen? wie sind die messungen zeitlich voneinander verschieden?
-//TODO: klapp das wirklich mit den 4ms. passt der buffer? gehen auch höhere zeiten?
 //TODO: einfügen der logik autoCorr + FindPeaks
 public class BpmEstimator implements SensorEventListener {
 
diff --git a/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/WorkerFragment.java b/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/WorkerFragment.java
index 9a0cae6..34f267a 100644
--- a/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/WorkerFragment.java
+++ b/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/WorkerFragment.java
@@ -111,12 +111,6 @@ public class WorkerFragment extends Fragment implements Metronome.OnMetronomeLis
 
         // stop the worker thread for bpm estimator
         mBpmEstimator.stop();
-        /*mBpmThread.interrupt();
-        try {
-            mBpmThread.join();
-        } catch (InterruptedException e) {
-            e.printStackTrace();
-        }*/
 
         // stop the worker thread for metronom
         mMetronome.stop();
diff --git a/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/utilities/Utils.java b/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/utilities/Utils.java
index f71b8fa..f4b0909 100644
--- a/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/utilities/Utils.java
+++ b/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/utilities/Utils.java
@@ -3,11 +3,9 @@ package de.tonifetzer.conductorswatch.utilities;
 import android.content.Context;
 import android.content.res.Resources;
 import android.util.DisplayMetrics;
-
-import org.jtransforms.fft.DoubleFFT_1D;
-
+import org.jtransforms.fft.FloatFFT_1D;
 import java.util.ArrayList;
-import java.util.Queue;
+import java.util.Arrays;
 
 
 public class Utils {
@@ -40,18 +38,25 @@ public class Utils {
         }
     }
 
-    //TODO: implement methods providing x,y,z and ts as solo vectors
-    //TODO: implement sliding window counter
     public static class AccelerometerWindowBuffer extends ArrayList<AccelerometerData> {
 
         private int mWindowSize;
         private int mOverlapSize;
         private int mOverlapCounter;
+        private float[] mX;
+        private float[] mY;
+        private float[] mZ;
+        private long[] mTs;
 
         public AccelerometerWindowBuffer(int windowSize, int overlap){
             this.mWindowSize = windowSize;
             this.mOverlapSize = overlap;
             mOverlapCounter = 1;
+
+            mX = new float[this.mWindowSize];
+            mY = new float[this.mWindowSize];
+            mZ = new float[this.mWindowSize];
+            mTs = new long[this.mWindowSize];
         }
 
         public boolean add(AccelerometerData ad){
@@ -60,6 +65,14 @@ public class Utils {
                 removeRange(0, size() - mWindowSize);
             }
 
+            //update the double arrays.
+            for (int i = 0; i < size(); ++i) {
+                mX[i] = get(i).x;
+                mY[i] = get(i).y;
+                mZ[i] = get(i).z;
+                mTs[i] = get(i).ts;
+            }
+
             ++mOverlapCounter;
             return r;
         }
@@ -79,34 +92,107 @@ public class Utils {
         public AccelerometerData getOldest() {
             return get(0);
         }
+
+        public float[] getX(){
+            return mX;
+        }
+
+        public float[] getY(){
+            return mY;
+        }
+
+        public float[] getZ(){
+            return mZ;
+        }
+
+        public long[] getTs(){
+            return mTs;
+        }
     }
 
 
-    public static double sqr(double x) {
+    public static float sqr(float x) {
         return x * x;
     }
 
-    //TODO: implement maxLag as input
-    //TODO: implement positive and negative lag output
-    public void fftAutoCorrelation(double [] x, double [] ac) {
-        int n = x.length;
-        // Assumes n is even.
-        DoubleFFT_1D fft = new DoubleFFT_1D(n);
-        fft.realForward(x);
-        //ac[0] = sqr(x[0]); // For normal xcov
-        ac[0] = 0;  // For statistical convention, zero out the mean
-        ac[1] = sqr(x[1]);
-        for (int i = 2; i < n; i += 2) {
-            ac[i] = sqr(x[i]) + sqr(x[i+1]);
-            ac[i+1] = 0;
-        }
-        DoubleFFT_1D ifft = new DoubleFFT_1D(n);
-        ifft.realInverse(ac, true);
-        //For statistical convention, normalize by dividing through with variance
-        for (int i = 1; i < n; i++){
-            ac[i] /= ac[0];
-        }
-        ac[0] = 1;
+    public static int nextPow2(int a){
+        return a == 0 ? 0 : 32 - Integer.numberOfLeadingZeros(a - 1);
     }
 
+    public static float mean(float[] data){
+        float sum = 0;
+        for (int i = 0; i < data.length; i++) {
+            sum += data[i];
+        }
+        return sum / data.length;
+    }
+
+    public static float[] removeZero(float[] array){
+        int j = 0;
+        for( int i=0;  i<array.length;  i++ )
+        {
+            if (array[i] != 0)
+                array[j++] = array[i];
+        }
+        float[] newArray = new float[j];
+        System.arraycopy( array, 0, newArray, 0, j );
+
+        return newArray;
+    }
+
+    //TODO: errorhandling maxLag = 0;
+    //TODO: größeren Testcase schreiben
+    public static float[] fftAutoCorrelation(float[] data, int maxLag) {
+
+        int n = data.length;
+        float[] x = Arrays.copyOf(data, n);
+        int mxl = Math.min(maxLag, n - 1);
+        int ceilLog2 = nextPow2(2*n -1);
+        int n2 = (int) Math.pow(2,ceilLog2);
+
+        // x - mean(x) (pointwise)
+        float x_mean = mean(x);
+        for(int i = 0; i < x.length; ++i){
+            x[i] -= x_mean;
+        }
+
+        // double the size of x and fill up with zeros. if x is not even, add additional 0
+        float[] x2 = new float[n2 * 2]; //need double the size for fft.realForwardFull (look into method description)
+        Arrays.fill(x2, 0);
+        System.arraycopy(x,0, x2, 0, x.length);
+
+        // x_fft calculate fft 1D
+        FloatFFT_1D fft = new FloatFFT_1D(n2);
+        fft.realForwardFull(x2);
+
+        // Cr = abs(x_fft).^2 (absolute with complex numbers is (r^2) + (i^2)
+        float[] Cr = new float[n2 * 2];
+        int j = 0;
+        for(int i = 0; i < x2.length; ++i){
+            Cr[j++] = sqr(x2[i]) + sqr(x2[i+1]);
+            ++i;  //skip the complex part
+        }
+
+        // ifft(Cr,[],1)
+        FloatFFT_1D ifft = new FloatFFT_1D(n2);
+        ifft.realInverseFull(Cr, true);
+
+        // remove complex part and scale/normalize
+        float[] c1 = new float[n2];
+        j = 0;
+        for(int i = 0; i < Cr.length; ++i){
+            c1[j++] = Cr[i] / Cr[0];
+            ++i;  //skip the complex part
+        }
+
+        // Keep only the lags we want and move negative lags before positive lags.
+        float[] c = new float[(mxl * 2) + 1];
+        System.arraycopy(c1, 0, c, mxl, mxl + 1); // +1 to place the 1.0 in the middle of correlation
+        System.arraycopy(c1, n2 - mxl, c, 0, mxl);
+
+        return c;
+    }
+
+    //TODO: findPeaks
+
 }
diff --git a/android/ConductorsWatch/build.gradle b/android/ConductorsWatch/build.gradle
index 020eae2..e6b32bc 100644
--- a/android/ConductorsWatch/build.gradle
+++ b/android/ConductorsWatch/build.gradle
@@ -7,7 +7,7 @@ buildscript {
         jcenter()
     }
     dependencies {
-        classpath 'com.android.tools.build:gradle:3.0.0'
+        classpath 'com.android.tools.build:gradle:3.0.1'
         
 
         // NOTE: Do not place your application dependencies here; they belong
diff --git a/octave/AutoCorrMethodNew_Watch.m b/octave/AutoCorrMethodNew_Watch.m
index f55a51f..5d06daf 100644
--- a/octave/AutoCorrMethodNew_Watch.m
+++ b/octave/AutoCorrMethodNew_Watch.m
@@ -84,9 +84,9 @@ for i = window_size+1:length(data)
 	if(mod(i,overlap) == 0)
 
 		%measure periodicity of window and use axis with best periodicity
-		[corr_x, lag_x] = xcov(m(i-window_size:i,3), (window_size/4), "coeff");
-		[corr_y, lag_y] = xcov(m(i-window_size:i,4), (window_size/4), "coeff");
-		[corr_z, lag_z] = xcov(m(i-window_size:i,5), (window_size/4), "coeff");
+		[corr_x, lag_x] = xcov(m(i-window_size:i,3), (window_size/2), "coeff");
+		[corr_y, lag_y] = xcov(m(i-window_size:i,4), (window_size/2), "coeff");
+		[corr_z, lag_z] = xcov(m(i-window_size:i,5), (window_size/2), "coeff");
 
 		corr_x_pos = corr_x;
 		corr_y_pos = corr_y;