toni/Dirigent
1
0
Fork 0
Code Issues 17 Pull Requests Releases 2 Wiki Activity

added distance based correlation function in matlab and java

Browse Source
This commit is contained in:
toni
2018-12-14 16:29:34 +01:00
parent 42e38bd929
commit 6bb8bb6b4f
11 changed files with 1681 additions and 289 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
android/ConductorsWatch/app/src/main/java/de/tonifetzer/conductorswatch/MainActivity.java
View File

@@ -247,6 +247,7 @@ public class MainActivity extends WearableActivity implements WorkerFragment.OnF
mLastSendBPM = currentBPM; mLastSendBPM = currentBPM;
} }
//künstlicher delay.
if(mReadyToSend){ if(mReadyToSend){
mReadyToSend = false; mReadyToSend = false;
new java.util.Timer().schedule( new java.util.Timer().schedule(

2
android/ConductorsWatch/build.gradle
View File

@@ -10,7 +10,7 @@ buildscript {
google() google()
} }
dependencies { 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 // NOTE: Do not place your application dependencies here; they belong

158
java/src/main/java/Main.java
View File

@@ -21,28 +21,19 @@ public class Main {
public static void main(String [ ] args) { 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/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/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(); File[] listOfFiles = folder.listFiles();
Arrays.sort(listOfFiles); Arrays.sort(listOfFiles);
//calc results
// //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();
// iterate trough files in measurements folder // iterate trough files in measurements folder
for (File file : listOfFiles) { for (File file : listOfFiles) {
if (file.isFile() && file.getName().contains(".csv")) { 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); AccelerometerWindowBuffer accWindowBuffer = new AccelerometerWindowBuffer(6000, 750);
BpmEstimator bpmEstimator = new BpmEstimator(accWindowBuffer, 0, 5000); BpmEstimator bpmEstimator = new BpmEstimator(accWindowBuffer, 0, 5000);
Butterworth butterLowpass = new Butterworth();
//read the file line by line //read the file line by line
try (BufferedReader br = new BufferedReader(new FileReader(file))) { 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! //read the first three lines and print out what file it is!
String comment = br.readLine(); String comment = br.readLine();
br.readLine(); br.readLine();
br.readLine(); String groundTruth = br.readLine();
System.out.println(comment); System.out.println(comment);
//long startTs = Long.parseLong(br.readLine().split(";")[0]);
for (String line; (line = br.readLine()) != null; ) { for (String line; (line = br.readLine()) != null; ) {
// process the line. // process the line.
String[] measurement = line.split(";"); String[] measurement = line.split(";");
//if linear acc //if linear acc
long ts = 0;
if(measurement[1].equals("3")){ if(measurement[1].equals("3")){
long ts = Long.parseLong(measurement[0]); ts = Long.parseLong(measurement[0]);
double x = Double.parseDouble(measurement[2]); double x = Double.parseDouble(measurement[2]);
double y = Double.parseDouble(measurement[3]); double y = Double.parseDouble(measurement[3]);
double z = Double.parseDouble(measurement[4]); double z = Double.parseDouble(measurement[4]);
@@ -88,120 +81,18 @@ public class Main {
bpmList.add(bpm160); bpmList.add(bpm160);
bpmList.add(bpm200); bpmList.add(bpm200);
// //TODO: plot the different autocorrelation of different window sizes while(bpmList.remove(Double.valueOf(-1))) {}
// // Draw the data for the complete 6 second window Utils.removeOutliersZScore(bpmList, 3.4);
// 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<Integer> 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<Integer> 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<Integer> 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());
//fill hols improve peaks double bpmMean = Utils.mean(bpmList);
double magMean = bpmEstimator.getMagnitudeMean();
double bpmDist = bpmEstimator.getDistEstimation();
//estimate bpm between detected peaks //double bpmSingle = bpmEstimator.getBestSingleAxis();
//System.out.println("BPM-X: " + pX.getBPM(bpmEstimator.getSampleRate_ms())); double bpmAllAverage = bpmEstimator.getAverageOfAllWindows();
//System.out.println("BPM-Y: " + pY.getBPM(bpmEstimator.getSampleRate_ms()));
//System.out.println("BPM-Z: " + pZ.getBPM(bpmEstimator.getSampleRate_ms())); System.out.println( ts + " all: " + Math.round(bpmMean) + " avg_all: " + Math.round(bpmAllAverage) + " 3D: " + Math.round(bpmDist));
System.out.println(" ");
int dummyForBreakpoint = 0; int dummyForBreakpoint = 0;
} }
@@ -213,17 +104,24 @@ public class Main {
//System.out.println("MEAN BPM: " + Math.round(meanBPM)); //System.out.println("MEAN BPM: " + Math.round(meanBPM));
//System.out.println("MEDIAN BPM: " + Math.round(medianBPM)); //System.out.println("MEDIAN BPM: " + Math.round(medianBPM));
if(Utils.DEBUG_MODE){
bpmEstimator.closeDebugWindows();
}
} catch (IOException e) { } catch (IOException e) {
e.printStackTrace(); e.printStackTrace();
} }
} }
try {
System.in.read();
} catch (IOException e) { // try {
e.printStackTrace(); // System.in.read();
} // } catch (IOException e) {
// e.printStackTrace();
// }
} }
} }
} }

9
java/src/main/java/bpmEstimation/AccelerometerInterpolator.java
View File

@@ -11,11 +11,12 @@ public class AccelerometerInterpolator {
private double[] mY; private double[] mY;
private double[] mZ; private double[] mZ;
private long[] mTsInterp; private long[] mTsInterp;
private int size;
public AccelerometerInterpolator(AccelerometerWindowBuffer ab, double sampleRate_ms){ public AccelerometerInterpolator(AccelerometerWindowBuffer ab, double sampleRate_ms){
long size = (ab.getYongest().ts - (ab.getOldest().ts - (long) sampleRate_ms)) / (long) sampleRate_ms; this.size = (int) ((ab.getYongest().ts - (ab.getOldest().ts - (long) sampleRate_ms)) / (long) sampleRate_ms);
mTsInterp = new long[(int)size]; mTsInterp = new long[size];
int j = 0; int j = 0;
for(long i = ab.getOldest().ts; i <= ab.getYongest().ts; i += sampleRate_ms){ for(long i = ab.getOldest().ts; i <= ab.getYongest().ts; i += sampleRate_ms){
mTsInterp[j++] = i; mTsInterp[j++] = i;
@@ -42,6 +43,10 @@ public class AccelerometerInterpolator {
return mZ; return mZ;
} }
public int size(){
return size;
}
private static double[] interpLinear(double[] x, double[] y, double[] xi) throws IllegalArgumentException { private static double[] interpLinear(double[] x, double[] y, double[] xi) throws IllegalArgumentException {
if (x.length != y.length) { if (x.length != y.length) {
throw new IllegalArgumentException("X and Y must be the same length"); throw new IllegalArgumentException("X and Y must be the same length");

192
java/src/main/java/bpmEstimation/BpmEstimator.java
View File

@@ -3,7 +3,6 @@ package bpmEstimation;
import static utilities.Utils.DEBUG_MODE; import static utilities.Utils.DEBUG_MODE;
import utilities.MovingFilter; import utilities.MovingFilter;
import utilities.SimpleKalman;
import utilities.Utils; import utilities.Utils;
import uk.me.berndporr.iirj.*; import uk.me.berndporr.iirj.*;
@@ -18,10 +17,11 @@ public class BpmEstimator {
private AccelerometerWindowBuffer mBuffer; private AccelerometerWindowBuffer mBuffer;
private static double mSampleRate_ms; private static double mSampleRate_ms;
private LinkedList<Double> mBpmHistory_X; private BpmHistory mBpmHistory_X;
private LinkedList<Double> mBpmHistory_Y; private BpmHistory mBpmHistory_Y;
private LinkedList<Double> mBpmHistory_Z; private BpmHistory mBpmHistory_Z;
private LinkedList<Double> mBpmHistory_Mag; private BpmHistory mBpmHistory_Mag;
private BpmHistory mBpmHistory_Dist;
private LinkedList<Double> mBpmHistory; private LinkedList<Double> mBpmHistory;
private int mResetCounter; private int mResetCounter;
@@ -35,6 +35,10 @@ public class BpmEstimator {
private Butterworth mButter_Z; private Butterworth mButter_Z;
private Butterworth mButter_Mag; private Butterworth mButter_Mag;
//todo: hack
private double magnitudeMean = 1.0;
//Debugging stuff //Debugging stuff
private Utils.DebugPlotter plotter; private Utils.DebugPlotter plotter;
@@ -42,12 +46,13 @@ public class BpmEstimator {
mBuffer = windowBuffer; mBuffer = windowBuffer;
mSampleRate_ms = sampleRate_ms; mSampleRate_ms = sampleRate_ms;
mBpmHistory_X = new LinkedList<>(); mBpmHistory_X = new BpmHistory();
mBpmHistory_Y = new LinkedList<>(); mBpmHistory_Y = new BpmHistory();
mBpmHistory_Z = new LinkedList<>(); mBpmHistory_Z = new BpmHistory();
mBpmHistory_Mag = new LinkedList<>(); mBpmHistory_Mag = new BpmHistory();
mBpmHistory_Dist = new BpmHistory();
mBpmHistory = new LinkedList<>(); mBpmHistory = new BpmHistory();
mResetCounter = 0; mResetCounter = 0;
mResetLimit_ms = resetAfter_ms; mResetLimit_ms = resetAfter_ms;
@@ -70,6 +75,8 @@ public class BpmEstimator {
public double estimate(int length_ms, int overlap_ms){ 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); AccelerometerWindowBuffer tmpBuffer = AccelerometerWindowBuffer.getNewInstance(mBuffer, length_ms, overlap_ms);
double sampleRate = mSampleRate_ms; double sampleRate = mSampleRate_ms;
@@ -89,25 +96,29 @@ public class BpmEstimator {
mButter_X.lowPass(1,(1/(sampleRate / 1000))/2, 1); mButter_X.lowPass(1,(1/(sampleRate / 1000))/2, 1);
mButter_Y.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_Z.lowPass(1,(1/(sampleRate / 1000))/2, 1);
mButter_Mag.lowPass(1,(1/(sampleRate / 1000))/2, 1);
int n = interp.getX().length; int n = interp.getX().length;
double[] xButter = new double[n]; double[] xButter = new double[n];
double[] yButter = new double[n];; double[] yButter = new double[n];
double[] zButter = new double[n];; double[] zButter = new double[n];
double[] magButter = new double[n];
for(int i = 0; i < interp.getX().length; ++i){ for(int i = 0; i < interp.getX().length; ++i){
xButter[i] = mButter_X.filter(interp.getX()[i]); xButter[i] = mButter_X.filter(interp.getX()[i]);
yButter[i] = mButter_Y.filter(interp.getY()[i]); yButter[i] = mButter_Y.filter(interp.getY()[i]);
zButter[i] = mButter_Z.filter(interp.getZ()[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 //auto correlation
double[] xAutoCorr = new AutoCorrelation(xButter, tmpBuffer.size()).getCorr(); double[] xAutoCorr = new AutoCorrelation(xButter, tmpBuffer.size()).getCorr();
double[] yAutoCorr = new AutoCorrelation(yButter, tmpBuffer.size()).getCorr(); double[] yAutoCorr = new AutoCorrelation(yButter, tmpBuffer.size()).getCorr();
double[] zAutoCorr = new AutoCorrelation(zButter, 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 //find a peak within range of 250 ms
int peakWidth = (int) Math.round(250 / sampleRate); 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 zPeaks = new Peaks(zAutoCorr, peakWidth, 0.1f, 0, false);
Peaks magPeaks = new Peaks(magAutoCorr, 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_X.add(xPeaks.getBPM(sampleRate));
mBpmHistory_Y.add(yPeaks.getBPM(sampleRate)); mBpmHistory_Y.add(yPeaks.getBPM(sampleRate));
mBpmHistory_Z.add(zPeaks.getBPM(sampleRate)); mBpmHistory_Z.add(zPeaks.getBPM(sampleRate));
mBpmHistory_Mag.add(magPeaks.getBPM(sampleRate)); mBpmHistory_Mag.add(magPeaks.getBPM(sampleRate));
mBpmHistory_Dist.add(distPeaks.getBPM(sampleRate));
if(DEBUG_MODE){ if(DEBUG_MODE){
plotter.setPlotRawX(interp.getTs(), interp.getX()); //plotter.setPlotRawX(interp.getTs(), interp.getX());
plotter.setPlotRawY(interp.getTs(), interp.getY()); //plotter.setPlotRawY(interp.getTs(), interp.getY());
plotter.setPlotRawZ(interp.getTs(), interp.getZ()); //plotter.setPlotRawZ(interp.getTs(), interp.getZ());
plotter.setPlotRawMag(interp.getTs(), magRaw); //plotter.setPlotRawMag(interp.getTs(), magRaw);
plotter.setPlotButterX(interp.getTs(), xButter); //plotter.setPlotButterX(interp.getTs(), xButter);
plotter.setPlotButterY(interp.getTs(), yButter); //plotter.setPlotButterY(interp.getTs(), yButter);
plotter.setPlotButterZ(interp.getTs(), zButter); //plotter.setPlotButterZ(interp.getTs(), zButter);
plotter.setPlotButterMag(interp.getTs(), magButter); //plotter.setPlotButterMag(interp.getTs(), magButter);
plotter.setPlotCorrX(xAutoCorr, xPeaks); plotter.setPlotCorrX(xAutoCorr, xPeaks);
plotter.setPlotCorrY(yAutoCorr, yPeaks); plotter.setPlotCorrY(yAutoCorr, yPeaks);
plotter.setPlotCorrZ(zAutoCorr, zPeaks); plotter.setPlotCorrZ(zAutoCorr, zPeaks);
plotter.setPlotCorrMag(magAutoCorr, magPeaks); 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){ if(estimatedBPM != -1){
//moving avg (lohnt dann, wenn wir viele daten haben) //moving avg (lohnt dann, wenn wir viele daten haben)
@@ -178,6 +197,22 @@ public class BpmEstimator {
return mBpmHistory.getLast(); 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(){ public double getMeanBpm(){
return Utils.mean(mBpmHistory); return Utils.mean(mBpmHistory);
} }
@@ -186,9 +221,84 @@ public class BpmEstimator {
return Utils.median(mBpmHistory); 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: 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 //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; int cntNumAxis = 0;
double sumCorr = 0; //to prevent division by zero double sumCorr = 0; //to prevent division by zero
@@ -234,6 +344,19 @@ public class BpmEstimator {
sumNumInter += corrNumInterZ; 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 //no peaks, reject
if(cntNumAxis == 0){ if(cntNumAxis == 0){
//throw new IllegalArgumentException("All Peaks are empty! -> Reject Estimation"); //throw new IllegalArgumentException("All Peaks are empty! -> Reject Estimation");
@@ -256,23 +379,28 @@ public class BpmEstimator {
//values to low, reject //values to low, reject
//TODO: this is a pretty simple assumption. first shot! //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; return -1;
} }
double quantityX = ((corrMeanX / sumCorr) + (corrRmsX / sumRms) + (corrNumInterX / sumNumInter)) / cntNumAxis; double quantityX = ((corrMeanX / sumCorr) + (corrRmsX / sumRms) + (corrNumInterX / sumNumInter)) / cntNumAxis;
double quantityY = ((corrMeanY / sumCorr) + (corrRmsY / sumRms) + (corrNumInterY / sumNumInter)) / cntNumAxis; double quantityY = ((corrMeanY / sumCorr) + (corrRmsY / sumRms) + (corrNumInterY / sumNumInter)) / cntNumAxis;
double quantityZ = ((corrMeanZ / sumCorr) + (corrRmsZ / sumRms) + (corrNumInterZ / 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 //get best axis by quantity and estimate bpm
if(quantityX > quantityY && quantityX > quantityZ){ if(quantityX > quantityY && quantityX > quantityZ){
return mBpmHistory_X.getLast(); return mBpmHistory_X.getLast();
} }
else if(quantityY > quantityZ){ else if(quantityY > quantityZ ){
return mBpmHistory_Y.getLast(); return mBpmHistory_Y.getLast();
} }
else { else {
return mBpmHistory_Z.getLast(); return mBpmHistory_Z.getLast();
} }
} }
public void closeDebugWindows(){
plotter.close();
}
} }

53
java/src/main/java/bpmEstimation/BpmHistory.java Normal file
View File

@@ -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<Double> {
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);
}
}

66
java/src/main/java/bpmEstimation/DistanceCorrelation.java Normal file
View File

@@ -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;
}
}

1029
java/src/main/java/utilities/Plot.java Normal file
View File

File diff suppressed because it is too large Load Diff

108
java/src/main/java/utilities/Utils.java
View File

@@ -5,12 +5,10 @@ import bpmEstimation.AccelerometerWindowBuffer;
import bpmEstimation.Peaks; import bpmEstimation.Peaks;
import javax.swing.*; import javax.swing.*;
import java.awt.*; import java.awt.*;
import java.awt.event.WindowEvent;
import java.awt.image.BufferedImage; import java.awt.image.BufferedImage;
import java.text.Collator; import java.util.*;
import java.util.Arrays; import java.util.List;
import java.util.Collections;
import java.util.Comparator;
import java.util.LinkedList;
import java.util.stream.IntStream; import java.util.stream.IntStream;
//TODO: change from double to generic type //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)); 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<Double> 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<Double> data){
return Math.sqrt(var(data));
}
public static double sqr(double x) { public static double sqr(double x) {
return x * x; return x * x;
} }
@@ -94,7 +124,7 @@ public class Utils {
return sum(data) / data.size(); return sum(data) / data.size();
} }
public static double median(LinkedList<Double> data){ public static double median(List<Double> data){
data.sort(Comparator.naturalOrder()); data.sort(Comparator.naturalOrder());
double median; double median;
@@ -106,6 +136,18 @@ public class Utils {
return median; return median;
} }
public static double mad(List<Double> data){
double median = median(data);
java.util.List<Double> tmpList = new ArrayList<>();
for(double value : data){
tmpList.add(Math.abs(value - median));
}
return median(tmpList);
}
//TODO: Could be slow.. faster method? //TODO: Could be slow.. faster method?
public static double geometricMean(double[] data) { public static double geometricMean(double[] data) {
double sum = data[0]; double sum = data[0];
@@ -115,6 +157,17 @@ public class Utils {
return Math.pow(sum, 1.0 / data.length); 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){ public static int intersectionNumber(double[] signal, double border){
int cnt = 0; int cnt = 0;
boolean isSmallerValue = false; boolean isSmallerValue = false;
@@ -170,6 +223,20 @@ public class Utils {
return newArray; return newArray;
} }
public static void removeOutliersZScore(List<Double> data, double score) {
if(!data.isEmpty()){
double median = median(data);
double mad = mad(data);
for(Iterator<Double> 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){ public static double magnitude(double x, double y, double z){
return Math.sqrt((x*x) + (y*y) + (z*z)); return Math.sqrt((x*x) + (y*y) + (z*z));
} }
@@ -287,6 +354,7 @@ public class Utils {
private Utils.ShowPNG plotCorrY; private Utils.ShowPNG plotCorrY;
private Utils.ShowPNG plotCorrZ; private Utils.ShowPNG plotCorrZ;
private Utils.ShowPNG plotCorrMag; private Utils.ShowPNG plotCorrMag;
private Utils.ShowPNG plotCorr3D;
public DebugPlotter(){ public DebugPlotter(){
@@ -302,6 +370,7 @@ public class Utils {
plotCorrY = new Utils.ShowPNG(); plotCorrY = new Utils.ShowPNG();
plotCorrZ = new Utils.ShowPNG(); plotCorrZ = new Utils.ShowPNG();
plotCorrMag = new Utils.ShowPNG(); plotCorrMag = new Utils.ShowPNG();
plotCorr3D = new Utils.ShowPNG();
} }
public void setPlotRawX(long[] ts, double[] data){ public void setPlotRawX(long[] ts, double[] data){
@@ -316,8 +385,7 @@ public class Utils {
plotRawZ.plotData("Raw Data Z", "z", ts, data); plotRawZ.plotData("Raw Data Z", "z", ts, data);
} }
public void setPlotRawMag(long[] ts, double[] data){ public void setPlotRawMag(long[] ts, double[] data){plotRawMag.plotData("Raw Data Magnitude", "mag", ts, data);
plotRawMag.plotData("Raw Data Magnitude", "mag", ts, data);
} }
public void setPlotButterX(long[] ts, double[] data){ public void setPlotButterX(long[] ts, double[] data){
@@ -325,11 +393,11 @@ public class Utils {
} }
public void setPlotButterY(long[] ts, double[] data){ 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){ 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){ public void setPlotButterMag(long[] ts, double[] data){
@@ -352,6 +420,26 @@ public class Utils {
plotCorrMag.plotCorrWithPeaks("Autocorr Mag", "mag", data,2, peaks, true); 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));
}
} }

299
matlab/AutoCorrMethodNew_Watch.m
View File

@@ -1,58 +1,47 @@
%We are using a threshold-based version for bpm estimation %using autocorrelation to estimate the current bmp within some fixed window
%Only the z axis of the acc is used
%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 sensor files
%files = dir(fullfile('../../measurements/2017.06/mSensor/', '*.csv'));
%load file provided by the sensor readout app %files = natsortfiles(dir(fullfile('../../measurements/2017.06/lgWear/', '*.csv')));
% 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'));
%files = dir(fullfile('../../measurements/wearR/', '*.csv')); %files = dir(fullfile('../../measurements/wearR/', '*.csv'));
%files = dir(fullfile('../../measurements/peter_failed/', '*.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' for file = files'
filename = [file.folder '/' file.name]; 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 %draw the raw acc data
m_idx = []; 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, :); m = measurements(m_idx, :);
%Interpolate to generate a constant sample rate to 250hz (4ms per sample) %Interpolate to generate a constant sample rate to 250hz (4ms per sample)
@@ -66,52 +55,73 @@ for file = files'
%put all together again %put all together again
m = [t_interp', t_interp', m_interp]; m = [t_interp', t_interp', m_interp];
figure(1); % figure(1);
plot(m(:,1),m(:,3)) %x % plot(m(:,1),m(:,3)) %x
legend("x", "location", "eastoutside"); % 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");
figure(2); %waitforbuttonpress();
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");
waitforbuttonpress();
%save timestamps %save timestamps
timestamps = m(:,1); timestamps = m(:,1);
data = m(:,3); %only z data = m(:,3); %only z
%TODO: Different window sizes for periods under 16.3 s %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; overlap = 256;
bpm_per_window_ms = []; bpm_per_window_ms = [];
bpm_per_window = []; 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 %wait until window is filled with new data
if(mod(i,overlap) == 0) 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 %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_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/4), "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/4), "coeff"); [corr_z, lag_z] = xcov(m(i-window_size:i,5), (window_size/2), "coeff");
[corr_mag, lag_mag] = xcov(magnitude(i-window_size:i), (window_size/4), "coeff");
%magnitude
[corr_mag, lag_mag] = xcov(magnitude(i-window_size:i), (window_size/2), "coeff");
%autocorrelation of the autocorrelation?! %TODO: stichwort spatial autocorrelation
%[corr_corr_x, lag_lag_x] = xcov(corr_x, length(corr_x), "coeff"); %figure(77);
%[corr_corr_y, lag_lag_y] = xcov(corr_y, length(corr_x), "coeff"); %scatter3(timestamps(i-window_size:i), m(i-window_size:i,4), m(i-window_size:i,5));
%[corr_corr_z, lag_lag_z] = xcov(corr_z, length(corr_x), "coeff");
%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_x_pos = corr_x;
corr_y_pos = corr_y; corr_y_pos = corr_y;
@@ -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_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_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_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_x_raw = sort(idx_x_raw);
idx_y_raw = sort(idx_y_raw); idx_y_raw = sort(idx_y_raw);
idx_z_raw = sort(idx_z_raw); idx_z_raw = sort(idx_z_raw);
idx_mag_raw = sort(idx_mag_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_x = findFalseDetectedPeaks(idx_x_raw, lag_x, corr_x);
idx_y = findFalseDetectedPeaks(idx_y_raw, lag_y, corr_y); idx_y = findFalseDetectedPeaks(idx_y_raw, lag_y, corr_y);
idx_z = findFalseDetectedPeaks(idx_z_raw, lag_z, corr_z); idx_z = findFalseDetectedPeaks(idx_z_raw, lag_z, corr_z);
idx_mag = findFalseDetectedPeaks(idx_mag_raw, lag_mag, corr_mag); 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 = 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_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)); Xwindow_mean_bpm = (60000 / (Xwindow_mean_ts_diff));
figure(11); % 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 % 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") % hold ("on")
m_label_ms = strcat(" mean ms: ", num2str(Xwindow_mean_ts_diff)); % m_label_ms = strcat(" mean ms: ", num2str(Xwindow_mean_ts_diff));
m_label_bpm = strcat(" mean bpm: ", num2str(Xwindow_mean_bpm)); % m_label_bpm = strcat(" mean bpm: ", num2str(Xwindow_mean_bpm));
title(strcat(" ", m_label_ms, " ", m_label_bpm)); % title(strcat(" ", m_label_ms, " ", m_label_bpm));
hold ("off"); % hold ("off");
Ywindow = m(i-window_size:i,4); Ywindow = m(i-window_size:i,4);
Ywindow_mean_ts_diff = mean(diff(lag_y(idx_y) * sample_rate_ms)); Ywindow_mean_ts_diff = mean(diff(lag_y(idx_y) * sample_rate_ms));
Ywindow_mean_bpm = (60000 / (Ywindow_mean_ts_diff)); Ywindow_mean_bpm = (60000 / (Ywindow_mean_ts_diff));
figure(12); % 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 % 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") % hold ("on")
m_label_ms = strcat(" mean ms: ", num2str(Ywindow_mean_ts_diff)); % m_label_ms = strcat(" mean ms: ", num2str(Ywindow_mean_ts_diff));
m_label_bpm = strcat(" mean bpm: ", num2str(Ywindow_mean_bpm)); % m_label_bpm = strcat(" mean bpm: ", num2str(Ywindow_mean_bpm));
title(strcat(" ", m_label_ms, " ", m_label_bpm)); % title(strcat(" ", m_label_ms, " ", m_label_bpm));
hold ("off"); % hold ("off");
Zwindow = m(i-window_size:i,5); Zwindow = m(i-window_size:i,5);
Zwindow_mean_ts_diff = mean(diff(lag_z(idx_z)* sample_rate_ms)); Zwindow_mean_ts_diff = mean(diff(lag_z(idx_z)* sample_rate_ms));
Zwindow_mean_bpm = (60000 / (Zwindow_mean_ts_diff)); Zwindow_mean_bpm = (60000 / (Zwindow_mean_ts_diff));
figure(13); % 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 % 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") % hold ("on")
m_label_ms = strcat(" mean ms: ", num2str(Zwindow_mean_ts_diff)); % m_label_ms = strcat(" mean ms: ", num2str(Zwindow_mean_ts_diff));
m_label_bpm = strcat(" mean bpm: ", num2str(Zwindow_mean_bpm)); % m_label_bpm = strcat(" mean bpm: ", num2str(Zwindow_mean_bpm));
title(strcat(" ", m_label_ms, " ", m_label_bpm)); % title(strcat(" ", m_label_ms, " ", m_label_bpm));
hold ("off"); % hold ("off");
%magnitude %magnitude
Mwindow = magnitude(i-window_size:i); Mwindow = magnitude(i-window_size:i);
Mwindow_mean_ts_diff = mean(diff(lag_mag(idx_mag)* sample_rate_ms)); Mwindow_mean_ts_diff = mean(diff(lag_mag(idx_mag)* sample_rate_ms));
Mwindow_mean_bpm = (60000 / (Mwindow_mean_ts_diff)); Mwindow_mean_bpm = (60000 / (Mwindow_mean_ts_diff));
figure(14); % 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 % 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") % hold ("on")
m_label_ms = strcat(" mean ms: ", num2str(Mwindow_mean_ts_diff)); % m_label_ms = strcat(" mean ms: ", num2str(Mwindow_mean_ts_diff));
m_label_bpm = strcat(" mean bpm: ", num2str(Mwindow_mean_bpm)); % m_label_bpm = strcat(" mean bpm: ", num2str(Mwindow_mean_bpm));
title(strcat(" ", m_label_ms, " ", m_label_bpm)); % title(strcat(" ", m_label_ms, " ", m_label_bpm));
hold ("off"); % hold ("off");
%breakpoints dummy for testing %breakpoints dummy for testing
if(length(idx_x) > length(idx_x_raw)) if(length(idx_x) > length(idx_x_raw))
@@ -278,10 +302,21 @@ for file = files'
if(isnan(window_mean_ts_diff) || isnan(window_mean_bpm)) if(isnan(window_mean_ts_diff) || isnan(window_mean_bpm))
%do nothing %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_ms = [bpm_per_window_ms, window_mean_ts_diff];
bpm_per_window = [bpm_per_window, window_mean_bpm]; bpm_per_window = [bpm_per_window, window_mean_bpm];
end 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. %TODO: if correlation value is lower then a treshhold, we are not conducting TODO: change to a real classification instead of a treshhold.
end end
@@ -290,10 +325,10 @@ 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) %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. %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))); %number_to_remove = round(abs(0.1 * length(bpm_per_window_ms)));
num_all = 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_ms = bpm_per_window_ms(number_to_remove:num_all);
bpm_per_window = bpm_per_window(number_to_remove:num_all); %bpm_per_window = bpm_per_window(number_to_remove:num_all);
mean_final_ms = mean(bpm_per_window_ms); mean_final_ms = mean(bpm_per_window_ms);
std_final_ms = std(bpm_per_window_ms); std_final_ms = std(bpm_per_window_ms);
@@ -301,11 +336,59 @@ for file = files'
mean_final_bpm = mean(bpm_per_window); mean_final_bpm = mean(bpm_per_window);
std_final_bpm = std(bpm_per_window); std_final_bpm = std(bpm_per_window);
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); 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 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 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)));

41
matlab/distCorr.m Normal file
View File

@@ -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