diff --git a/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/MainActivity.java b/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/MainActivity.java index 7d0d796..55f274c 100644 --- a/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/MainActivity.java +++ b/android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/MainActivity.java @@ -247,6 +247,7 @@ public class MainActivity extends WearableActivity implements WorkerFragment.OnF mLastSendBPM = currentBPM; } + //künstlicher delay. if(mReadyToSend){ mReadyToSend = false; new java.util.Timer().schedule( diff --git a/android/ConductorsWatch/build.gradle b/android/ConductorsWatch/build.gradle index 9cbd907..2ea4df9 100644 --- a/android/ConductorsWatch/build.gradle +++ b/android/ConductorsWatch/build.gradle @@ -10,7 +10,7 @@ buildscript { google() } dependencies { - classpath 'com.android.tools.build:gradle:3.1.2' + classpath 'com.android.tools.build:gradle:3.1.3' // NOTE: Do not place your application dependencies here; they belong diff --git a/java/src/main/java/Main.java b/java/src/main/java/Main.java index d84ac21..2bd650b 100644 --- a/java/src/main/java/Main.java +++ b/java/src/main/java/Main.java @@ -21,28 +21,19 @@ public class Main { public static void main(String [ ] args) { //File folder = new File("/home/toni/Documents/programme/dirigent/measurements/2017.06/lgWear"); //File folder = new File("/home/toni/Documents/programme/dirigent/measurements/peter_failed"); - File folder = new File("/home/toni/Documents/programme/dirigent/measurements/2018.06/frank/mSensor"); + File folder = new File("/home/toni/Documents/programme/dirigent/measurements/2018.06/leon/mSensor"); File[] listOfFiles = folder.listFiles(); Arrays.sort(listOfFiles); - -// //TODO: write debug class, that is able to simply draw images... -// Utils.ShowPNG windowRaw = new Utils.ShowPNG(); -// Utils.ShowPNG windowAuto = new Utils.ShowPNG(); -// Utils.ShowPNG windowAutoButter = new Utils.ShowPNG(); -// Utils.ShowPNG windowPeaksX = new Utils.ShowPNG(); -// Utils.ShowPNG windowPeaksY = new Utils.ShowPNG(); -// Utils.ShowPNG windowPeaksZ = new Utils.ShowPNG(); + //calc results // iterate trough files in measurements folder for (File file : listOfFiles) { if (file.isFile() && file.getName().contains(".csv")) { - //TODO: Die Raw Sensordaten sollte man vielleicht etwas glätten. Sieh Aufnahmen von Frank! AccelerometerWindowBuffer accWindowBuffer = new AccelerometerWindowBuffer(6000, 750); BpmEstimator bpmEstimator = new BpmEstimator(accWindowBuffer, 0, 5000); - Butterworth butterLowpass = new Butterworth(); //read the file line by line try (BufferedReader br = new BufferedReader(new FileReader(file))) { @@ -50,16 +41,18 @@ public class Main { //read the first three lines and print out what file it is! String comment = br.readLine(); br.readLine(); - br.readLine(); + String groundTruth = br.readLine(); System.out.println(comment); + //long startTs = Long.parseLong(br.readLine().split(";")[0]); for (String line; (line = br.readLine()) != null; ) { // process the line. String[] measurement = line.split(";"); //if linear acc + long ts = 0; if(measurement[1].equals("3")){ - long ts = Long.parseLong(measurement[0]); + ts = Long.parseLong(measurement[0]); double x = Double.parseDouble(measurement[2]); double y = Double.parseDouble(measurement[3]); double z = Double.parseDouble(measurement[4]); @@ -88,120 +81,18 @@ public class Main { bpmList.add(bpm160); bpmList.add(bpm200); -// //TODO: plot the different autocorrelation of different window sizes -// // Draw the data for the complete 6 second window -// double sampleRate = 4; -// AccelerometerInterpolator acInterp = new AccelerometerInterpolator(accWindowBuffer, sampleRate); -// int peakWidth = (int) Math.round(250 / sampleRate); -// -// //print raw x,y,z -// double[] dTs = IntStream.range(0, accWindowBuffer.getTs().length).mapToDouble(i -> accWindowBuffer.getTs()[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(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 -// double[] xAutoCorr = new AutoCorrelation(acInterp.getX(), accWindowBuffer.size()).getCorr(); -// double[] yAutoCorr = new AutoCorrelation(acInterp.getY(), accWindowBuffer.size()).getCorr(); -// double[] zAutoCorr = new AutoCorrelation(acInterp.getZ(), accWindowBuffer.size()).getCorr(); -// -// //print autocorr raw -// 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(); -// Plot plotCorr = Plot.plot(Plot.plotOpts(). -// title("Auto Correlation"). -// legend(Plot.LegendFormat.BOTTOM)). -// 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()); -// -// //print autocorr butter -// // butterworth lowpass filter, cutoff at 3 hz (~180 bpm) -// butterLowpass.lowPass(1,(1/(sampleRate / 1000))/2, 1); -// -// int n = acInterp.getX().length; -// double[] xButter = new double[n]; -// double[] yButter = new double[n];; -// double[] zButter = new double[n];; -// for(int i = 0; i < acInterp.getX().length; ++i){ -// //xButter[i] = butterLowpass.filter(acInterp.getX()[i]); -// //yButter[i] = butterLowpass.filter(acInterp.getY()[i]); -// zButter[i] = butterLowpass.filter(acInterp.getZ()[i]); -// } -// -// //double[] xAutoCorrButter = new AutoCorrelation(xButter, accWindowBuffer.size()).getCorr(); -// //double[] yAutoCorrButter = new AutoCorrelation(yButter, accWindowBuffer.size()).getCorr(); -// double[] zAutoCorrButter = new AutoCorrelation(zButter, accWindowBuffer.size()).getCorr(); -// -// Plot plotCorrButter = Plot.plot(Plot.plotOpts(). -// title("Auto Correlation Butter"). -// legend(Plot.LegendFormat.BOTTOM)). -// //series("x", utilities.Plot.data().xy(rangeAuto, xAutoCorrButter), utilities.Plot.seriesOpts().color(Color.RED)). -// //series("y", utilities.Plot.data().xy(rangeAuto, yAutoCorrButter), utilities.Plot.seriesOpts().color(Color.BLUE)). -// series("z", Plot.data().xy(rangeAuto, zAutoCorrButter), Plot.seriesOpts().color(Color.GREEN)); -// -// windowAutoButter.set(plotCorrButter.draw()); -// -// //Print peaks -// Peaks pX = new Peaks(xAutoCorr, peakWidth, 0.1f, 0, false); -// int xOffset = xAutoCorr.length / 2; -// LinkedList peaksX = pX.getPeaksIdx(); -// -// double[] dPeaksXX = IntStream.range(0, peaksX.size()).mapToDouble(i -> (peaksX.get(i) - xOffset)).toArray();//peaks.stream().mapToDouble(i->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, 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()); -// -// Peaks pY = new Peaks(yAutoCorr, peakWidth, 0.1f, 0, false); -// int yOffset = yAutoCorr.length / 2; -// LinkedList peaksY = pY.getPeaksIdx(); -// -// double[] dPeaksYX = IntStream.range(0, peaksY.size()).mapToDouble(i -> (peaksY.get(i) - yOffset)).toArray();//peaks.stream().mapToDouble(i->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, 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()); -// -// Peaks pZ = new Peaks(zAutoCorr, peakWidth, 0.1f, 0, false); -// int zOffset = zAutoCorr.length / 2; -// LinkedList peaksZ = pZ.getPeaksIdx(); -// -// double[] dPeaksZX = IntStream.range(0, peaksZ.size()).mapToDouble(i -> (peaksZ.get(i) - zOffset)).toArray();//peaks.stream().mapToDouble(i->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, 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)); -// -// windowPeaksZ.set(plotPeaksZ.draw()); + while(bpmList.remove(Double.valueOf(-1))) {} + Utils.removeOutliersZScore(bpmList, 3.4); - //fill hols improve peaks + double bpmMean = Utils.mean(bpmList); + double magMean = bpmEstimator.getMagnitudeMean(); + double bpmDist = bpmEstimator.getDistEstimation(); - //estimate bpm between detected peaks - //System.out.println("BPM-X: " + pX.getBPM(bpmEstimator.getSampleRate_ms())); - //System.out.println("BPM-Y: " + pY.getBPM(bpmEstimator.getSampleRate_ms())); - //System.out.println("BPM-Z: " + pZ.getBPM(bpmEstimator.getSampleRate_ms())); + //double bpmSingle = bpmEstimator.getBestSingleAxis(); + double bpmAllAverage = bpmEstimator.getAverageOfAllWindows(); + + System.out.println( ts + " all: " + Math.round(bpmMean) + " avg_all: " + Math.round(bpmAllAverage) + " 3D: " + Math.round(bpmDist)); + System.out.println(" "); int dummyForBreakpoint = 0; } @@ -213,17 +104,24 @@ public class Main { //System.out.println("MEAN BPM: " + Math.round(meanBPM)); //System.out.println("MEDIAN BPM: " + Math.round(medianBPM)); + if(Utils.DEBUG_MODE){ + bpmEstimator.closeDebugWindows(); + } + + } catch (IOException e) { e.printStackTrace(); } } - try { - System.in.read(); - } catch (IOException e) { - e.printStackTrace(); - } + + +// try { +// System.in.read(); +// } catch (IOException e) { +// e.printStackTrace(); +// } } } } diff --git a/java/src/main/java/bpmEstimation/AccelerometerInterpolator.java b/java/src/main/java/bpmEstimation/AccelerometerInterpolator.java index e270eea..3940a6d 100644 --- a/java/src/main/java/bpmEstimation/AccelerometerInterpolator.java +++ b/java/src/main/java/bpmEstimation/AccelerometerInterpolator.java @@ -11,11 +11,12 @@ public class AccelerometerInterpolator { private double[] mY; private double[] mZ; private long[] mTsInterp; + private int size; public AccelerometerInterpolator(AccelerometerWindowBuffer ab, double sampleRate_ms){ - long size = (ab.getYongest().ts - (ab.getOldest().ts - (long) sampleRate_ms)) / (long) sampleRate_ms; - mTsInterp = new long[(int)size]; + this.size = (int) ((ab.getYongest().ts - (ab.getOldest().ts - (long) sampleRate_ms)) / (long) sampleRate_ms); + mTsInterp = new long[size]; int j = 0; for(long i = ab.getOldest().ts; i <= ab.getYongest().ts; i += sampleRate_ms){ mTsInterp[j++] = i; @@ -42,6 +43,10 @@ public class AccelerometerInterpolator { return mZ; } + public int size(){ + return size; + } + private 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"); diff --git a/java/src/main/java/bpmEstimation/BpmEstimator.java b/java/src/main/java/bpmEstimation/BpmEstimator.java index 75163d5..7531e6e 100644 --- a/java/src/main/java/bpmEstimation/BpmEstimator.java +++ b/java/src/main/java/bpmEstimation/BpmEstimator.java @@ -3,7 +3,6 @@ package bpmEstimation; import static utilities.Utils.DEBUG_MODE; import utilities.MovingFilter; -import utilities.SimpleKalman; import utilities.Utils; import uk.me.berndporr.iirj.*; @@ -18,10 +17,11 @@ public class BpmEstimator { private AccelerometerWindowBuffer mBuffer; private static double mSampleRate_ms; - private LinkedList mBpmHistory_X; - private LinkedList mBpmHistory_Y; - private LinkedList mBpmHistory_Z; - private LinkedList mBpmHistory_Mag; + private BpmHistory mBpmHistory_X; + private BpmHistory mBpmHistory_Y; + private BpmHistory mBpmHistory_Z; + private BpmHistory mBpmHistory_Mag; + private BpmHistory mBpmHistory_Dist; private LinkedList mBpmHistory; private int mResetCounter; @@ -35,6 +35,10 @@ public class BpmEstimator { private Butterworth mButter_Z; private Butterworth mButter_Mag; + + //todo: hack + private double magnitudeMean = 1.0; + //Debugging stuff private Utils.DebugPlotter plotter; @@ -42,12 +46,13 @@ public class BpmEstimator { mBuffer = windowBuffer; mSampleRate_ms = sampleRate_ms; - mBpmHistory_X = new LinkedList<>(); - mBpmHistory_Y = new LinkedList<>(); - mBpmHistory_Z = new LinkedList<>(); - mBpmHistory_Mag = new LinkedList<>(); + mBpmHistory_X = new BpmHistory(); + mBpmHistory_Y = new BpmHistory(); + mBpmHistory_Z = new BpmHistory(); + mBpmHistory_Mag = new BpmHistory(); + mBpmHistory_Dist = new BpmHistory(); - mBpmHistory = new LinkedList<>(); + mBpmHistory = new BpmHistory(); mResetCounter = 0; mResetLimit_ms = resetAfter_ms; @@ -70,6 +75,8 @@ public class BpmEstimator { public double estimate(int length_ms, int overlap_ms){ + //todo: interpolator könnte man direkt vom buffer zurück geben lassen. + //todo: anstatt hier eine zweite instance zu öffnen, einfach ein globales format wie bpsw: AccelerometerWindow AccelerometerWindowBuffer tmpBuffer = AccelerometerWindowBuffer.getNewInstance(mBuffer, length_ms, overlap_ms); double sampleRate = mSampleRate_ms; @@ -89,25 +96,29 @@ public class BpmEstimator { mButter_X.lowPass(1,(1/(sampleRate / 1000))/2, 1); mButter_Y.lowPass(1,(1/(sampleRate / 1000))/2, 1); mButter_Z.lowPass(1,(1/(sampleRate / 1000))/2, 1); + mButter_Mag.lowPass(1,(1/(sampleRate / 1000))/2, 1); int n = interp.getX().length; double[] xButter = new double[n]; - double[] yButter = new double[n];; - double[] zButter = new double[n];; + double[] yButter = new double[n]; + double[] zButter = new double[n]; + double[] magButter = new double[n]; for(int i = 0; i < interp.getX().length; ++i){ xButter[i] = mButter_X.filter(interp.getX()[i]); yButter[i] = mButter_Y.filter(interp.getY()[i]); zButter[i] = mButter_Z.filter(interp.getZ()[i]); + magButter[i] = mButter_Mag.filter(magRaw[i]); } - //TODO: compare with own mButter_Mag.filter - double[] magButter = Utils.magnitude(xButter, yButter, zButter); //auto correlation double[] xAutoCorr = new AutoCorrelation(xButter, tmpBuffer.size()).getCorr(); double[] yAutoCorr = new AutoCorrelation(yButter, tmpBuffer.size()).getCorr(); double[] zAutoCorr = new AutoCorrelation(zButter, tmpBuffer.size()).getCorr(); - double[] magAutoCorr = new AutoCorrelation(magRaw, tmpBuffer.size()).getCorr(); + double[] magAutoCorr = new AutoCorrelation(magButter, tmpBuffer.size()).getCorr(); + + //dist correlation + double[] distCorr = new DistanceCorrelation(interp, (int) (interp.size() * 0.8)).getCorr(); //find a peak within range of 250 ms int peakWidth = (int) Math.round(250 / sampleRate); @@ -116,36 +127,44 @@ public class BpmEstimator { Peaks zPeaks = new Peaks(zAutoCorr, peakWidth, 0.1f, 0, false); Peaks magPeaks = new Peaks(magAutoCorr, peakWidth, 0.1f, 0, false); + Peaks distPeaks = new Peaks(distCorr, peakWidth, 0.1f, 0, false); + mBpmHistory_X.add(xPeaks.getBPM(sampleRate)); mBpmHistory_Y.add(yPeaks.getBPM(sampleRate)); mBpmHistory_Z.add(zPeaks.getBPM(sampleRate)); mBpmHistory_Mag.add(magPeaks.getBPM(sampleRate)); + mBpmHistory_Dist.add(distPeaks.getBPM(sampleRate)); if(DEBUG_MODE){ - plotter.setPlotRawX(interp.getTs(), interp.getX()); - plotter.setPlotRawY(interp.getTs(), interp.getY()); - plotter.setPlotRawZ(interp.getTs(), interp.getZ()); - plotter.setPlotRawMag(interp.getTs(), magRaw); + //plotter.setPlotRawX(interp.getTs(), interp.getX()); + //plotter.setPlotRawY(interp.getTs(), interp.getY()); + //plotter.setPlotRawZ(interp.getTs(), interp.getZ()); + //plotter.setPlotRawMag(interp.getTs(), magRaw); - plotter.setPlotButterX(interp.getTs(), xButter); - plotter.setPlotButterY(interp.getTs(), yButter); - plotter.setPlotButterZ(interp.getTs(), zButter); - plotter.setPlotButterMag(interp.getTs(), magButter); + //plotter.setPlotButterX(interp.getTs(), xButter); + //plotter.setPlotButterY(interp.getTs(), yButter); + //plotter.setPlotButterZ(interp.getTs(), zButter); + //plotter.setPlotButterMag(interp.getTs(), magButter); plotter.setPlotCorrX(xAutoCorr, xPeaks); plotter.setPlotCorrY(yAutoCorr, yPeaks); plotter.setPlotCorrZ(zAutoCorr, zPeaks); plotter.setPlotCorrMag(magAutoCorr, magPeaks); + plotter.setPlotCorr3D(distCorr, distPeaks); - //printout the current BPM - System.out.println(length_ms + "; x: " + mBpmHistory_X.getLast() - + "; y: " + mBpmHistory_Y.getLast() - + "; z: " + mBpmHistory_Z.getLast() - + "; mag: " + mBpmHistory_Mag.getLast()); } - double estimatedBPM = getBestBpmEstimation(xPeaks, yPeaks, zPeaks); + //printout the current BPM + System.out.println(length_ms + "; x: " + Math.round(xPeaks.getBPM(sampleRate)) + + "; y: " + Math.round(yPeaks.getBPM(sampleRate)) + + "; z: " + Math.round(zPeaks.getBPM(sampleRate)) + + "; mag: " + Math.round(magPeaks.getBPM(sampleRate)) + + "; 3D: " + Math.round(distPeaks.getBPM(sampleRate))); + + + + double estimatedBPM = getBestBpmEstimation(xPeaks, yPeaks, zPeaks, magPeaks); if(estimatedBPM != -1){ //moving avg (lohnt dann, wenn wir viele daten haben) @@ -178,6 +197,22 @@ public class BpmEstimator { return mBpmHistory.getLast(); } + public double getDistEstimation(){ + BpmHistory tmp = (BpmHistory) mBpmHistory_Dist.clone(); + tmp.removeOutliers(); + return tmp.getMean(); + } + + public double getMagnitudeMean(){ + //while(mBpmHistory_Mag.remove(Double.valueOf(-1))) {} + //Utils.removeOutliersZScore(mBpmHistory_Mag, 3.4); + //double mean = Utils.mean(mBpmHistory_Mag); + //mBpmHistory_Mag.clear(); + BpmHistory tmp = (BpmHistory) mBpmHistory_Mag.clone(); + tmp.removeOutliers(); + return tmp.getMean(); + } + public double getMeanBpm(){ return Utils.mean(mBpmHistory); } @@ -186,9 +221,84 @@ public class BpmEstimator { return Utils.median(mBpmHistory); } + //call this after the different estimations are done + //1) nimm die achse mit der geringsten varianz der bpm schätzung über alle fenster + //2) schmeiße outliers raus + //3) berechne mean -> profit?! + //=> ist bisschen willkürlich... da natürlich auch eine falsche achse stabil sein kann. meist ist es dann genau die hälfte. + public double getBestSingleAxis(){ + + mBpmHistory_X.removeOutliers(); + mBpmHistory_Y.removeOutliers(); + mBpmHistory_Z.removeOutliers(); + mBpmHistory_Mag.removeOutliers(); + + double varX = mBpmHistory_X.getVariance(); + double varY = mBpmHistory_Y.getVariance(); + double varZ = mBpmHistory_Z.getVariance(); + double varM = mBpmHistory_Mag.getVariance(); + double mean; + + + //TODO: nimm den der am wenigsten streut und der, der sich am wenigsten von der letzten messung unterscheidet. + //TODO: gewichte bpm schätzung von größeren fenstern höher. + + if(varX < varY && varX < varZ && varX < varM){ + + mean = mBpmHistory_X.getMean(); + } else if (varY < varZ && varY < varM) { + + mean = mBpmHistory_Y.getMean(); + } else if (varZ < varM){ + + mean = mBpmHistory_Z.getMean(); + } else { + + mean = mBpmHistory_Mag.getMean(); + } + + //todo: this is so ugly... needs to be refactored in app. + mBpmHistory_X.clear(); + mBpmHistory_Y.clear(); + mBpmHistory_Z.clear(); + mBpmHistory_Mag.clear(); + + return mean; + } + + + public double getAverageOfAllWindows(){ + + //remove outliers solo + //mBpmHistory_X.removeOutliers(); + //mBpmHistory_Y.removeOutliers(); + //mBpmHistory_Z.removeOutliers(); + //mBpmHistory_Mag.removeOutliers(); + + //write everything in a single vector + BpmHistory tmpHistory = new BpmHistory(); + tmpHistory.add(mBpmHistory_X); + tmpHistory.add(mBpmHistory_Y); + tmpHistory.add(mBpmHistory_Z); + tmpHistory.add(mBpmHistory_Mag); + + //remove outliers again + tmpHistory.removeOutliers(); + + //clear + mBpmHistory_X.clear(); + mBpmHistory_Y.clear(); + mBpmHistory_Z.clear(); + mBpmHistory_Mag.clear(); + mBpmHistory_Dist.clear(); + + return tmpHistory.getMean(); + + } + //TODO: die einzelnen achsen cleverer kombinieren. bspw. bei Peter brauchen wir zwei Achsen in einer Bewegung. //TODO: Vielleicht die einzelnen Kombinationen / Magnitudes der Achsen noch mit einbeziehen. Also xy, xz, yz - private double getBestBpmEstimation(Peaks peaksX, Peaks peaksY, Peaks peaksZ) throws IllegalArgumentException { + private double getBestBpmEstimation(Peaks peaksX, Peaks peaksY, Peaks peaksZ, Peaks peaksMag) throws IllegalArgumentException { int cntNumAxis = 0; double sumCorr = 0; //to prevent division by zero @@ -234,6 +344,19 @@ public class BpmEstimator { sumNumInter += corrNumInterZ; } +// double corrMeanMag = 0, corrRmsMag = 0; +// int corrNumInterMag = 0; +// if(peaksMag.hasPeaks()){ +// corrMeanMag = Utils.geometricMean(peaksMag.getPeaksValueWithoutZeroIdxAndNegativeValues()); +// corrRmsMag = Utils.rms(peaksMag.getPeaksValueWithoutZeroIdx()); +// corrNumInterMag = Utils.intersectionNumber(peaksMag.getData(), 0.2f); +// +// ++cntNumAxis; +// sumCorr += corrMeanMag; +// sumRms += corrRmsMag; +// sumNumInter += corrNumInterMag; +// } + //no peaks, reject if(cntNumAxis == 0){ //throw new IllegalArgumentException("All Peaks are empty! -> Reject Estimation"); @@ -256,23 +379,28 @@ public class BpmEstimator { //values to low, reject //TODO: this is a pretty simple assumption. first shot! - if(corrRmsX < 0.25 && corrRmsY < 0.25 && corrRmsZ < 0.25){ + if(corrRmsX < 0.25 && corrRmsY < 0.25 && corrRmsZ < 0.25 ){ return -1; } 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; + //double quantityMag = ((corrMeanMag / sumCorr) + (corrRmsMag / sumRms) + (corrNumInterMag / sumNumInter)) / cntNumAxis; //get best axis by quantity and estimate bpm if(quantityX > quantityY && quantityX > quantityZ){ return mBpmHistory_X.getLast(); } - else if(quantityY > quantityZ){ + else if(quantityY > quantityZ ){ return mBpmHistory_Y.getLast(); } else { return mBpmHistory_Z.getLast(); } } + + public void closeDebugWindows(){ + plotter.close(); + } } diff --git a/java/src/main/java/bpmEstimation/BpmHistory.java b/java/src/main/java/bpmEstimation/BpmHistory.java new file mode 100644 index 0000000..12896a7 --- /dev/null +++ b/java/src/main/java/bpmEstimation/BpmHistory.java @@ -0,0 +1,53 @@ +package bpmEstimation; + +import utilities.Utils; + +import java.util.LinkedList; + +/** + * Created by toni on 03/08/18. + */ +public class BpmHistory extends LinkedList { + + private double mMean = 0.0; + private double mVar = 0.0; + private double mStdDev = 0.0; + + public boolean add(double val){ + if(val != -1){ + super.add(val); + return true; + } + return false; + } + + public boolean add(BpmHistory vals){ + if(!vals.isEmpty()){ + for(double val : vals){ + super.add(val); + } + return true; + } + return false; + } + + public double getMean(){ + if(this.size() > 2){ + return Utils.mean(this); + } else { + return 333; //TODO: das ist natürlich quatsch und faulheit. mal schaun wie man das am besten löst. + } + } + + public double getVariance(){ + if(this.size() > 2){ + return Utils.var(this); + } else { + return 666; //TODO: das ist natürlich quatsch und faulheit. mal schaun wie man das am besten löst. + } + } + + public void removeOutliers(){ + Utils.removeOutliersZScore(this, 3.4); + } +} diff --git a/java/src/main/java/bpmEstimation/DistanceCorrelation.java b/java/src/main/java/bpmEstimation/DistanceCorrelation.java new file mode 100644 index 0000000..3b84419 --- /dev/null +++ b/java/src/main/java/bpmEstimation/DistanceCorrelation.java @@ -0,0 +1,66 @@ +package bpmEstimation; + +import utilities.Utils; + +import java.util.Arrays; +import java.util.Collections; +import java.util.DoubleSummaryStatistics; + +/** + * Created by toni on 06/12/18. + */ +public class DistanceCorrelation { + + private static int mMaxLag; + private double[] mCorr; + + public DistanceCorrelation(AccelerometerInterpolator data, int maxLag){ + + mMaxLag = maxLag; + mCorr = calc(data); + } + + public double[] getCorr(){ + return mCorr; + } + + private double[] calc(AccelerometerInterpolator data){ + + if(mMaxLag < 1){ + throw new RuntimeException("maxlag has to be greater 1"); + } + + //init + int n = data.size(); + int lag_size = Math.min(mMaxLag, n - 1); + double[] corr = new double[lag_size + 1]; + + //do the math + for(int j = 0; j <= lag_size; ++j){ + double[] dist = new double[n - Math.abs(j)]; + int idx = 0; + + for (int i = j; i < n; ++i){ + dist[idx] = Utils.getDistance(data.getX()[i], data.getY()[i], data.getZ()[i], data.getX()[i-j], data.getY()[i-j], data.getZ()[i-j]); + ++idx; + } + + corr[j] = Utils.geometricMeanLog(dist); + } + + //to [0, 1] + DoubleSummaryStatistics corrStat = Arrays.stream(corr).summaryStatistics(); + double corMaxVal = corrStat.getMax(); + for(int k = 0; k < corr.length; ++k){ + corr[k] = ((corr[k] * (-1)) / corMaxVal) + 1; + } + + // mirror corr(2:512) and put it in front + double[] output = new double[(2 * lag_size) + 1]; + System.arraycopy(corr, 0, output, lag_size, lag_size + 1); // +1 to place the 1.0 in the middle of correlation + Utils.reverse(corr); + System.arraycopy(corr, 0, output, 0, lag_size); + + return output; + } +} diff --git a/java/src/main/java/utilities/Plot.java b/java/src/main/java/utilities/Plot.java new file mode 100644 index 0000000..51f4002 --- /dev/null +++ b/java/src/main/java/utilities/Plot.java @@ -0,0 +1,1029 @@ +package utilities; + +import java.awt.BasicStroke; +import java.awt.Color; +import java.awt.Font; +import java.awt.FontMetrics; +import java.awt.Graphics2D; +import java.awt.Point; +import java.awt.Polygon; +import java.awt.Rectangle; +import java.awt.Stroke; +import java.awt.geom.Rectangle2D; +import java.awt.image.BufferedImage; +import java.io.File; +import java.io.IOException; +import java.util.ArrayList; +import java.util.HashMap; +import java.util.Iterator; +import java.util.LinkedHashMap; +import java.util.List; +import java.util.Map; + +import javax.imageio.ImageIO; + +/** + * Simple implementation of plot. Minimal features, no dependencies besides standard libraries. + * Options are self-descriptive, see also samples. + * + * @author Yuriy Guskov + */ +public class Plot { + + public enum Line { NONE, SOLID, DASHED } + public enum Marker { NONE, CIRCLE, SQUARE, DIAMOND, COLUMN, BAR } + public enum AxisFormat { NUMBER, NUMBER_KGM, NUMBER_INT, TIME_HM, TIME_HMS, DATE, DATETIME_HM, DATETIME_HMS } + public enum LegendFormat { NONE, TOP, RIGHT, BOTTOM } + + private enum HorizAlign { LEFT, CENTER, RIGHT } + private enum VertAlign { TOP, CENTER, BOTTOM } + + private PlotOptions opts = new PlotOptions(); + + private Rectangle boundRect; + private PlotArea plotArea; + private Map xAxes = new HashMap(3); + private Map yAxes = new HashMap(3); + private Map dataSeriesMap = new LinkedHashMap(5); + + public static Plot plot(PlotOptions opts) { + return new Plot(opts); + } + + public static PlotOptions plotOpts() { + return new PlotOptions(); + } + + public static class PlotOptions { + + private String title = ""; + private int width = 800; + private int height = 600; + private Color backgroundColor = Color.WHITE; + private Color foregroundColor = Color.BLACK; + private Font titleFont = new Font("Arial", Font.BOLD, 16); + private int padding = 10; // padding for the entire image + private int plotPadding = 5; // padding for plot area (to have min and max values padded) + private int labelPadding = 10; + private int defaultLegendSignSize = 10; + private int legendSignSize = 10; + private Point grids = new Point(10 ,10); // grid lines by x and y + private Color gridColor = Color.GRAY; + private Stroke gridStroke = new BasicStroke(1.0f, BasicStroke.CAP_BUTT, + BasicStroke.JOIN_MITER, 10.0f, new float[] { 5.0f }, 0.0f); + private int tickSize = 5; + private Font labelFont = new Font("Arial", 0, 12); + private LegendFormat legend = LegendFormat.NONE; + + private PlotOptions() {} + + public PlotOptions title(String title) { + this.title = title; + return this; + } + + public PlotOptions width(int width) { + this.width = width; + return this; + } + + public PlotOptions height(int height) { + this.height = height; + return this; + } + + public PlotOptions bgColor(Color color) { + this.backgroundColor = color; + return this; + } + + public PlotOptions fgColor(Color color) { + this.foregroundColor = color; + return this; + } + + public PlotOptions titleFont(Font font) { + this.titleFont = font; + return this; + } + + public PlotOptions padding(int padding) { + this.padding = padding; + return this; + } + + public PlotOptions plotPadding(int padding) { + this.plotPadding = padding; + return this; + } + + public PlotOptions labelPadding(int padding) { + this.labelPadding = padding; + return this; + } + + public PlotOptions labelFont(Font font) { + this.labelFont = font; + return this; + } + + public PlotOptions grids(int byX, int byY) { + this.grids = new Point(byX, byY); + return this; + } + + public PlotOptions gridColor(Color color) { + this.gridColor = color; + return this; + } + + public PlotOptions gridStroke(Stroke stroke) { + this.gridStroke = stroke; + return this; + } + + public PlotOptions tickSize(int value) { + this.tickSize = value; + return this; + } + + public PlotOptions legend(LegendFormat legend) { + this.legend = legend; + return this; + } + + } + + private Plot(PlotOptions opts) { + if (opts != null) + this.opts = opts; + boundRect = new Rectangle(0, 0, this.opts.width, this.opts.height); + plotArea = new PlotArea(); + } + + public PlotOptions opts() { + return opts; + } + + public Plot xAxis(String name, AxisOptions opts) { + xAxes.put(name, new Axis(name, opts)); + return this; + } + + public Plot yAxis(String name, AxisOptions opts) { + yAxes.put(name, new Axis(name, opts)); + return this; + } + + public Plot series(String name, Data data, DataSeriesOptions opts) { + DataSeries series = dataSeriesMap.get(name); + if (opts != null) + opts.setPlot(this); + if (series == null) { + series = new DataSeries(name, data, opts); + dataSeriesMap.put(name, series); + } else { + series.data = data; + series.opts = opts; + } + return this; + } + + public Plot series(String name, DataSeriesOptions opts) { + DataSeries series = dataSeriesMap.get(name); + if (opts != null) + opts.setPlot(this); + if (series != null) + series.opts = opts; + return this; + } + + private void calc(Graphics2D g) { + plotArea.calc(g); + } + + private void clear() { + plotArea.clear(); + for (DataSeries series : dataSeriesMap.values()) + series.clear(); + } + + public BufferedImage draw() { + BufferedImage image = new BufferedImage(opts.width, opts.height, BufferedImage.TYPE_INT_RGB); + Graphics2D g = image.createGraphics(); + try { + calc(g); + drawBackground(g); + plotArea.draw(g); + for (DataSeries series : dataSeriesMap.values()) + series.draw(g); + return image; + } finally { + g.dispose(); + } + } + + private void drawBackground(Graphics2D g) { + g.setColor(opts.backgroundColor); + g.fillRect(0, 0, opts.width, opts.height); + } + + public void save(String fileName, String type) throws IOException { + clear(); + BufferedImage bi = draw(); + File outputFile = new File(fileName + "." + type); + ImageIO.write(bi, type, outputFile); + } + + private class Legend { + Rectangle rect; + Rectangle2D labelRect; + public int entryWidth; + public int entryWidthPadded; + public int entryCount; + public int xCount; + public int yCount; + } + + private class PlotArea { + + private Rectangle plotBorderRect = new Rectangle(); // boundRect | labels/legend | plotBorderRect | plotPadding | plotRect/clipRect + private Rectangle plotRect = new Rectangle(); + private Rectangle plotClipRect = new Rectangle(); + private Legend legend = new Legend(); + + private Range xPlotRange = new Range(0, 0); + private Range yPlotRange = new Range(0, 0); + + public PlotArea() { + clear(); + } + + private void clear() { + plotBorderRect.setBounds(boundRect); + plotRectChanged(); + } + + private void offset(int dx, int dy, int dw, int dh) { + plotBorderRect.translate(dx, dy); + plotBorderRect.setSize(plotBorderRect.width - dx - dw, plotBorderRect.height - dy - dh); + plotRectChanged(); + } + + private void plotRectChanged() { + plotRect.setBounds(plotBorderRect.x + opts.plotPadding, plotBorderRect.y + opts.plotPadding, + plotBorderRect.width - opts.plotPadding * 2, plotBorderRect.height - opts.plotPadding * 2); + xPlotRange.setMin(plotRect.getX()); + xPlotRange.setMax(plotRect.getX() + plotRect.getWidth()); + yPlotRange.setMin(plotRect.getY()); + yPlotRange.setMax(plotRect.getY() + plotRect.getHeight()); + + plotClipRect.setBounds(plotBorderRect.x + 1, plotBorderRect.y + 1, plotBorderRect.width - 1, plotBorderRect.height - 1); + } + + private void calc(Graphics2D g) { + calcAxes(g); + calcRange(true); + calcRange(false); + calcAxisLabels(g, true); + calcAxisLabels(g, false); + g.setFont(opts.titleFont); + FontMetrics fm = g.getFontMetrics(); + Rectangle2D titleRect = fm.getStringBounds(opts.title, g); + g.setFont(opts.labelFont); + fm = g.getFontMetrics(); + int xAxesHeight = 0, xAxesHalfWidth = 0; + for (Map.Entry entry : xAxes.entrySet()) { + Axis xAxis = entry.getValue(); + xAxesHeight += toInt(xAxis.labelRect.getHeight()) + opts.labelPadding * 2; + if (xAxis.labelRect.getWidth() > xAxesHalfWidth) + xAxesHalfWidth = toInt(xAxis.labelRect.getWidth()); + } + int yAxesWidth = 0; + for (Map.Entry entry : yAxes.entrySet()) + yAxesWidth += toInt(entry.getValue().labelRect.getWidth()) + opts.labelPadding * 2; + int dx = opts.padding + yAxesWidth; + int dy = opts.padding + toInt(titleRect.getHeight() + opts.labelPadding); + int dw = opts.padding; + if (opts.legend != LegendFormat.RIGHT) + dw += xAxesHalfWidth; // half of label goes beyond a plot in right bottom corner + int dh = opts.padding + xAxesHeight; + // offset for legend + Rectangle temp = new Rectangle(plotBorderRect); // save plotRect + offset(dx, dy, dw, dh); + calcLegend(g); // use plotRect + plotBorderRect.setBounds(temp); // restore plotRect + switch (opts.legend) { + case TOP: dy += legend.rect.height + opts.labelPadding; break; + case RIGHT: dw += legend.rect.width + opts.labelPadding; break; + case BOTTOM: dh += legend.rect.height; break; + default: + } + offset(dx, dy, dw, dh); + } + + private void draw(Graphics2D g) { + drawPlotArea(g); + drawGrid(g); + drawAxes(g); + drawLegend(g); + // if check needed that content is inside padding + //g.setColor(Color.GRAY); + //g.drawRect(boundRect.x + opts.padding, boundRect.y + opts.padding, boundRect.width - opts.padding * 2, boundRect.height - opts.padding * 2); + } + + private void drawPlotArea(Graphics2D g) { + g.setColor(opts.foregroundColor); + g.drawRect(plotBorderRect.x, plotBorderRect.y, plotBorderRect.width, plotBorderRect.height); + g.setFont(opts.titleFont); + drawLabel(g, opts.title, plotBorderRect.x + toInt(plotBorderRect.getWidth() / 2), opts.padding, HorizAlign.CENTER, VertAlign.TOP); + } + + private void drawGrid(Graphics2D g) { + Stroke stroke = g.getStroke(); + g.setStroke(opts.gridStroke); + g.setColor(opts.gridColor); + + int leftX = plotBorderRect.x + 1; + int rightX = plotBorderRect.x + plotBorderRect.width - 1; + int topY = plotBorderRect.y + 1; + int bottomY = plotBorderRect.y + plotBorderRect.height - 1; + + for (int i = 0; i < opts.grids.x + 1; i++) { + int x = toInt(plotRect.x + (plotRect.getWidth() / opts.grids.x) * i); + g.drawLine(x, topY, x, bottomY); + } + + for (int i = 0; i < opts.grids.y + 1; i++) { + int y = toInt(plotRect.y + (plotRect.getHeight() / opts.grids.y) * i); + g.drawLine(leftX, y, rightX, y); + } + + g.setStroke(stroke); + } + + private void calcAxes(Graphics2D g) { + Axis xAxis = xAxes.isEmpty() ? new Axis("", null) : xAxes.values().iterator().next(); + Axis yAxis = yAxes.isEmpty() ? new Axis("", null) : yAxes.values().iterator().next(); + int xCount = 0, yCount = 0; + for (DataSeries series : dataSeriesMap.values()) { + if (series.opts.xAxis == null) { + series.opts.xAxis = xAxis; + xCount++; + } + if (series.opts.yAxis == null) { + series.opts.yAxis = yAxis; + yCount++; + } + series.addAxesToName(); + } + if (xAxes.isEmpty() && xCount > 0) + xAxes.put("x", xAxis); + if (yAxes.isEmpty() && yCount > 0) + yAxes.put("y", yAxis); + } + + private void calcAxisLabels(Graphics2D g, boolean isX) { + FontMetrics fm = g.getFontMetrics(); + Rectangle2D rect = null; + double w = 0, h = 0; + Map axes = isX ? xAxes : yAxes; + int grids = isX ? opts.grids.x : opts.grids.y; + for (Map.Entry entry : axes.entrySet()) { + Axis axis = entry.getValue(); + axis.labels = new String[grids + 1]; + axis.labelRect = fm.getStringBounds("", g); + double xStep = axis.opts.range.diff / grids; + for (int j = 0; j < grids + 1; j++) { + axis.labels[j] = formatDouble(axis.opts.range.min + xStep * j, axis.opts.format); + rect = fm.getStringBounds(axis.labels[j], g); + if (rect.getWidth() > w) + w = rect.getWidth(); + if (rect.getHeight() > h) + h = rect.getHeight(); + } + axis.labelRect.setRect(0, 0, w, h); + } + } + + private void calcRange(boolean isX) { + for (DataSeries series : dataSeriesMap.values()) { + Axis axis = isX ? series.opts.xAxis : series.opts.yAxis; + if (axis.opts.dynamicRange) { + Range range = isX ? series.xRange() : series.yRange(); + if (axis.opts.range == null) + axis.opts.range = range; + else { + if (range.max > axis.opts.range.max) + axis.opts.range.setMax(range.max); + if (range.min < axis.opts.range.min) + axis.opts.range.setMin(range.min); + } + } + } + Map axes = isX ? xAxes : yAxes; + for (Iterator it = axes.values().iterator(); it.hasNext(); ) { + Axis axis = it.next(); + if (axis.opts.range == null) + it.remove(); + } + } + + private void drawAxes(Graphics2D g) { + g.setFont(opts.labelFont); + g.setColor(opts.foregroundColor); + + int leftXPadded = plotBorderRect.x - opts.labelPadding; + int rightX = plotBorderRect.x + plotBorderRect.width; + int bottomY = plotBorderRect.y + plotBorderRect.height; + int bottomYPadded = bottomY + opts.labelPadding; + + int axisOffset = 0; + for (Map.Entry entry : xAxes.entrySet()) { + Axis axis = entry.getValue(); + double xStep = axis.opts.range.diff / opts.grids.x; + + drawLabel(g, axis.name, rightX + opts.labelPadding, bottomY + axisOffset, HorizAlign.LEFT, VertAlign.CENTER); + g.drawLine(plotRect.x, bottomY + axisOffset, plotRect.x + plotRect.width, bottomY + axisOffset); + + for (int j = 0; j < opts.grids.x + 1; j++) { + int x = toInt(plotRect.x + (plotRect.getWidth() / opts.grids.x) * j); + drawLabel(g, formatDouble(axis.opts.range.min + xStep * j, axis.opts.format), x, bottomYPadded + axisOffset, HorizAlign.CENTER, VertAlign.TOP); + g.drawLine(x, bottomY + axisOffset, x, bottomY + opts.tickSize + axisOffset); + } + axisOffset += toInt(axis.labelRect.getHeight() + opts.labelPadding * 2); + } + + axisOffset = 0; + for (Map.Entry entry : yAxes.entrySet()) { + Axis axis = entry.getValue(); + double yStep = axis.opts.range.diff / opts.grids.y; + + drawLabel(g, axis.name, leftXPadded - axisOffset, plotBorderRect.y - toInt(axis.labelRect.getHeight() + opts.labelPadding), HorizAlign.RIGHT, VertAlign.CENTER); + g.drawLine(plotBorderRect.x - axisOffset, plotRect.y + plotRect.height, plotBorderRect.x - axisOffset, plotRect.y); + + for (int j = 0; j < opts.grids.y + 1; j++) { + int y = toInt(plotRect.y + (plotRect.getHeight() / opts.grids.y) * j); + drawLabel(g, formatDouble(axis.opts.range.max - yStep * j, axis.opts.format), leftXPadded - axisOffset, y, HorizAlign.RIGHT, VertAlign.CENTER); + g.drawLine(plotBorderRect.x - axisOffset, y, plotBorderRect.x - opts.tickSize - axisOffset, y); + } + axisOffset += toInt(axis.labelRect.getWidth() + opts.labelPadding * 2); + } + } + + private void calcLegend(Graphics2D g) { + legend.rect = new Rectangle(0, 0); + if (opts.legend == LegendFormat.NONE) + return; + int size = dataSeriesMap.size(); + if (size == 0) + return; + + FontMetrics fm = g.getFontMetrics(); + Iterator it = dataSeriesMap.values().iterator(); + legend.labelRect = fm.getStringBounds(it.next().nameWithAxes, g); + int legendSignSize = opts.defaultLegendSignSize; + while (it.hasNext()) { + DataSeries series = it.next(); + Rectangle2D rect = fm.getStringBounds(series.nameWithAxes, g); + if (rect.getWidth() > legend.labelRect.getWidth()) + legend.labelRect.setRect(0, 0, rect.getWidth(), legend.labelRect.getHeight()); + if (rect.getHeight() > legend.labelRect.getHeight()) + legend.labelRect.setRect(0, 0, legend.labelRect.getWidth(), rect.getHeight()); + switch (series.opts.marker) { + case CIRCLE: case SQUARE: + if (series.opts.markerSize + opts.defaultLegendSignSize > legendSignSize) + legendSignSize = series.opts.markerSize + opts.defaultLegendSignSize; + break; + case DIAMOND: + if (series.getDiagMarkerSize() + opts.defaultLegendSignSize > legendSignSize) + legendSignSize = series.getDiagMarkerSize() + opts.defaultLegendSignSize; + break; + default: + } + } + opts.legendSignSize = legendSignSize; + + legend.entryWidth = legendSignSize + opts.labelPadding + toInt(legend.labelRect.getWidth()); + legend.entryWidthPadded = legend.entryWidth + opts.labelPadding; + + switch (opts.legend) { + case TOP: case BOTTOM: + legend.entryCount = (int) Math.floor((double) (plotBorderRect.width - opts.labelPadding) / legend.entryWidthPadded); + legend.xCount = size <= legend.entryCount ? size : legend.entryCount; + legend.yCount = size <= legend.entryCount ? 1 : (int) Math.ceil((double) size / legend.entryCount); + legend.rect.width = opts.labelPadding + (legend.xCount * legend.entryWidthPadded); + legend.rect.height = opts.labelPadding + toInt(legend.yCount * (opts.labelPadding + legend.labelRect.getHeight())); + legend.rect.x = plotBorderRect.x + (plotBorderRect.width - legend.rect.width) / 2; + if (opts.legend == LegendFormat.TOP) + legend.rect.y = plotBorderRect.y; + else + legend.rect.y = boundRect.height - legend.rect.height - opts.padding; + break; + case RIGHT: + legend.rect.width = opts.labelPadding * 3 + legendSignSize + toInt(legend.labelRect.getWidth()); + legend.rect.height = opts.labelPadding * (size + 1) + toInt(legend.labelRect.getHeight() * size); + legend.rect.x = boundRect.width - legend.rect.width - opts.padding; + legend.rect.y = plotBorderRect.y + plotBorderRect.height / 2 - legend.rect.height / 2; + break; + default: + } + } + + private void drawLegend(Graphics2D g) { + if (opts.legend == LegendFormat.NONE) + return; + + g.drawRect(legend.rect.x, legend.rect.y, legend.rect.width, legend.rect.height); + int labelHeight = toInt(legend.labelRect.getHeight()); + int x = legend.rect.x + opts.labelPadding; + int y = legend.rect.y + opts.labelPadding + labelHeight / 2; + + switch (opts.legend) { + case TOP: case BOTTOM: + int i = 0; + for (DataSeries series : dataSeriesMap.values()) { + drawLegendEntry(g, series, x, y); + x += legend.entryWidthPadded; + if ((i + 1) % legend.xCount == 0) { + x = legend.rect.x + opts.labelPadding; + y += opts.labelPadding + labelHeight; + } + i++; + } + break; + case RIGHT: + for (DataSeries series : dataSeriesMap.values()) { + drawLegendEntry(g, series, x, y); + y += opts.labelPadding + labelHeight; + } + break; + default: + } + } + + private void drawLegendEntry(Graphics2D g, DataSeries series, int x, int y) { + series.fillArea(g, x, y, x + opts.legendSignSize, y, y + opts.legendSignSize / 2); + series.drawLine(g, x, y, x + opts.legendSignSize, y); + series.drawMarker(g, x + opts.legendSignSize / 2, y, x, y + opts.legendSignSize / 2); + g.setColor(opts.foregroundColor); + drawLabel(g, series.nameWithAxes, x + opts.legendSignSize + opts.labelPadding, y, HorizAlign.LEFT, VertAlign.CENTER); + } + + } + + public static class Range { + + private double min; + private double max; + private double diff; + + public Range(double min, double max) { + this.min = min; + this.max = max; + this.diff = max - min; + } + + public Range(Range range) { + this.min = range.min; + this.max = range.max; + this.diff = max - min; + } + + public void setMin(double min) { + this.min = min; + this.diff = max - min; + } + + public void setMax(double max) { + this.max = max; + this.diff = max - min; + } + + @Override + public String toString() { + return "Range [min=" + min + ", max=" + max + "]"; + } + + } + + public static AxisOptions axisOpts() { + return new AxisOptions(); + } + + public static class AxisOptions { + + private AxisFormat format = AxisFormat.NUMBER; + private boolean dynamicRange = true; + private Range range; + + public AxisOptions format(AxisFormat format) { + this.format = format; + return this; + } + + public AxisOptions range(double min, double max) { + this.range = new Range(min, max); + this.dynamicRange = false; + return this; + } + + } + + private class Axis { + + private String name; + private AxisOptions opts = new AxisOptions(); + private Rectangle2D labelRect; + private String[] labels; + + public Axis(String name, AxisOptions opts) { + this.name = name; + if (opts != null) + this.opts = opts; + } + + @Override + public String toString() { + return "Axis [name=" + name + ", opts=" + opts + "]"; + } + + } + + public static DataSeriesOptions seriesOpts() { + return new DataSeriesOptions(); + } + + public static class DataSeriesOptions { + + private Color seriesColor = Color.BLUE; + private Line line = Line.SOLID; + private int lineWidth = 2; + private float[] lineDash = new float[] { 3.0f, 3.0f }; + private Marker marker = Marker.NONE; + private int markerSize = 10; + private Color markerColor = Color.WHITE; + private Color areaColor = null; + private String xAxisName; + private String yAxisName; + private Axis xAxis; + private Axis yAxis; + + public DataSeriesOptions color(Color seriesColor) { + this.seriesColor = seriesColor; + return this; + } + + public DataSeriesOptions line(Line line) { + this.line = line; + return this; + } + + public DataSeriesOptions lineWidth(int width) { + this.lineWidth = width; + return this; + } + + public DataSeriesOptions lineDash(float[] dash) { + this.lineDash = dash; + return this; + } + + public DataSeriesOptions marker(Marker marker) { + this.marker = marker; + return this; + } + + public DataSeriesOptions markerSize(int markerSize) { + this.markerSize = markerSize; + return this; + } + + public DataSeriesOptions markerColor(Color color) { + this.markerColor = color; + return this; + } + + public DataSeriesOptions areaColor(Color color) { + this.areaColor = color; + return this; + } + + public DataSeriesOptions xAxis(String name) { + this.xAxisName = name; + return this; + } + + public DataSeriesOptions yAxis(String name) { + this.yAxisName = name; + return this; + } + + private void setPlot(Plot plot) { + if (plot != null) + this.xAxis = plot.xAxes.get(xAxisName); + if (plot != null) + this.yAxis = plot.yAxes.get(yAxisName); + } + + } + + public static Data data() { + return new Data(); + } + + public static class Data { + + private double[] x1; + private double[] y1; + private List x2; + private List y2; + + private Data() {} + + public Data xy(double[] x, double[] y) { + this.x1 = x; + this.y1 = y; + return this; + } + + public Data xy(double x, double y) { + if (this.x2 == null || this.y2 == null) { + this.x2 = new ArrayList(10); + this.y2 = new ArrayList(10); + } + x2.add(x); + y2.add(y); + return this; + } + + public Data xy(List x, List y) { + this.x2 = x; + this.y2 = y; + return this; + } + + public int size() { + if (x1 != null) + return x1.length; + if (x2 != null) + return x2.size(); + return 0; + } + + public double x(int i) { + if (x1 != null) + return x1[i]; + if (x2 != null) + return x2.get(i); + return 0; + } + + public double y(int i) { + if (y1 != null) + return y1[i]; + if (y2 != null) + return y2.get(i); + return 0; + } + + } + + public class DataSeries { + + private String name; + private String nameWithAxes; + private DataSeriesOptions opts = new DataSeriesOptions(); + private Data data; + + public DataSeries(String name, Data data, DataSeriesOptions opts) { + if (opts != null) + this.opts = opts; + this.name = name; + this.data = data; + if (this.data == null) + this.data = data(); + } + + public void clear() { + } + + private void addAxesToName() { + this.nameWithAxes = this.name + " (" + opts.yAxis.name + "/" + opts.xAxis.name + ")"; + } + + private Range xRange() { + Range range = new Range(0, 0); + if (data != null && data.size() > 0) { + range = new Range(data.x(0), data.x(0)); + for (int i = 1; i < data.size(); i++) { + if (data.x(i) > range.max) + range.setMax(data.x(i)); + if (data.x(i) < range.min) + range.setMin(data.x(i)); + } + } + return range; + } + + private Range yRange() { + Range range = new Range(0, 0); + if (data != null && data.size() > 0) { + range = new Range(data.y(0), data.y(0)); + for (int i = 1; i < data.size(); i++) { + if (data.y(i) > range.max) + range.setMax(data.y(i)); + if (data.y(i) < range.min) + range.setMin(data.y(i)); + } + } + return range; + } + + private void draw(Graphics2D g) { + g.setClip(plotArea.plotClipRect); + if (data != null) { + double x1 = 0, y1 = 0; + int size = data.size(); + if (opts.line != Line.NONE) + for (int j = 0; j < size; j++) { + double x2 = x2x(data.x(j), opts.xAxis.opts.range, plotArea.xPlotRange); + double y2 = y2y(data.y(j), opts.yAxis.opts.range, plotArea.yPlotRange); + int ix1 = toInt(x1), iy1 = toInt(y1), ix2 = toInt(x2), iy2 = toInt(y2); + int iy3 = plotArea.plotRect.y + plotArea.plotRect.height; + // special case for the case when only the first point present + if (size == 1) { + ix1 = ix2; + iy1 = iy2; + } + if (j != 0 || size == 1) { + fillArea(g, ix1, iy1, ix2, iy2, iy3); + drawLine(g, ix1, iy1, ix2, iy2); + } + x1 = x2; + y1 = y2; + } + + int halfMarkerSize = opts.markerSize / 2; + int halfDiagMarkerSize = getDiagMarkerSize() / 2; + g.setStroke(new BasicStroke(2)); + if (opts.marker != Marker.NONE) + for (int j = 0; j < size; j++) { + double x2 = x2x(data.x(j), opts.xAxis.opts.range, plotArea.xPlotRange); + double y2 = y2y(data.y(j), opts.yAxis.opts.range, plotArea.yPlotRange); + drawMarker(g, halfMarkerSize, halfDiagMarkerSize, x2, y2, + plotArea.plotRect.x, plotArea.plotRect.y + plotArea.plotRect.height); + } + } + } + + private int getDiagMarkerSize() { + return (int) Math.round(Math.sqrt(2 * opts.markerSize * opts.markerSize)); + } + + private void fillArea(Graphics2D g, int ix1, int iy1, int ix2, int iy2, int iy3) { + if (opts.areaColor != null) { + g.setColor(opts.areaColor); + g.fill(new Polygon( + new int[] { ix1, ix2, ix2, ix1 }, + new int[] { iy1, iy2, iy3, iy3 }, + 4)); + g.setColor(opts.seriesColor); + } + } + + private void drawLine(Graphics2D g, int ix1, int iy1, int ix2, int iy2) { + if (opts.line != Line.NONE) { + g.setColor(opts.seriesColor); + setStroke(g); + g.drawLine(ix1, iy1, ix2, iy2); + } + } + + private void setStroke(Graphics2D g) { + switch (opts.line) { + case SOLID: + g.setStroke(new BasicStroke(opts.lineWidth)); + break; + case DASHED: + g.setStroke(new BasicStroke(opts.lineWidth, BasicStroke.CAP_ROUND, + BasicStroke.JOIN_ROUND, 10.0f, opts.lineDash, 0.0f)); + break; + default: + } + } + + private void drawMarker(Graphics2D g, int x2, int y2, int x3, int y3) { + int halfMarkerSize = opts.markerSize / 2; + int halfDiagMarkerSize = getDiagMarkerSize() / 2; + g.setStroke(new BasicStroke(2)); + drawMarker(g, halfMarkerSize, halfDiagMarkerSize, x2, y2, x3, y3); + } + + private void drawMarker(Graphics2D g, int halfMarkerSize, int halfDiagMarkerSize, double x2, double y2, double x3, double y3) { + switch (opts.marker) { + case CIRCLE: + g.setColor(opts.markerColor); + g.fillOval(toInt(x2 - halfMarkerSize), toInt(y2 - halfMarkerSize), opts.markerSize, opts.markerSize); + g.setColor(opts.seriesColor); + g.drawOval(toInt(x2 - halfMarkerSize), toInt(y2 - halfMarkerSize), opts.markerSize, opts.markerSize); + break; + case SQUARE: + g.setColor(opts.markerColor); + g.fillRect(toInt(x2 - halfMarkerSize), toInt(y2 - halfMarkerSize), opts.markerSize, opts.markerSize); + g.setColor(opts.seriesColor); + g.drawRect(toInt(x2 - halfMarkerSize), toInt(y2 - halfMarkerSize), opts.markerSize, opts.markerSize); + break; + case DIAMOND: + int[] xpts = { toInt(x2), toInt(x2 + halfDiagMarkerSize), toInt(x2), toInt(x2 - halfDiagMarkerSize) }; + int[] ypts = { toInt(y2 - halfDiagMarkerSize), toInt(y2), toInt(y2 + halfDiagMarkerSize), toInt(y2) }; + g.setColor(opts.markerColor); + g.fillPolygon(xpts, ypts, 4); + g.setColor(opts.seriesColor); + g.drawPolygon(xpts, ypts, 4); + break; + case COLUMN: + g.setColor(opts.markerColor); + g.fillRect(toInt(x2), toInt(y2), opts.markerSize, toInt(y3 - y2)); + g.setColor(opts.seriesColor); + g.drawRect(toInt(x2), toInt(y2), opts.markerSize, toInt(y3 - y2)); + break; + case BAR: + g.setColor(opts.markerColor); + g.fillRect(toInt(x3), toInt(y2), toInt(x2 - x3), opts.markerSize); + g.setColor(opts.seriesColor); + g.drawRect(toInt(x3), toInt(y2), toInt(x2 - x3), opts.markerSize); + break; + default: + } + } + + } + + private static void drawLabel(Graphics2D g, String s, int x, int y, HorizAlign hAlign, VertAlign vAlign) { + FontMetrics fm = g.getFontMetrics(); + Rectangle2D rect = fm.getStringBounds(s, g); + + // by default align by left + if (hAlign == HorizAlign.RIGHT) + x -= rect.getWidth(); + else if (hAlign == HorizAlign.CENTER) + x -= rect.getWidth() / 2; + + // by default align by bottom + if (vAlign == VertAlign.TOP) + y += rect.getHeight(); + else if (vAlign == VertAlign.CENTER) + y += rect.getHeight() / 2; + + g.drawString(s, x, y); + } + + public static String formatDouble(double d, AxisFormat format) { + switch (format) { + case TIME_HM: return String.format("%tR", new java.util.Date((long) d)); + case TIME_HMS: return String.format("%tT", new java.util.Date((long) d)); + case DATE: return String.format("%tF", new java.util.Date((long) d)); + case DATETIME_HM: return String.format("%tF %1$tR", new java.util.Date((long) d)); + case DATETIME_HMS: return String.format("%tF %1$tT", new java.util.Date((long) d)); + case NUMBER_KGM: return formatDoubleAsNumber(d, true); + case NUMBER_INT: return Integer.toString((int) d); + default: return formatDoubleAsNumber(d, false); + } + } + + private static String formatDoubleAsNumber(double d, boolean useKGM) { + if (useKGM && d > 1000 && d < 1000000000000l) { + long[] numbers = new long[] { 1000l, 1000000l, 1000000000l }; + char[] suffix = new char[] { 'K', 'M', 'G' }; + + int i = 0; + double r = 0; + for (long number : numbers) { + r = d / number; + if (r < 1000) + break; + i++; + } + if (i == suffix.length) + i--; + return String.format("%1$,.2f%2$c", r, suffix[i]); + } + else + return String.format("%1$.3G", d); + } + + private static double x2x(double x, Range xr1, Range xr2) { + return xr1.diff == 0 ? xr2.min + xr2.diff / 2 : xr2.min + (x - xr1.min) / xr1.diff * xr2.diff; + } + + // y axis is reverse in Graphics + private static double y2y(double x, Range xr1, Range xr2) { + return xr1.diff == 0 ? xr2.min + xr2.diff / 2 : xr2.max - (x - xr1.min) / xr1.diff * xr2.diff; + } + + private static int toInt(double d) { + return (int) Math.round(d); + } + +} diff --git a/java/src/main/java/utilities/Utils.java b/java/src/main/java/utilities/Utils.java index f2b85ab..190197e 100644 --- a/java/src/main/java/utilities/Utils.java +++ b/java/src/main/java/utilities/Utils.java @@ -5,12 +5,10 @@ import bpmEstimation.AccelerometerWindowBuffer; import bpmEstimation.Peaks; import javax.swing.*; import java.awt.*; +import java.awt.event.WindowEvent; import java.awt.image.BufferedImage; -import java.text.Collator; -import java.util.Arrays; -import java.util.Collections; -import java.util.Comparator; -import java.util.LinkedList; +import java.util.*; +import java.util.List; import java.util.stream.IntStream; //TODO: change from double to generic type @@ -22,6 +20,38 @@ public class Utils { return (double) Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); } + public static double getDistance(double x1, double y1, double z1, double x2, double y2, double z2) { + return (double) Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2)); + } + + public static void reverse(double[] array) { + if (array == null) { + return; + } + int i = 0; + int j = array.length - 1; + double tmp; + while (j > i) { + tmp = array[j]; + array[j] = array[i]; + array[i] = tmp; + j--; + i++; + } + } + + public static double var(LinkedList data){ + double mean = mean(data); + double temp = 0; + for(double a : data) + temp += (a-mean)*(a-mean); + return temp/(data.size()-1); + } + + public static double stdDev(LinkedList data){ + return Math.sqrt(var(data)); + } + public static double sqr(double x) { return x * x; } @@ -94,7 +124,7 @@ public class Utils { return sum(data) / data.size(); } - public static double median(LinkedList data){ + public static double median(List data){ data.sort(Comparator.naturalOrder()); double median; @@ -106,6 +136,18 @@ public class Utils { return median; } + public static double mad(List data){ + + double median = median(data); + + java.util.List tmpList = new ArrayList<>(); + for(double value : data){ + tmpList.add(Math.abs(value - median)); + } + + return median(tmpList); + } + //TODO: Could be slow.. faster method? public static double geometricMean(double[] data) { double sum = data[0]; @@ -115,6 +157,17 @@ public class Utils { return Math.pow(sum, 1.0 / data.length); } + public static double geometricMeanLog(double[] data){ + double GM_log = 0.0d; + for (int i = 0; i < data.length; ++i) { + if (data[i] == 0.0d) { + return 0.0d; + } + GM_log += Math.log(data[i]); + } + return Math.exp(GM_log / data.length); + } + public static int intersectionNumber(double[] signal, double border){ int cnt = 0; boolean isSmallerValue = false; @@ -170,6 +223,20 @@ public class Utils { return newArray; } + public static void removeOutliersZScore(List data, double score) { + + if(!data.isEmpty()){ + double median = median(data); + double mad = mad(data); + + for(Iterator it = data.iterator(); it.hasNext(); ){ + + double M = Math.abs((0.6745 * (it.next() - median)) / mad); + if (M > score){ it.remove(); } + } + } + } + public static double magnitude(double x, double y, double z){ return Math.sqrt((x*x) + (y*y) + (z*z)); } @@ -287,6 +354,7 @@ public class Utils { private Utils.ShowPNG plotCorrY; private Utils.ShowPNG plotCorrZ; private Utils.ShowPNG plotCorrMag; + private Utils.ShowPNG plotCorr3D; public DebugPlotter(){ @@ -302,6 +370,7 @@ public class Utils { plotCorrY = new Utils.ShowPNG(); plotCorrZ = new Utils.ShowPNG(); plotCorrMag = new Utils.ShowPNG(); + plotCorr3D = new Utils.ShowPNG(); } public void setPlotRawX(long[] ts, double[] data){ @@ -316,8 +385,7 @@ public class Utils { plotRawZ.plotData("Raw Data Z", "z", ts, data); } - public void setPlotRawMag(long[] ts, double[] data){ - plotRawMag.plotData("Raw Data Magnitude", "mag", ts, data); + public void setPlotRawMag(long[] ts, double[] data){plotRawMag.plotData("Raw Data Magnitude", "mag", ts, data); } public void setPlotButterX(long[] ts, double[] data){ @@ -325,11 +393,11 @@ public class Utils { } public void setPlotButterY(long[] ts, double[] data){ - plotButterY.plotData("Butter Data Y", "z", ts, data); + plotButterY.plotData("Butter Data Y", "y", ts, data); } public void setPlotButterZ(long[] ts, double[] data){ - plotButterZ.plotData("Butter Data Z", "y", ts, data); + plotButterZ.plotData("Butter Data Z", "z", ts, data); } public void setPlotButterMag(long[] ts, double[] data){ @@ -352,6 +420,26 @@ public class Utils { plotCorrMag.plotCorrWithPeaks("Autocorr Mag", "mag", data,2, peaks, true); } + public void setPlotCorr3D(double[] data, Peaks peaks){ + plotCorr3D.plotCorrWithPeaks("DistCorr - 3D", "3D", data,2, peaks, true); + }; + + public void close(){ + plotRawX.dispatchEvent(new WindowEvent(plotRawX, WindowEvent.WINDOW_CLOSING)); + plotRawY.dispatchEvent(new WindowEvent(plotRawY, WindowEvent.WINDOW_CLOSING)); + plotRawZ.dispatchEvent(new WindowEvent(plotRawZ, WindowEvent.WINDOW_CLOSING)); + plotRawMag.dispatchEvent(new WindowEvent(plotRawMag, WindowEvent.WINDOW_CLOSING)); + plotButterX.dispatchEvent(new WindowEvent(plotButterX, WindowEvent.WINDOW_CLOSING)); + plotButterY.dispatchEvent(new WindowEvent(plotButterY, WindowEvent.WINDOW_CLOSING)); + plotButterZ.dispatchEvent(new WindowEvent(plotButterZ, WindowEvent.WINDOW_CLOSING)); + plotButterMag.dispatchEvent(new WindowEvent(plotButterMag, WindowEvent.WINDOW_CLOSING)); + plotCorrX.dispatchEvent(new WindowEvent(plotCorrX, WindowEvent.WINDOW_CLOSING)); + plotCorrY.dispatchEvent(new WindowEvent(plotCorrY, WindowEvent.WINDOW_CLOSING)); + plotCorrZ.dispatchEvent(new WindowEvent(plotCorrZ, WindowEvent.WINDOW_CLOSING)); + plotCorrMag.dispatchEvent(new WindowEvent(plotCorrMag, WindowEvent.WINDOW_CLOSING)); + plotCorr3D.dispatchEvent(new WindowEvent(plotCorr3D, WindowEvent.WINDOW_CLOSING)); + + } } diff --git a/matlab/AutoCorrMethodNew_Watch.m b/matlab/AutoCorrMethodNew_Watch.m index b13ee1b..f5002b5 100644 --- a/matlab/AutoCorrMethodNew_Watch.m +++ b/matlab/AutoCorrMethodNew_Watch.m @@ -1,58 +1,47 @@ -%We are using a threshold-based version for bpm estimation -%Only the z axis of the acc is used +%using autocorrelation to estimate the current bmp within some fixed window -%NOTE: depending on the measurement device we have a highly different sample rate. the smartwatches are not capable of providing a constant sample rate. the xsens on the other hand is able to do this. - -%load file provided by the sensor readout app - -% SMARTWATCH LG WEAR ------> 100 hz - 1000hz -%measurements = dlmread('../../measurements/lgWear/PR_recording_80bpm_4-4_177596720.csv', ';'); % -%measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176527527.csv', ';'); -%measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176606785.csv', ';'); -%measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176696356.csv', ';'); -%measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176820066.csv', ';'); % -%measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176931941.csv', ';'); %double -%measurements = dlmread('../../measurements/lgWear/recording_72bpm_4-4_176381633.csv', ';'); -%measurements = dlmread('../../measurements/lgWear/recording_72bpm_4-4_176453327.csv', ';'); % -%measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_176073767.csv', ';'); %* -%measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_176165357.csv', ';'); -%measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_176230146.csv', ';'); -%measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_176284687.csv', ';'); % -%measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_177368860.csv', ';'); %(besonders) -%measurements = dlmread('../../measurements/lgWear/recording_180bpm_4-4_177011641.csv', ';'); % -%measurements = dlmread('../../measurements/lgWear/recording_180bpm_4-4_177064915.csv', ';'); % - - -% SMARTWATCH G WATCH WEAR R ----> 100hz - 250hz -%measurements = dlmread('../measurements/wearR/PR_recording_80bpm_4-4_177596720.csv', ';'); %* -%measurements = dlmread('../measurements/wearR/recording_48bpm_4-4_176527527.csv', ';'); -%measurements = dlmread('../measurements/wearR/recording_48bpm_4-4_176606785.csv', ';'); -%measurements = dlmread('../../measurements/wearR/recording_48bpm_4-4_176696356.csv', ';'); %* -%measurements = dlmread('../measurements/wearR/recording_48bpm_4-4_176820066.csv', ';'); -%measurements = dlmread('../measurements/wearR/recording_48bpm_4-4_176931941.csv', ';'); %double -%measurements = dlmread('../measurements/wearR/recording_72bpm_4-4_176381633.csv', ';'); -%measurements = dlmread('../measurements/wearR/recording_72bpm_4-4_176453327.csv', ';'); -%measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_176073767.csv', ';'); * -%measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_176165357.csv', ';'); -%measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_176230146.csv', ';'); %* 72? -%measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_176284687.csv', ';'); %*48? -%measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_177368860.csv', ';'); %(besonders) -%measurements = dlmread('../measurements/wearR/recording_180bpm_4-4_177011641.csv', ';'); * -%measurements = dlmread('../measurements/wearR/recording_180bpm_4-4_177064915.csv', ';'); * - -files = dir(fullfile('../../measurements/mSensor/', '*.csv')); -%files = dir(fullfile('../../measurements/lgWear/', '*.csv')); +%load sensor files +%files = dir(fullfile('../../measurements/2017.06/mSensor/', '*.csv')); +%files = natsortfiles(dir(fullfile('../../measurements/2017.06/lgWear/', '*.csv'))); %files = dir(fullfile('../../measurements/wearR/', '*.csv')); %files = dir(fullfile('../../measurements/peter_failed/', '*.csv')); +%files = dir(fullfile('../../measurements/2018.06/manfred/LGWatchR/', '*.csv')); +%files = dir(fullfile('../../measurements/2018.06/peter/Huawai/', '*.csv')); +%files = dir(fullfile('../../measurements/2018.06/peter/mSensor/', '*.csv')); +%files = dir(fullfile('../../measurements/2018.06/frank/mSensor/', '*.csv')); +files = dir(fullfile('../../measurements/2018.06/leon/mSensor/', '*.csv')); +%files_sorted = natsortfiles({files.name}); for file = files' filename = [file.folder '/' file.name]; - measurements = dlmread(filename, ';'); + measurements = dlmread(filename, ';', 3, 0); + + %load ground truth file + fid = fopen(filename); + fgetl(fid); + Str = fgetl(fid); + Key = 'Metronom: '; + Index = strfind(Str, Key); + gtDataRaw = sscanf(Str(Index(1) + length(Key):end), '%g', 1); + gtData = []; + gtFile = []; + if(isempty(gtDataRaw)) + gtFile = extractAfter(Str, Key); + gtFile = strcat('../../measurements/2018.06/gt_toni/', gtFile); + f = fopen(gtFile); + gtDataRaw = textscan(f, '%f %s', 'Delimiter', ' '); + fclose(f); + [~,~,~,hours,minutes,seconds] = datevec(gtDataRaw{2}, 'HH:MM:SS.FFF'); + gtData(:,1) = 1000*(60*minutes + seconds); %we do not use hours! + gtData(:,2) = gtDataRaw{1}; + else + gtData = gtDataRaw; + end %draw the raw acc data m_idx = []; - m_idx = (measurements(:,2)==3); %Android App: 10, Normal Data: 2 + m_idx = (measurements(:,2)==3); %Android App: 10, Sensor: 3, Normal Data: 2 m = measurements(m_idx, :); %Interpolate to generate a constant sample rate to 250hz (4ms per sample) @@ -66,53 +55,74 @@ for file = files' %put all together again m = [t_interp', t_interp', m_interp]; - figure(1); - plot(m(:,1),m(:,3)) %x - legend("x", "location", "eastoutside"); - - figure(2); - plot(m(:,1),m(:,4)) %yt - legend("y", "location", "eastoutside"); - - figure(3); - plot(m(:,1),m(:,5)) %z - legend("z", "location", "eastoutside"); +% figure(1); +% plot(m(:,1),m(:,3)) %x +% legend("x", "location", "eastoutside"); +% +% figure(2); +% plot(m(:,1),m(:,4)) %yt +% legend("y", "location", "eastoutside"); +% +% figure(3); +% plot(m(:,1),m(:,5)) %z +% legend("z", "location", "eastoutside"); +% +% %magnitude + magnitude = sqrt(sum(m(:,3:5).^2,2)); +% figure(5); +% plot(m(:,1), magnitude); +% legend("magnitude", "location", "eastoutside"); - %magnitude - magnitude = sqrt(sum(m(:,3:5).^2,2)); - figure(5); - plot(m(:,1), magnitude); - legend("magnitude", "location", "eastoutside"); - - waitforbuttonpress(); + %waitforbuttonpress(); %save timestamps timestamps = m(:,1); data = m(:,3); %only z %TODO: Different window sizes for periods under 16.3 s - window_size = 2048; %about 2 seconds using 2000hz, 16.3 s using 250hz + window_size = 1024; %about 2 seconds using 2000hz, 16.3 s using 250hz overlap = 256; bpm_per_window_ms = []; bpm_per_window = []; - for i = window_size+1:length(data) + bpm_3D = []; + ms_3D = []; + + gtIdx = 1; + gtError_3D = []; + gtError_1D = []; + + for i = window_size+1:1:length(data) %wait until window is filled with new data if(mod(i,overlap) == 0) - + + %set cur ground truth + if(length(gtData) > 1) + curTimestamp = timestamps(i); + while(curTimestamp > gtData(gtIdx,1) && gtIdx < length(gtData)) + curGtBpm = gtData(gtIdx,2); + gtIdx = gtIdx + 1; + end + else + curGtBpm = gtData; + end %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_mag, lag_mag] = xcov(magnitude(i-window_size:i), (window_size/4), "coeff"); - - - %autocorrelation of the autocorrelation?! - %[corr_corr_x, lag_lag_x] = xcov(corr_x, length(corr_x), "coeff"); - %[corr_corr_y, lag_lag_y] = xcov(corr_y, length(corr_x), "coeff"); - %[corr_corr_z, lag_lag_z] = xcov(corr_z, length(corr_x), "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"); + + %magnitude + [corr_mag, lag_mag] = xcov(magnitude(i-window_size:i), (window_size/2), "coeff"); + %TODO: stichwort spatial autocorrelation + %figure(77); + %scatter3(timestamps(i-window_size:i), m(i-window_size:i,4), m(i-window_size:i,5)); + %distanz zwischen den vektoren nehmen und in eine normale autocorrelation zu packen + %aufpassen wegen der norm, dass die richtung quasi nicht verloren geht. + %https://en.wikipedia.org/wiki/Lp_space + [corr_3D, lag_3D] = distCorr(m(i-window_size:i, 3:5), (round(window_size * 0.8))); + corr_x_pos = corr_x; corr_y_pos = corr_y; corr_z_pos = corr_z; @@ -127,67 +137,81 @@ for file = files' [peak_y, idx_y_raw] = findpeaks(corr_y_pos, 'MinPeakHeight', 0.1,'MinPeakDistance', 50, 'MinPeakProminence', 0.1); [peak_z, idx_z_raw] = findpeaks(corr_z_pos, 'MinPeakHeight', 0.1,'MinPeakDistance', 50, 'MinPeakProminence', 0.1); [peak_mag, idx_mag_raw] = findpeaks(corr_mag_pos, 'MinPeakHeight', 0.1,'MinPeakDistance', 50, 'MinPeakProminence', 0.1); + [peak_3D, idx_3D_raw] = findpeaks(corr_3D, 'MinPeakHeight', 0.1,'MinPeakDistance', 50, 'MinPeakProminence', 0.1); + idx_x_raw = sort(idx_x_raw); idx_y_raw = sort(idx_y_raw); idx_z_raw = sort(idx_z_raw); idx_mag_raw = sort(idx_mag_raw); + idx_3D_raw = sort(idx_3D_raw); idx_x = findFalseDetectedPeaks(idx_x_raw, lag_x, corr_x); idx_y = findFalseDetectedPeaks(idx_y_raw, lag_y, corr_y); idx_z = findFalseDetectedPeaks(idx_z_raw, lag_z, corr_z); idx_mag = findFalseDetectedPeaks(idx_mag_raw, lag_mag, corr_mag); + idx_3D = findFalseDetectedPeaks(idx_3D_raw, lag_3D', corr_3D); + + Dwindow = m(i-window_size:i,3); + Dwindow_mean_ts_diff = mean(diff(lag_3D(idx_3D) * sample_rate_ms)); %2.5 ms is the time between two samples at 400hz + Dwindow_mean_bpm = (60000 / (Dwindow_mean_ts_diff)); + + figure(10); + plot(lag_3D, corr_3D, lag_3D(idx_3D), corr_3D(idx_3D), 'r*', lag_3D(idx_3D_raw), corr_3D(idx_3D_raw), 'g*') + hold ("on") + m_label_ms = strcat(" mean ms: ", num2str(Dwindow_mean_ts_diff)); + m_label_bpm = strcat(" mean bpm: ", num2str(Dwindow_mean_bpm)); + title(strcat(" ", m_label_ms, " ", m_label_bpm)); + hold ("off"); Xwindow = m(i-window_size:i,3); Xwindow_mean_ts_diff = mean(diff(lag_x(idx_x) * sample_rate_ms)); %2.5 ms is the time between two samples at 400hz Xwindow_mean_bpm = (60000 / (Xwindow_mean_ts_diff)); - - figure(11); - plot(lag_x, corr_x, lag_x(idx_x), corr_x(idx_x), 'r*', lag_x(idx_x_raw), corr_x(idx_x_raw), 'g*') %z - hold ("on") - m_label_ms = strcat(" mean ms: ", num2str(Xwindow_mean_ts_diff)); - m_label_bpm = strcat(" mean bpm: ", num2str(Xwindow_mean_bpm)); - title(strcat(" ", m_label_ms, " ", m_label_bpm)); - hold ("off"); + +% figure(11); +% plot(lag_x, corr_x, lag_x(idx_x), corr_x(idx_x), 'r*', lag_x(idx_x_raw), corr_x(idx_x_raw), 'g*') %z +% hold ("on") +% m_label_ms = strcat(" mean ms: ", num2str(Xwindow_mean_ts_diff)); +% m_label_bpm = strcat(" mean bpm: ", num2str(Xwindow_mean_bpm)); +% title(strcat(" ", m_label_ms, " ", m_label_bpm)); +% hold ("off"); Ywindow = m(i-window_size:i,4); Ywindow_mean_ts_diff = mean(diff(lag_y(idx_y) * sample_rate_ms)); Ywindow_mean_bpm = (60000 / (Ywindow_mean_ts_diff)); - figure(12); - plot(lag_y, corr_y, lag_y(idx_y), corr_y(idx_y), 'r*', lag_y(idx_y_raw), corr_y(idx_y_raw), 'g*') %z - hold ("on") - m_label_ms = strcat(" mean ms: ", num2str(Ywindow_mean_ts_diff)); - m_label_bpm = strcat(" mean bpm: ", num2str(Ywindow_mean_bpm)); - title(strcat(" ", m_label_ms, " ", m_label_bpm)); - hold ("off"); +% figure(12); +% plot(lag_y, corr_y, lag_y(idx_y), corr_y(idx_y), 'r*', lag_y(idx_y_raw), corr_y(idx_y_raw), 'g*') %z +% hold ("on") +% m_label_ms = strcat(" mean ms: ", num2str(Ywindow_mean_ts_diff)); +% m_label_bpm = strcat(" mean bpm: ", num2str(Ywindow_mean_bpm)); +% title(strcat(" ", m_label_ms, " ", m_label_bpm)); +% hold ("off"); Zwindow = m(i-window_size:i,5); Zwindow_mean_ts_diff = mean(diff(lag_z(idx_z)* sample_rate_ms)); Zwindow_mean_bpm = (60000 / (Zwindow_mean_ts_diff)); - figure(13); - plot(lag_z, corr_z, lag_z(idx_z), corr_z(idx_z), 'r*', lag_z(idx_z_raw), corr_z(idx_z_raw), 'g*') %z - hold ("on") - m_label_ms = strcat(" mean ms: ", num2str(Zwindow_mean_ts_diff)); - m_label_bpm = strcat(" mean bpm: ", num2str(Zwindow_mean_bpm)); - title(strcat(" ", m_label_ms, " ", m_label_bpm)); - hold ("off"); +% figure(13); +% plot(lag_z, corr_z, lag_z(idx_z), corr_z(idx_z), 'r*', lag_z(idx_z_raw), corr_z(idx_z_raw), 'g*') %z +% hold ("on") +% m_label_ms = strcat(" mean ms: ", num2str(Zwindow_mean_ts_diff)); +% m_label_bpm = strcat(" mean bpm: ", num2str(Zwindow_mean_bpm)); +% title(strcat(" ", m_label_ms, " ", m_label_bpm)); +% hold ("off"); %magnitude Mwindow = magnitude(i-window_size:i); Mwindow_mean_ts_diff = mean(diff(lag_mag(idx_mag)* sample_rate_ms)); Mwindow_mean_bpm = (60000 / (Mwindow_mean_ts_diff)); - figure(14); - plot(lag_mag, corr_mag, lag_mag(idx_mag), corr_mag(idx_mag), 'r*', lag_mag(idx_mag_raw), corr_mag(idx_mag_raw), 'g*') %z - hold ("on") - m_label_ms = strcat(" mean ms: ", num2str(Mwindow_mean_ts_diff)); - m_label_bpm = strcat(" mean bpm: ", num2str(Mwindow_mean_bpm)); - title(strcat(" ", m_label_ms, " ", m_label_bpm)); - hold ("off"); - - +% figure(14); +% plot(lag_mag, corr_mag, lag_mag(idx_mag), corr_mag(idx_mag), 'r*', lag_mag(idx_mag_raw), corr_mag(idx_mag_raw), 'g*') %z +% hold ("on") +% m_label_ms = strcat(" mean ms: ", num2str(Mwindow_mean_ts_diff)); +% m_label_bpm = strcat(" mean bpm: ", num2str(Mwindow_mean_bpm)); +% title(strcat(" ", m_label_ms, " ", m_label_bpm)); +% hold ("off"); %breakpoints dummy for testing if(length(idx_x) > length(idx_x_raw)) @@ -277,11 +301,22 @@ for file = files' if(isnan(window_mean_ts_diff) || isnan(window_mean_bpm)) %do nothing - else + else + gtError_1D = [gtError_1D, abs(window_mean_bpm - curGtBpm)]; bpm_per_window_ms = [bpm_per_window_ms, window_mean_ts_diff]; bpm_per_window = [bpm_per_window, window_mean_bpm]; end + + %3D mean + if(isnan(Dwindow_mean_bpm)) + %nothing + else + gtError_3D = [gtError_3D, abs(Dwindow_mean_bpm - curGtBpm)]; + bpm_3D = [bpm_3D, Dwindow_mean_bpm]; + ms_3D = [ms_3D, Dwindow_mean_ts_diff]; + end + %TODO: if correlation value is lower then a treshhold, we are not conducting TODO: change to a real classification instead of a treshhold. end @@ -290,22 +325,70 @@ for file = files' %TODO: smooth the results using a moving avg or 1d kalman filter.(transition for kalman could be adding the last measured value) %remove the first 40% of the results, due to starting delays while recording. - number_to_remove = round(abs(0.1 * length(bpm_per_window_ms))); - num_all = length(bpm_per_window_ms); - bpm_per_window_ms = bpm_per_window_ms(number_to_remove:num_all); - bpm_per_window = bpm_per_window(number_to_remove:num_all); + %number_to_remove = round(abs(0.1 * length(bpm_per_window_ms))); + %num_all = length(bpm_per_window_ms); + %bpm_per_window_ms = bpm_per_window_ms(number_to_remove:num_all); + %bpm_per_window = bpm_per_window(number_to_remove:num_all); mean_final_ms = mean(bpm_per_window_ms); std_final_ms = std(bpm_per_window_ms); mean_final_bpm = mean(bpm_per_window); std_final_bpm = std(bpm_per_window); + + mean_final_error_1D = mean(gtError_1D); + std_final_error_1D = std(gtError_1D); + + mean_final_ms_3D = mean(ms_3D); + std_final_ms_3D = std(ms_3D); + + mean_final_bpm_3D = mean(bpm_3D); + std_final_bpm_3D = std(bpm_3D); + + mean_final_error_3D = mean(gtError_3D); + std_final_error_3D = std(gtError_3D); - fprintf('%s: mean = %f bpm (%f ms) stddev = %f bpm (%f ms)\n', strrep(regexprep(filename,'^.*recording_',''),'.txt',''), mean_final_bpm, mean_final_ms, std_final_bpm, std_final_ms); - + fprintf('%s: mean = %f bpm (%f bpm) stddev = %f bpm (%f bpm) --- 1D\n', strrep(regexprep(filename,'^.*recording_',''),'.txt',''), mean_final_error_1D, mean_final_bpm, std_final_error_1D, std_final_bpm); + fprintf('%s: mean = %f bpm (%f bpm) stddev = %f bpm (%f bpm) --- 3D\n', strrep(regexprep(filename,'^.*recording_',''),'.txt',''), mean_final_error_3D, mean_final_bpm_3D, std_final_error_3D, std_final_bpm_3D); + end +% %1D fft - nicht so der brüller +% z_fft = fft(m(i-window_size:i,5)); +% L = length(z_fft); +% Fs = 250; +% P2 = abs(z_fft/L); +% P1 = P2(1:L/2+1); +% P1(2:end-1) = 2*P1(2:end-1); +% f = Fs*(0:(L/2))/L; %nyquist frequence +% +% figure(66); +% plot(f, P1); +% +% %3D fft +% m_3D = m(i-window_size:i, 3); +% m_3D(:,:,2) = m(i-window_size:i, 4); +% m_3D(:,:,3) = m(i-window_size:i, 5); +% +% fft_3D = fftn(m_3D); +% +% %2D fft +% fft_xy = fft2(m(i-window_size:i, 3:4)); +% fft_yz = fft2(m(i-window_size:i, 3:4)); +% +% fft_test = fft(m(i-window_size:i, 3:5),[],2); +% +% figure(60); +% imagesc(abs(fftshift(fft_test))); +% +% figure(61); +% imagesc(abs(fftshift(fft_xy))); + + + + + diff --git a/matlab/distCorr.m b/matlab/distCorr.m new file mode 100644 index 0000000..a072192 --- /dev/null +++ b/matlab/distCorr.m @@ -0,0 +1,41 @@ +% Autocorrelation for Points based on the distance between those points +% data is the sensor data. one point per row. +% lag_size is the number of required lags +function [corr, lag] = distCorr(data, lag_size) + + n = length(data); + + %if lag size is bigger as data, then use data length - 1 + % -1 because the max. lag is -1 of the data length.. + if(lag_size >= n) + lag_size = n - 1; % + end + + %init + corr = zeros(lag_size + 1, 1); % +1, because the first index is lag 0. + lag = (-lag_size:1:lag_size)'; + + %mean shift and l2 normalization + %data = data - mean(data); + %data = data / norm(data); + + for j = 1:lag_size + 1 % +1, because the first index is lag 0. + dist = zeros(n - abs(j), 1); + idx = 1; + + for i = j:n + %x_i * x_i-(j-1) + dist(idx) = norm(data(i, :) - data(i-(j-1), :)); + idx = idx + 1; + end + + corr(j) = geomean(dist); + end + + %mirror corr(2:512) and put it infront + corr = [flipud(corr(2:end)); corr]; + + %to [0, 1] + corr = ((corr .* -1) / max(corr)) + 1; + +end \ No newline at end of file