From 709a846a91e6027f6a27be954fd758318df62ee9 Mon Sep 17 00:00:00 2001 From: toni Date: Sun, 17 Dec 2017 00:35:57 +0100 Subject: [PATCH] ref #6 -refactoring - added interp1 - added peak finder - added bpm estimation - added findbestaxis --- .../conductorswatch/BpmEstimator.java | 2 + java/src/main/java/Main.java | 72 ++--- java/src/main/java/Utils.java | 272 ------------------ .../java/bpmEstimation/AccelerometerData.java | 31 ++ .../AccelerometerInterpolator.java | 114 ++++++++ .../AccelerometerWindowBuffer.java | 89 ++++++ .../java/bpmEstimation/AutoCorrelation.java | 79 +++++ java/src/main/java/bpmEstimation/Peaks.java | 207 +++++++++++++ java/src/main/java/utilities/Utils.java | 232 +++++++++++++++ 9 files changed, 793 insertions(+), 305 deletions(-) delete mode 100644 java/src/main/java/Utils.java create mode 100644 java/src/main/java/bpmEstimation/AccelerometerData.java create mode 100644 java/src/main/java/bpmEstimation/AccelerometerInterpolator.java create mode 100644 java/src/main/java/bpmEstimation/AccelerometerWindowBuffer.java create mode 100644 java/src/main/java/bpmEstimation/AutoCorrelation.java create mode 100644 java/src/main/java/bpmEstimation/Peaks.java create mode 100644 java/src/main/java/utilities/Utils.java 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 df414fd..7f2c2ad 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 @@ -34,6 +34,8 @@ public class BpmEstimator implements SensorEventListener { if(mAccelerometerWindowBuffer.isNextWindowReady()){ + //TODO: calculate average samplerate using the window - interpolation of sensordata would be better i think + //int randomNum = ThreadLocalRandom.current().nextInt(0, 240 + 1); long diff = mAccelerometerWindowBuffer.getYongest().ts - mAccelerometerWindowBuffer.getOldest().ts; diff --git a/java/src/main/java/Main.java b/java/src/main/java/Main.java index a2b9c6f..ef82c27 100644 --- a/java/src/main/java/Main.java +++ b/java/src/main/java/Main.java @@ -1,3 +1,6 @@ +import bpmEstimation.*; +import utilities.Utils; + import java.awt.*; import java.io.BufferedReader; import java.io.File; @@ -24,8 +27,7 @@ public class Main { for (File file : listOfFiles) { if (file.isFile() && file.getName().contains(".csv")) { - Utils.AccelerometerWindowBuffer accWindowBuffer = new Utils.AccelerometerWindowBuffer(4096, 256); - + AccelerometerWindowBuffer accWindowBuffer = new AccelerometerWindowBuffer(4096, 256); //read the file line by line try (BufferedReader br = new BufferedReader(new FileReader(file))) { @@ -36,86 +38,83 @@ public class Main { //if linear acc if(measurement[1].equals("2")){ long ts = Long.parseLong(measurement[0]); - float x = Float.parseFloat(measurement[2]); - float y = Float.parseFloat(measurement[3]); - float z = Float.parseFloat(measurement[4]); - accWindowBuffer.add(new Utils.AccelerometerData(ts, x, y, z)); + double x = Double.parseDouble(measurement[2]); + double y = Double.parseDouble(measurement[3]); + double z = Double.parseDouble(measurement[4]); + accWindowBuffer.add(new AccelerometerData(ts, x, y, z)); } //do calculation stuff if(accWindowBuffer.isNextWindowReady()){ + AccelerometerInterpolator acInterp = new AccelerometerInterpolator(accWindowBuffer, 6); + //print raw x,y,z double[] dTs = IntStream.range(0, accWindowBuffer.getTs().length).mapToDouble(i -> accWindowBuffer.getTs()[i]).toArray(); - double[] dX = IntStream.range(0, accWindowBuffer.getX().length).mapToDouble(i -> accWindowBuffer.getX()[i]).toArray(); - double[] dY = IntStream.range(0, accWindowBuffer.getY().length).mapToDouble(i -> accWindowBuffer.getY()[i]).toArray(); - double[] dZ = IntStream.range(0, accWindowBuffer.getZ().length).mapToDouble(i -> accWindowBuffer.getZ()[i]).toArray(); - + double[] dTsInterp = IntStream.range(0, acInterp.getTs().length).mapToDouble(i -> acInterp.getTs()[i]).toArray(); Plot plotRaw = Plot.plot(Plot.plotOpts(). title("Raw Acc Data"). legend(Plot.LegendFormat.BOTTOM)). - series("x", Plot.data().xy(dTs, dX), Plot.seriesOpts().color(Color.RED)). - series("y", Plot.data().xy(dTs, dY), Plot.seriesOpts().color(Color.BLUE)). - series("z", Plot.data().xy(dTs, dZ), Plot.seriesOpts().color(Color.GREEN)); + series("x", Plot.data().xy(dTsInterp, acInterp.getX()), Plot.seriesOpts().color(Color.RED)). + series("y", Plot.data().xy(dTsInterp, acInterp.getY()), Plot.seriesOpts().color(Color.BLUE)). + series("z", Plot.data().xy(dTsInterp, acInterp.getZ()), Plot.seriesOpts().color(Color.GREEN)); windowRaw.set(plotRaw.draw()); //auto corr - float[] xAutoCorr = Utils.fftAutoCorrelation(accWindowBuffer.getX(), 1024); - float[] yAutoCorr = Utils.fftAutoCorrelation(accWindowBuffer.getY(), 1024); - float[] zAutoCorr = Utils.fftAutoCorrelation(accWindowBuffer.getZ(), 1024); + double[] xAutoCorr = new AutoCorrelation(acInterp.getX(), 1024).getCorr(); + double[] yAutoCorr = new AutoCorrelation(acInterp.getY(), 1024).getCorr(); + double[] zAutoCorr = new AutoCorrelation(acInterp.getZ(), 1024).getCorr(); //print autocorr int[] tmp = IntStream.rangeClosed(-((xAutoCorr.length - 1)/2), ((xAutoCorr.length - 1)/2)).toArray(); double[] rangeAuto = IntStream.range(0, tmp.length).mapToDouble(i -> tmp[i]).toArray(); - double[] dXAuto = IntStream.range(0, xAutoCorr.length).mapToDouble(i -> xAutoCorr[i]).toArray(); - double[] dYAuto = IntStream.range(0, yAutoCorr.length).mapToDouble(i -> yAutoCorr[i]).toArray(); - double[] dZAuto = IntStream.range(0, zAutoCorr.length).mapToDouble(i -> zAutoCorr[i]).toArray(); - Plot plotCorr = Plot.plot(Plot.plotOpts(). title("Auto Correlation"). legend(Plot.LegendFormat.BOTTOM)). - series("x", Plot.data().xy(rangeAuto, dXAuto), Plot.seriesOpts().color(Color.RED)). - series("y", Plot.data().xy(rangeAuto, dYAuto), Plot.seriesOpts().color(Color.BLUE)). - series("z", Plot.data().xy(rangeAuto, dZAuto), Plot.seriesOpts().color(Color.GREEN)); + series("x", Plot.data().xy(rangeAuto, xAutoCorr), Plot.seriesOpts().color(Color.RED)). + series("y", Plot.data().xy(rangeAuto, yAutoCorr), Plot.seriesOpts().color(Color.BLUE)). + series("z", Plot.data().xy(rangeAuto, zAutoCorr), Plot.seriesOpts().color(Color.GREEN)); windowAuto.set(plotCorr.draw()); - //find peaks - LinkedList peaksX = Utils.findPeaks(xAutoCorr, 50, 0.1f, 0, false); + Peaks pX = new Peaks(xAutoCorr, 50, 0.1f, 0, false); + LinkedList peaksX = pX.getPeaksIdx(); double[] dPeaksXX = IntStream.range(0, peaksX.size()).mapToDouble(i -> (peaksX.get(i) - 1024)).toArray();//peaks.stream().mapToDouble(i->i).toArray(); - double[] dPeaksXY = IntStream.range(0, peaksX.size()).mapToDouble(i -> (dXAuto[peaksX.get(i)])).toArray(); + double[] dPeaksXY = IntStream.range(0, peaksX.size()).mapToDouble(i -> (xAutoCorr[peaksX.get(i)])).toArray(); Plot plotPeaksX = Plot.plot(Plot.plotOpts(). title("Peak Detection on X"). legend(Plot.LegendFormat.BOTTOM)). - series("x", Plot.data().xy(rangeAuto, dXAuto), Plot.seriesOpts().color(Color.RED)). + series("x", Plot.data().xy(rangeAuto, xAutoCorr), Plot.seriesOpts().color(Color.RED)). series("Peaks", Plot.data().xy(dPeaksXX, dPeaksXY), Plot.seriesOpts().color(Color.CYAN). marker(Plot.Marker.DIAMOND).line(Plot.Line.NONE)); windowPeaksX.set(plotPeaksX.draw()); - LinkedList peaksY = Utils.findPeaks(yAutoCorr, 50, 0.1f, 0, false); + Peaks pY = new Peaks(yAutoCorr, 50, 0.1f, 0, false); + LinkedList peaksY = pY.getPeaksIdx(); double[] dPeaksYX = IntStream.range(0, peaksY.size()).mapToDouble(i -> (peaksY.get(i) - 1024)).toArray();//peaks.stream().mapToDouble(i->i).toArray(); - double[] dPeaksYY = IntStream.range(0, peaksY.size()).mapToDouble(i -> (dYAuto[peaksY.get(i)])).toArray(); + double[] dPeaksYY = IntStream.range(0, peaksY.size()).mapToDouble(i -> (yAutoCorr[peaksY.get(i)])).toArray(); Plot plotPeaksY = Plot.plot(Plot.plotOpts(). title("Peak Detection on Y"). legend(Plot.LegendFormat.BOTTOM)). - series("x", Plot.data().xy(rangeAuto, dYAuto), Plot.seriesOpts().color(Color.RED)). + series("x", Plot.data().xy(rangeAuto, yAutoCorr), Plot.seriesOpts().color(Color.RED)). series("Peaks", Plot.data().xy(dPeaksYX, dPeaksYY), Plot.seriesOpts().color(Color.CYAN). marker(Plot.Marker.DIAMOND).line(Plot.Line.NONE)); windowPeaksY.set(plotPeaksY.draw()); - LinkedList peaksZ = Utils.findPeaks(zAutoCorr, 50, 0.1f, 0, false); + Peaks pZ = new Peaks(zAutoCorr, 50, 0.1f, 0, false); + LinkedList peaksZ = pZ.getPeaksIdx(); double[] dPeaksZX = IntStream.range(0, peaksZ.size()).mapToDouble(i -> (peaksZ.get(i) - 1024)).toArray();//peaks.stream().mapToDouble(i->i).toArray(); - double[] dPeaksZY = IntStream.range(0, peaksZ.size()).mapToDouble(i -> (dZAuto[peaksZ.get(i)])).toArray(); + double[] dPeaksZY = IntStream.range(0, peaksZ.size()).mapToDouble(i -> (zAutoCorr[peaksZ.get(i)])).toArray(); Plot plotPeaksZ = Plot.plot(Plot.plotOpts(). title("Peak Detection on Z"). legend(Plot.LegendFormat.BOTTOM)). - series("x", Plot.data().xy(rangeAuto, dZAuto), Plot.seriesOpts().color(Color.RED)). + series("x", Plot.data().xy(rangeAuto, zAutoCorr), Plot.seriesOpts().color(Color.RED)). series("Peaks", Plot.data().xy(dPeaksZX, dPeaksZY), Plot.seriesOpts().color(Color.CYAN). marker(Plot.Marker.DIAMOND).line(Plot.Line.NONE)); @@ -124,6 +123,13 @@ public class Main { //fill hols improve peaks //estimate bpm between detected peaks + System.out.println("BPM-X: " + pX.getBPM(6)); + System.out.println("BPM-Y: " + pY.getBPM(6)); + System.out.println("BPM-Z: " + pZ.getBPM(6)); + + //todo: statistikzeuch und mit matlab vergleichen. + //todo: todos machen. lol + //todo: kleiner fenstergrößen testen. so ist doch etwas langsam auf der Uhr. int dummyForBreakpoint = 0; } diff --git a/java/src/main/java/Utils.java b/java/src/main/java/Utils.java deleted file mode 100644 index fb6cbe7..0000000 --- a/java/src/main/java/Utils.java +++ /dev/null @@ -1,272 +0,0 @@ -import org.jtransforms.fft.FloatFFT_1D; - -import javax.swing.*; -import java.awt.*; -import java.awt.image.BufferedImage; -import java.util.ArrayList; -import java.util.Arrays; -import java.util.LinkedList; - - -public class Utils { - public static float getDistance(float x1, float y1, float x2, float y2) { - return (float) Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); - } - - public static class AccelerometerData { - - public float x,y,z; - public long ts; - - public AccelerometerData(long ts, float x, float y, float z){ - this.ts = ts; - this.x = x; - this.y = y; - this.z = z; - } - } - - public static class AccelerometerWindowBuffer extends ArrayList { - - 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){ - boolean r = super.add(ad); - if (size() > mWindowSize){ - 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; - } - - public boolean isNextWindowReady(){ - if(size() == mWindowSize && mOverlapCounter > mOverlapSize){ - mOverlapCounter = 1; - return true; - } - return false; - } - - public AccelerometerData getYongest() { - return get(size() - 1); - } - - 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 float sqr(float x) { - return x * x; - } - - 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 findPeaks(float[] data, int width, float threshold, float decayRate, boolean isRelative) { - LinkedList peaks = new LinkedList<>(); - int maxp = 0; - int mid = 0; - int end = data.length; - float av = data[0]; - while (mid < end) { - av = decayRate * av + (1 - decayRate) * data[mid]; - if (av < data[mid]) - av = data[mid]; - int i = mid - width; - if (i < 0) - i = 0; - int stop = mid + width + 1; - if (stop > data.length) - stop = data.length; - maxp = i; - for (i++; i < stop; i++) - if (data[i] > data[maxp]) - maxp = i; - if (maxp == mid) { - if (overThreshold(data, maxp, width, threshold, isRelative, av)){ - peaks.add(new Integer(maxp)); - } - } - mid++; - } - return peaks; - } - - private static boolean overThreshold(float[] data, int index, int width, - float threshold, boolean isRelative, - float av) { - int pre = 3; - int post = 1; - - if (data[index] < av) - return false; - if (isRelative) { - int iStart = index - pre * width; - if (iStart < 0) - iStart = 0; - int iStop = index + post * width; - if (iStop > data.length) - iStop = data.length; - float sum = 0; - int count = iStop - iStart; - while (iStart < iStop) - sum += data[iStart++]; - return (data[index] > sum / count + threshold); - } else - return (data[index] > threshold); - } - - - @SuppressWarnings("serial") - public static class ShowPNG extends JFrame - { - JLabel mLabel; - ImageIcon mIcon; - - public ShowPNG(){ - mLabel = new JLabel(); - this.setLayout(new GridLayout(1,1)); - this.setSize(800, 640); - this.add(mLabel); - this.setVisible(true); - } - - public void set(BufferedImage bi){ - - mIcon = new ImageIcon(bi); - mLabel.setVisible(false); - mLabel.setIcon(mIcon); - mLabel.setVisible(true); - } - } - -} diff --git a/java/src/main/java/bpmEstimation/AccelerometerData.java b/java/src/main/java/bpmEstimation/AccelerometerData.java new file mode 100644 index 0000000..801c37a --- /dev/null +++ b/java/src/main/java/bpmEstimation/AccelerometerData.java @@ -0,0 +1,31 @@ +package bpmEstimation; + +/** + * Created by toni on 15/12/17. + */ +public class AccelerometerData { + + public double x,y,z; + public long ts; + + public AccelerometerData(long ts, double x, double y, double z){ + this.ts = ts; + this.x = x; + this.y = y; + this.z = z; + } + + @Override + public boolean equals(Object other){ + + if (this == other) + return true; + + if (!(other instanceof AccelerometerData)) { + return false; + } + + AccelerometerData ad = (AccelerometerData) other; + return ts == ad.ts; + } +} \ No newline at end of file diff --git a/java/src/main/java/bpmEstimation/AccelerometerInterpolator.java b/java/src/main/java/bpmEstimation/AccelerometerInterpolator.java new file mode 100644 index 0000000..7806995 --- /dev/null +++ b/java/src/main/java/bpmEstimation/AccelerometerInterpolator.java @@ -0,0 +1,114 @@ +package bpmEstimation; + +import utilities.Utils; + +import java.util.Arrays; + +/** + * Created by toni on 16/12/17. + */ +public class AccelerometerInterpolator { + + private double[] mX; + private double[] mY; + private double[] mZ; + private long[] mTsInterp; + + public AccelerometerInterpolator(AccelerometerWindowBuffer ab, double sampleRate_ms){ + + if(sampleRate_ms <= 0){ + sampleRate_ms = Utils.mean(Utils.diff(ab.getTs())); + } + + long size = (ab.getYongest().ts - (ab.getOldest().ts - (long) sampleRate_ms)) / (long) sampleRate_ms; + mTsInterp = new long[(int)size]; + int j = 0; + for(long i = ab.getOldest().ts; i <= ab.getYongest().ts; i += sampleRate_ms){ + mTsInterp[j++] = i; + } + + mX = interpLinear(ab.getTs(), ab.getX(), mTsInterp); + mY = interpLinear(ab.getTs(), ab.getY(), mTsInterp); + mZ = interpLinear(ab.getTs(), ab.getZ(), mTsInterp); + } + + public long[] getTs(){ + return mTsInterp; + } + + public double[] getX(){ + return mX; + } + + public double[] getY(){ + return mY; + } + + public double[] getZ(){ + return mZ; + } + + public static double[] interpLinear(double[] x, double[] y, double[] xi) throws IllegalArgumentException { + if (x.length != y.length) { + throw new IllegalArgumentException("X and Y must be the same length"); + } + if (x.length == 1) { + throw new IllegalArgumentException("X must contain more than one value"); + } + double[] dx = new double[x.length - 1]; + double[] dy = new double[x.length - 1]; + double[] slope = new double[x.length - 1]; + double[] intercept = new double[x.length - 1]; + + // Calculate the line equation (i.e. slope and intercept) between each point + for (int i = 0; i < x.length - 1; i++) { + dx[i] = x[i + 1] - x[i]; + if (dx[i] == 0) { + throw new IllegalArgumentException("X must be montotonic. A duplicate " + "x-value was found"); + } + if (dx[i] < 0) { + throw new IllegalArgumentException("X must be sorted"); + } + dy[i] = y[i + 1] - y[i]; + slope[i] = dy[i] / dx[i]; + intercept[i] = y[i] - x[i] * slope[i]; + } + + // Perform the interpolation here + double[] yi = new double[xi.length]; + for (int i = 0; i < xi.length; i++) { + if ((xi[i] > x[x.length - 1]) || (xi[i] < x[0])) { + yi[i] = Double.NaN; + } + else { + int loc = Arrays.binarySearch(x, xi[i]); + if (loc < -1) { + loc = -loc - 2; + yi[i] = slope[loc] * xi[i] + intercept[loc]; + } + else { + yi[i] = y[loc]; + } + } + } + + return yi; + } + + public static double[] interpLinear(long[] x, double[] y, long[] xi) throws IllegalArgumentException { + + double[] xd = new double[x.length]; + for (int i = 0; i < x.length; i++) { + xd[i] = x[i]; + } + + double[] xid = new double[xi.length]; + for (int i = 0; i < xi.length; i++) { + xid[i] = xi[i]; + } + + return interpLinear(xd, y, xid); + } + + +} diff --git a/java/src/main/java/bpmEstimation/AccelerometerWindowBuffer.java b/java/src/main/java/bpmEstimation/AccelerometerWindowBuffer.java new file mode 100644 index 0000000..86e5c36 --- /dev/null +++ b/java/src/main/java/bpmEstimation/AccelerometerWindowBuffer.java @@ -0,0 +1,89 @@ +package bpmEstimation; + +import utilities.Utils; +import java.util.ArrayList; + +/** + * Created by toni on 15/12/17. + */ +public class AccelerometerWindowBuffer extends ArrayList { + + private int mWindowSize; + private int mOverlapSize; + private int mOverlapCounter; + private double[] mX; + private double[] mY; + private double[] mZ; + private long[] mTs; + + public AccelerometerWindowBuffer(int windowSize, int overlap){ + this.mWindowSize = windowSize; + this.mOverlapSize = overlap; + mOverlapCounter = 1; + + mX = new double[this.mWindowSize]; + mY = new double[this.mWindowSize]; + mZ = new double[this.mWindowSize]; + mTs = new long[this.mWindowSize]; + } + + //TODO: add exception handling. falseArgument if ad has no numeric x,y,z + public boolean add(AccelerometerData ad){ + + //do not add duplicates! + if(!isEmpty() && getYongest().equals(ad)){ + return false; + } + + boolean r = super.add(ad); + if (size() > mWindowSize){ + removeRange(0, size() - mWindowSize); + } + + //TODO: bullshit... jedes mal wieder alle neue anlegen ist ja nur langsam + //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; + } + + public boolean isNextWindowReady(){ + if(size() == mWindowSize && mOverlapCounter > mOverlapSize){ + mOverlapCounter = 1; + + return true; + } + return false; + } + + public AccelerometerData getYongest() { + return get(size() - 1); + } + + public AccelerometerData getOldest() { + return get(0); + } + + public double[] getX(){ + return mX; + } + + public double[] getY(){ + return mY; + } + + public double[] getZ(){ + return mZ; + } + + public long[] getTs(){ + return mTs; + } + +} \ No newline at end of file diff --git a/java/src/main/java/bpmEstimation/AutoCorrelation.java b/java/src/main/java/bpmEstimation/AutoCorrelation.java new file mode 100644 index 0000000..fe9c216 --- /dev/null +++ b/java/src/main/java/bpmEstimation/AutoCorrelation.java @@ -0,0 +1,79 @@ +package bpmEstimation; + +import org.jtransforms.fft.DoubleFFT_1D; +import utilities.Utils; +import java.util.Arrays; + +/** + * Created by toni on 15/12/17. + */ +public class AutoCorrelation { + + int mMaxLag; + double[] mCorr; + + public AutoCorrelation(double[] data, int maxLag){ + + mMaxLag = maxLag; + mCorr = fft(data); + } + + public double[] getCorr(){ + return mCorr; + } + + private double[] fft(double[] data) { + + if(mMaxLag < 1){ + throw new RuntimeException("maxlag has to be greater 1"); + } + + int n = data.length; + double[] x = Arrays.copyOf(data, n); + int mxl = Math.min(mMaxLag, n - 1); + int ceilLog2 = Utils.nextPow2(2*n -1); + int n2 = (int) Math.pow(2,ceilLog2); + + // x - mean(x) (pointwise) + double x_mean = Utils.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 + double[] x2 = new double[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 + DoubleFFT_1D fft = new DoubleFFT_1D(n2); + fft.realForwardFull(x2); + + // Cr = abs(x_fft).^2 (absolute with complex numbers is (r^2) + (i^2) + double[] Cr = new double[n2 * 2]; + int j = 0; + for(int i = 0; i < x2.length; ++i){ + Cr[j++] = Utils.sqr(x2[i]) + Utils.sqr(x2[i+1]); + ++i; //skip the complex part + } + + // ifft(Cr,[],1) + DoubleFFT_1D ifft = new DoubleFFT_1D(n2); + ifft.realInverseFull(Cr, true); + + // remove complex part and scale/normalize + double[] c1 = new double[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. + double[] c = new double[(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; + } +} diff --git a/java/src/main/java/bpmEstimation/Peaks.java b/java/src/main/java/bpmEstimation/Peaks.java new file mode 100644 index 0000000..cd7c6b6 --- /dev/null +++ b/java/src/main/java/bpmEstimation/Peaks.java @@ -0,0 +1,207 @@ +package bpmEstimation; + +import utilities.Utils; + +import java.util.LinkedList; + +/** + * Created by toni on 15/12/17. + */ +public class Peaks { + + private LinkedList mPeaksIdx; //provide the idx within the given data array + private LinkedList mPeaksPos; //the real position within the data-rang e.g. lag -1024 to 1024 + private LinkedList mPeaksValue; //the value at mPeaksPos + + private double[] mData; + + /** + * Simple method for finding local maxima in an array + * @param data input + * @param width minimum distance between peaks + * @param threshold minimum value of peaks + * @param decayRate how quickly previous peaks are forgotten + * @param isRelative minimum value of peaks is relative to local average + * @return array of peaks + */ + public Peaks(double[] data, int width, double threshold, double decayRate, boolean isRelative){ + + this.mData = data; + this.mPeaksIdx = new LinkedList<>(); + this.mPeaksPos = new LinkedList<>(); + this.mPeaksValue = new LinkedList<>(); + + //create the peaks + simplePeakFinder(data, width, threshold, decayRate, isRelative); + + updateLists(); + } + + public LinkedList getPeaksIdx() { + return mPeaksIdx; + } + + public int[] getPeaksIdxAsArray() { + return mPeaksIdx.stream().mapToInt(i->i).toArray(); + } + + public LinkedList getPeaksPos() { + return mPeaksPos; + } + + public double[] getPeaksPosAsArray() { + return mPeaksPos.stream().mapToDouble(i -> i).toArray(); + } + + public LinkedList getPeaksValue() { + return mPeaksValue; + } + + public double[] getPeaksValueAsArray() { + return mPeaksValue.stream().mapToDouble(i -> i).toArray(); + } + + public double[] getData(){ + return mData; + } + + //TODO: implement findFalseDetectedPeaks + public void improveResults(double[] data){ + + updateLists(); + } + + public double[] getPeaksValueWithoutZeroIdx() { + + double[] values = new double[mPeaksIdx.size() - 1]; + int i = 0; + int mid = (mData.length / 2); + for(Integer idx : mPeaksIdx) { + if(!idx.equals(mid)){ + values[i] = mPeaksValue.get(idx); + ++i; + } + } + + return values; + } + + public double[] getPeaksValueWithoutNegativeValues() { + + double[] values = new double[mPeaksIdx.size() - 1]; + int i = 0; + for(Integer idx : mPeaksIdx) { + double curVal = mPeaksValue.get(idx); + if(curVal > 0){ + values[i] = curVal; + ++i; + } + } + + return values; + } + + public double[] getPeaksValueWithoutZeroIdxAndNegativeValues(){ + double[] values = new double[mPeaksIdx.size() - 1]; + int i = 0; + int mid = (mData.length / 2); + for(Integer idx : mPeaksIdx) { + double curVal = mPeaksValue.get(idx); + if(!idx.equals(mid) && curVal> 0){ + values[i] = curVal; + ++i; + } + } + + return values; + } + + public boolean hasPeaks() { + return mPeaksIdx.size() > 1; + } + + /** + * Provides an estimation of beats per minute given a samplerate in milliseconds + * @param sampleRate_ms + * @return bpm if peaks found and conducting activity recognized, else -1 + */ + public double getBPM(double sampleRate_ms){ + + //todo: rückweisungsklasse kann auch hier mit rein. + if(hasPeaks()){ + + //todo diff and mean method for linkedlists for speed + return 60000 / Utils.mean(Utils.diff(mPeaksPos.stream().mapToDouble(i -> i * sampleRate_ms).toArray())); + } + + return -1; + } + + /** + * updates the position and values of the found peaks. + * call this if peaks are somewhat changed. + */ + private void updateLists(){ + //fill the position and the value lists + for(Integer idx : mPeaksIdx){ + int mid = (mData.length / 2); + mPeaksPos.add(idx - mid); + + mPeaksValue.add(mData[idx]); + } + } + + //TODO: findPeaks method identical to Matlab... with PeakProminence + + private void simplePeakFinder(double[] data, int width, double threshold, double decayRate, boolean isRelative) { + int maxp = 0; + int mid = 0; + int end = data.length; + double av = data[0]; + while (mid < end) { + av = decayRate * av + (1 - decayRate) * data[mid]; + if (av < data[mid]) + av = data[mid]; + int i = mid - width; + if (i < 0) + i = 0; + int stop = mid + width + 1; + if (stop > data.length) + stop = data.length; + maxp = i; + for (i++; i < stop; i++) + if (data[i] > data[maxp]) + maxp = i; + if (maxp == mid) { + if (overThreshold(data, maxp, width, threshold, isRelative, av)){ + this.mPeaksIdx.add(new Integer(maxp)); + } + } + mid++; + } + } + + private boolean overThreshold(double[] data, int index, int width, + double threshold, boolean isRelative, + double av) { + int pre = 3; + int post = 1; + + if (data[index] < av) + return false; + if (isRelative) { + int iStart = index - pre * width; + if (iStart < 0) + iStart = 0; + int iStop = index + post * width; + if (iStop > data.length) + iStop = data.length; + double sum = 0; + int count = iStop - iStart; + while (iStart < iStop) + sum += data[iStart++]; + return (data[index] > sum / count + threshold); + } else + return (data[index] > threshold); + } +} diff --git a/java/src/main/java/utilities/Utils.java b/java/src/main/java/utilities/Utils.java new file mode 100644 index 0000000..8464e31 --- /dev/null +++ b/java/src/main/java/utilities/Utils.java @@ -0,0 +1,232 @@ +package utilities; + +import bpmEstimation.Peaks; +import javax.swing.*; +import java.awt.*; +import java.awt.image.BufferedImage; +import java.util.Arrays; + +//TODO: change from double to generic type +public class Utils { + public static double getDistance(double x1, double y1, double x2, double y2) { + return (double) Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); + } + + public static double sqr(double x) { + return x * x; + } + + public static int nextPow2(int a){ + return a == 0 ? 0 : 32 - Integer.numberOfLeadingZeros(a - 1); + } + + public static double sum(double[] data){ + double sum = 0; + for (int i = 0; i < data.length; i++) { + sum += data[i]; + } + return sum; + } + + public static long sum(long[] data){ + long sum = 0; + for (int i = 0; i < data.length; i++) { + sum += data[i]; + } + return sum; + } + + //TODO: Could be slow.. faster method? + public static double rms(double[] nums) { + double sum = 0.0f; + for (double num : nums) + sum += num * num; + return (double) Math.sqrt(sum / nums.length); + } + + public static double[] diff(double[] data){ + double[] diff = new double[data.length - 1]; + int i=0; + for(int j = 1; j < data.length; ++j){ + diff[i] = data[j] - data[i]; + ++i; + } + return diff; + } + + public static long[] diff(long[] data){ + long[] diff = new long[data.length - 1]; + int i=0; + for(int j = 1; j < data.length; ++j){ + diff[i] = data[j] - data[i]; + ++i; + } + return diff; + } + + public static double mean(double[] data){ + return sum(data) / data.length; + } + + public static long mean(long[] data){ + return sum(data) / data.length; + } + + //TODO: Could be slow.. faster method? + public static double geometricMean(double[] data) { + double sum = data[0]; + for (int i = 1; i < data.length; i++) { + sum *= data[i]; + } + return (double) Math.pow(sum, 1.0 / data.length); + } + + public static int intersectionNumber(double[] signal, double border){ + int cnt = 0; + boolean isSmallerValue = false; + boolean isBiggerValue = false; + + for(double value : signal){ + if(value < border){ + if(isBiggerValue){ + cnt++; + } + + isSmallerValue = true; + isBiggerValue = false; + } + else { + if(isSmallerValue){ + cnt++; + } + + isSmallerValue = false; + isBiggerValue = true; + } + } + + return cnt; + } + + public static double getBestBpmEstimation(Peaks peaksX, Peaks peaksY, Peaks peaksZ) throws IllegalArgumentException { + + int cntNumAxis = 0; + double sumCorr = 1; //to prevent division by zero + double sumRms = 1; + double sumNumInter = 1; + + double corrMeanX = 0, corrRmsX = 0; + int corrNumInterX = 0; + if(peaksX.hasPeaks()){ + corrMeanX = geometricMean(peaksX.getPeaksValueWithoutZeroIdxAndNegativeValues()); + corrRmsX = rms(peaksX.getPeaksValueWithoutZeroIdx()); + corrNumInterX = intersectionNumber(peaksX.getData(), 0.2f); + + ++cntNumAxis; + sumCorr += corrMeanX; + sumRms += corrRmsX; + sumNumInter += corrNumInterX; + } + + double corrMeanY = 0, corrRmsY = 0; + int corrNumInterY = 0; + if(peaksY.hasPeaks()){ + corrMeanY = geometricMean(peaksY.getPeaksValueWithoutZeroIdxAndNegativeValues()); + corrRmsY = rms(peaksY.getPeaksValueWithoutZeroIdx()); + corrNumInterY = intersectionNumber(peaksY.getData(), 0.2f); + + ++cntNumAxis; + sumCorr += corrMeanY; + sumRms += corrRmsY; + sumNumInter += corrNumInterY; + } + + double corrMeanZ = 0, corrRmsZ = 0; + int corrNumInterZ = 0; + if(peaksZ.hasPeaks()){ + corrMeanZ = geometricMean(peaksZ.getPeaksValueWithoutZeroIdxAndNegativeValues()); + corrRmsZ = rms(peaksZ.getPeaksValueWithoutZeroIdx()); + corrNumInterZ = intersectionNumber(peaksZ.getData(), 0.2f); + + ++cntNumAxis; + sumCorr += corrMeanZ; + sumRms += corrRmsZ; + sumNumInter += corrNumInterZ; + } + + if(cntNumAxis == 0){ + throw new IllegalArgumentException("All Peaks are empty! -> Reject Estimation"); + } + + double quantityX = ((corrMeanX / sumCorr) + (corrRmsX / sumRms) + (corrNumInterX / sumNumInter)) / cntNumAxis; + double quantityY = ((corrMeanY / sumCorr) + (corrRmsY / sumRms) + (corrNumInterY / sumNumInter)) / cntNumAxis; + double quantityZ = ((corrMeanZ / sumCorr) + (corrRmsZ / sumRms) + (corrNumInterZ / sumNumInter)) / cntNumAxis; + + //get best axis by quantity and estimate bpm + if(quantityX > quantityY && quantityX > quantityZ){ + + } + else if(quantityY > quantityZ){ + + } + else { + + } + + //check if values are empty. + + return 0.0f; + } + + public static double[] removeZero(double[] array){ + int j = 0; + for( int i=0; i 0){ + array[j++] = array[i]; + } + } + double[] newArray = new double[j]; + System.arraycopy( array, 0, newArray, 0, j ); + + return newArray; + } + + @SuppressWarnings("serial") + public static class ShowPNG extends JFrame + { + JLabel mLabel; + ImageIcon mIcon; + + public ShowPNG(){ + mLabel = new JLabel(); + this.setLayout(new GridLayout(1,1)); + this.setSize(800, 640); + this.add(mLabel); + this.setVisible(true); + } + + public void set(BufferedImage bi){ + + mIcon = new ImageIcon(bi); + mLabel.setVisible(false); + mLabel.setIcon(mIcon); + mLabel.setVisible(true); + } + } + +}