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

ref #6 -refactoring

Browse Source 
- added interp1
- added peak finder
- added bpm estimation
- added findbestaxis
This commit is contained in:
toni
2017-12-17 00:35:57 +01:00
parent efa6f95972
commit 709a846a91
9 changed files with 793 additions and 305 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -1,3 +1,6 @@
import bpmEstimation.*;
import utilities.Utils;
import java.awt.*;
import java.io.BufferedReader;
import java.io.File;
@@ -24,8 +27,7 @@ public class Main {
for (File file : listOfFiles) {
if (file.isFile() && file.getName().contains(".csv")) {
Utils.AccelerometerWindowBuffer accWindowBuffer = new Utils.AccelerometerWindowBuffer(4096, 256);
AccelerometerWindowBuffer accWindowBuffer = new AccelerometerWindowBuffer(4096, 256);
//read the file line by line
try (BufferedReader br = new BufferedReader(new FileReader(file))) {
@@ -36,86 +38,83 @@ public class Main {
//if linear acc
if(measurement[1].equals("2")){
long ts = Long.parseLong(measurement[0]);
float x = Float.parseFloat(measurement[2]);
float y = Float.parseFloat(measurement[3]);
float z = Float.parseFloat(measurement[4]);
accWindowBuffer.add(new Utils.AccelerometerData(ts, x, y, z));
double x = Double.parseDouble(measurement[2]);
double y = Double.parseDouble(measurement[3]);
double z = Double.parseDouble(measurement[4]);
accWindowBuffer.add(new AccelerometerData(ts, x, y, z));
}
//do calculation stuff
if(accWindowBuffer.isNextWindowReady()){
AccelerometerInterpolator acInterp = new AccelerometerInterpolator(accWindowBuffer, 6);
//print raw x,y,z
double[] dTs = IntStream.range(0, accWindowBuffer.getTs().length).mapToDouble(i -> accWindowBuffer.getTs()[i]).toArray();
double[] dX = IntStream.range(0, accWindowBuffer.getX().length).mapToDouble(i -> accWindowBuffer.getX()[i]).toArray();
double[] dY = IntStream.range(0, accWindowBuffer.getY().length).mapToDouble(i -> accWindowBuffer.getY()[i]).toArray();
double[] dZ = IntStream.range(0, accWindowBuffer.getZ().length).mapToDouble(i -> accWindowBuffer.getZ()[i]).toArray();
double[] dTsInterp = IntStream.range(0, acInterp.getTs().length).mapToDouble(i -> acInterp.getTs()[i]).toArray();
Plot plotRaw = Plot.plot(Plot.plotOpts().
title("Raw Acc Data").
legend(Plot.LegendFormat.BOTTOM)).
series("x", Plot.data().xy(dTs, dX), Plot.seriesOpts().color(Color.RED)).
series("y", Plot.data().xy(dTs, dY), Plot.seriesOpts().color(Color.BLUE)).
series("z", Plot.data().xy(dTs, dZ), Plot.seriesOpts().color(Color.GREEN));
series("x", Plot.data().xy(dTsInterp, acInterp.getX()), Plot.seriesOpts().color(Color.RED)).
series("y", Plot.data().xy(dTsInterp, acInterp.getY()), Plot.seriesOpts().color(Color.BLUE)).
series("z", Plot.data().xy(dTsInterp, acInterp.getZ()), Plot.seriesOpts().color(Color.GREEN));
windowRaw.set(plotRaw.draw());
//auto corr
float[] xAutoCorr = Utils.fftAutoCorrelation(accWindowBuffer.getX(), 1024);
float[] yAutoCorr = Utils.fftAutoCorrelation(accWindowBuffer.getY(), 1024);
float[] zAutoCorr = Utils.fftAutoCorrelation(accWindowBuffer.getZ(), 1024);
double[] xAutoCorr = new AutoCorrelation(acInterp.getX(), 1024).getCorr();
double[] yAutoCorr = new AutoCorrelation(acInterp.getY(), 1024).getCorr();
double[] zAutoCorr = new AutoCorrelation(acInterp.getZ(), 1024).getCorr();
//print autocorr
int[] tmp = IntStream.rangeClosed(-((xAutoCorr.length - 1)/2), ((xAutoCorr.length - 1)/2)).toArray();
double[] rangeAuto = IntStream.range(0, tmp.length).mapToDouble(i -> tmp[i]).toArray();
double[] dXAuto = IntStream.range(0, xAutoCorr.length).mapToDouble(i -> xAutoCorr[i]).toArray();
double[] dYAuto = IntStream.range(0, yAutoCorr.length).mapToDouble(i -> yAutoCorr[i]).toArray();
double[] dZAuto = IntStream.range(0, zAutoCorr.length).mapToDouble(i -> zAutoCorr[i]).toArray();
Plot plotCorr = Plot.plot(Plot.plotOpts().
title("Auto Correlation").
legend(Plot.LegendFormat.BOTTOM)).
series("x", Plot.data().xy(rangeAuto, dXAuto), Plot.seriesOpts().color(Color.RED)).
series("y", Plot.data().xy(rangeAuto, dYAuto), Plot.seriesOpts().color(Color.BLUE)).
series("z", Plot.data().xy(rangeAuto, dZAuto), Plot.seriesOpts().color(Color.GREEN));
series("x", Plot.data().xy(rangeAuto, xAutoCorr), Plot.seriesOpts().color(Color.RED)).
series("y", Plot.data().xy(rangeAuto, yAutoCorr), Plot.seriesOpts().color(Color.BLUE)).
series("z", Plot.data().xy(rangeAuto, zAutoCorr), Plot.seriesOpts().color(Color.GREEN));
windowAuto.set(plotCorr.draw());
//find peaks
LinkedList<Integer> peaksX = Utils.findPeaks(xAutoCorr, 50, 0.1f, 0, false);
Peaks pX = new Peaks(xAutoCorr, 50, 0.1f, 0, false);
LinkedList<Integer> peaksX = pX.getPeaksIdx();
double[] dPeaksXX = IntStream.range(0, peaksX.size()).mapToDouble(i -> (peaksX.get(i) - 1024)).toArray();//peaks.stream().mapToDouble(i->i).toArray();
double[] dPeaksXY = IntStream.range(0, peaksX.size()).mapToDouble(i -> (dXAuto[peaksX.get(i)])).toArray();
double[] dPeaksXY = IntStream.range(0, peaksX.size()).mapToDouble(i -> (xAutoCorr[peaksX.get(i)])).toArray();
Plot plotPeaksX = Plot.plot(Plot.plotOpts().
title("Peak Detection on X").
legend(Plot.LegendFormat.BOTTOM)).
series("x", Plot.data().xy(rangeAuto, dXAuto), Plot.seriesOpts().color(Color.RED)).
series("x", Plot.data().xy(rangeAuto, xAutoCorr), Plot.seriesOpts().color(Color.RED)).
series("Peaks", Plot.data().xy(dPeaksXX, dPeaksXY), Plot.seriesOpts().color(Color.CYAN).
marker(Plot.Marker.DIAMOND).line(Plot.Line.NONE));
windowPeaksX.set(plotPeaksX.draw());
LinkedList<Integer> peaksY = Utils.findPeaks(yAutoCorr, 50, 0.1f, 0, false);
Peaks pY = new Peaks(yAutoCorr, 50, 0.1f, 0, false);
LinkedList<Integer> peaksY = pY.getPeaksIdx();
double[] dPeaksYX = IntStream.range(0, peaksY.size()).mapToDouble(i -> (peaksY.get(i) - 1024)).toArray();//peaks.stream().mapToDouble(i->i).toArray();
double[] dPeaksYY = IntStream.range(0, peaksY.size()).mapToDouble(i -> (dYAuto[peaksY.get(i)])).toArray();
double[] dPeaksYY = IntStream.range(0, peaksY.size()).mapToDouble(i -> (yAutoCorr[peaksY.get(i)])).toArray();
Plot plotPeaksY = Plot.plot(Plot.plotOpts().
title("Peak Detection on Y").
legend(Plot.LegendFormat.BOTTOM)).
series("x", Plot.data().xy(rangeAuto, dYAuto), Plot.seriesOpts().color(Color.RED)).
series("x", Plot.data().xy(rangeAuto, yAutoCorr), Plot.seriesOpts().color(Color.RED)).
series("Peaks", Plot.data().xy(dPeaksYX, dPeaksYY), Plot.seriesOpts().color(Color.CYAN).
marker(Plot.Marker.DIAMOND).line(Plot.Line.NONE));
windowPeaksY.set(plotPeaksY.draw());
LinkedList<Integer> peaksZ = Utils.findPeaks(zAutoCorr, 50, 0.1f, 0, false);
Peaks pZ = new Peaks(zAutoCorr, 50, 0.1f, 0, false);
LinkedList<Integer> peaksZ = pZ.getPeaksIdx();
double[] dPeaksZX = IntStream.range(0, peaksZ.size()).mapToDouble(i -> (peaksZ.get(i) - 1024)).toArray();//peaks.stream().mapToDouble(i->i).toArray();
double[] dPeaksZY = IntStream.range(0, peaksZ.size()).mapToDouble(i -> (dZAuto[peaksZ.get(i)])).toArray();
double[] dPeaksZY = IntStream.range(0, peaksZ.size()).mapToDouble(i -> (zAutoCorr[peaksZ.get(i)])).toArray();
Plot plotPeaksZ = Plot.plot(Plot.plotOpts().
title("Peak Detection on Z").
legend(Plot.LegendFormat.BOTTOM)).
series("x", Plot.data().xy(rangeAuto, dZAuto), Plot.seriesOpts().color(Color.RED)).
series("x", Plot.data().xy(rangeAuto, zAutoCorr), Plot.seriesOpts().color(Color.RED)).
series("Peaks", Plot.data().xy(dPeaksZX, dPeaksZY), Plot.seriesOpts().color(Color.CYAN).
marker(Plot.Marker.DIAMOND).line(Plot.Line.NONE));
@@ -124,6 +123,13 @@ public class Main {
//fill hols improve peaks
//estimate bpm between detected peaks
System.out.println("BPM-X: " + pX.getBPM(6));
System.out.println("BPM-Y: " + pY.getBPM(6));
System.out.println("BPM-Z: " + pZ.getBPM(6));
//todo: statistikzeuch und mit matlab vergleichen.
//todo: todos machen. lol
//todo: kleiner fenstergrößen testen. so ist doch etwas langsam auf der Uhr.
int dummyForBreakpoint = 0;
}

272
java/src/main/java/Utils.java
View File

@@ -1,272 +0,0 @@
import org.jtransforms.fft.FloatFFT_1D;
import javax.swing.*;
import java.awt.*;
import java.awt.image.BufferedImage;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.LinkedList;
public class Utils {
public static float getDistance(float x1, float y1, float x2, float y2) {
return (float) Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
}
public static class AccelerometerData {
public float x,y,z;
public long ts;
public AccelerometerData(long ts, float x, float y, float z){
this.ts = ts;
this.x = x;
this.y = y;
this.z = z;
}
}
public static class AccelerometerWindowBuffer extends ArrayList<AccelerometerData> {
private int mWindowSize;
private int mOverlapSize;
private int mOverlapCounter;
private float[] mX;
private float[] mY;
private float[] mZ;
private long[] mTs;
public AccelerometerWindowBuffer(int windowSize, int overlap){
this.mWindowSize = windowSize;
this.mOverlapSize = overlap;
mOverlapCounter = 1;
mX = new float[this.mWindowSize];
mY = new float[this.mWindowSize];
mZ = new float[this.mWindowSize];
mTs = new long[this.mWindowSize];
}
public boolean add(AccelerometerData ad){
boolean r = super.add(ad);
if (size() > mWindowSize){
removeRange(0, size() - mWindowSize);
}
//update the double arrays.
for (int i = 0; i < size(); ++i) {
mX[i] = get(i).x;
mY[i] = get(i).y;
mZ[i] = get(i).z;
mTs[i] = get(i).ts;
}
++mOverlapCounter;
return r;
}
public boolean isNextWindowReady(){
if(size() == mWindowSize && mOverlapCounter > mOverlapSize){
mOverlapCounter = 1;
return true;
}
return false;
}
public AccelerometerData getYongest() {
return get(size() - 1);
}
public AccelerometerData getOldest() {
return get(0);
}
public float[] getX(){
return mX;
}
public float[] getY(){
return mY;
}
public float[] getZ(){
return mZ;
}
public long[] getTs(){
return mTs;
}
}
public static float sqr(float x) {
return x * x;
}
public static int nextPow2(int a){
return a == 0 ? 0 : 32 - Integer.numberOfLeadingZeros(a - 1);
}
public static float mean(float[] data){
float sum = 0;
for (int i = 0; i < data.length; i++) {
sum += data[i];
}
return sum / data.length;
}
public static float[] removeZero(float[] array){
int j = 0;
for( int i=0; i<array.length; i++ )
{
if (array[i] != 0)
array[j++] = array[i];
}
float[] newArray = new float[j];
System.arraycopy( array, 0, newArray, 0, j );
return newArray;
}
//TODO: errorhandling maxLag = 0;
//TODO: größeren Testcase schreiben
public static float[] fftAutoCorrelation(float[] data, int maxLag) {
int n = data.length;
float[] x = Arrays.copyOf(data, n);
int mxl = Math.min(maxLag, n - 1);
int ceilLog2 = nextPow2(2*n -1);
int n2 = (int) Math.pow(2,ceilLog2);
// x - mean(x) (pointwise)
float x_mean = mean(x);
for(int i = 0; i < x.length; ++i){
x[i] -= x_mean;
}
// double the size of x and fill up with zeros. if x is not even, add additional 0
float[] x2 = new float[n2 * 2]; //need double the size for fft.realForwardFull (look into method description)
Arrays.fill(x2, 0);
System.arraycopy(x,0, x2, 0, x.length);
// x_fft calculate fft 1D
FloatFFT_1D fft = new FloatFFT_1D(n2);
fft.realForwardFull(x2);
// Cr = abs(x_fft).^2 (absolute with complex numbers is (r^2) + (i^2)
float[] Cr = new float[n2 * 2];
int j = 0;
for(int i = 0; i < x2.length; ++i){
Cr[j++] = sqr(x2[i]) + sqr(x2[i+1]);
++i; //skip the complex part
}
// ifft(Cr,[],1)
FloatFFT_1D ifft = new FloatFFT_1D(n2);
ifft.realInverseFull(Cr, true);
// remove complex part and scale/normalize
float[] c1 = new float[n2];
j = 0;
for(int i = 0; i < Cr.length; ++i){
c1[j++] = Cr[i] / Cr[0];
++i; //skip the complex part
}
// Keep only the lags we want and move negative lags before positive lags.
float[] c = new float[(mxl * 2) + 1];
System.arraycopy(c1, 0, c, mxl, mxl + 1); // +1 to place the 1.0 in the middle of correlation
System.arraycopy(c1, n2 - mxl, c, 0, mxl);
return c;
}
//TODO: findPeaks method identical to Matlab... with PeakProminence
/**
* Simple method for finding local maxima in an array
* @param data input
* @param width minimum distance between peaks
* @param threshold minimum value of peaks
* @param decayRate how quickly previous peaks are forgotten
* @param isRelative minimum value of peaks is relative to local average
* @return array of peaks
*/
public static LinkedList<Integer> findPeaks(float[] data, int width, float threshold, float decayRate, boolean isRelative) {
LinkedList<Integer> peaks = new LinkedList<>();
int maxp = 0;
int mid = 0;
int end = data.length;
float av = data[0];
while (mid < end) {
av = decayRate * av + (1 - decayRate) * data[mid];
if (av < data[mid])
av = data[mid];
int i = mid - width;
if (i < 0)
i = 0;
int stop = mid + width + 1;
if (stop > data.length)
stop = data.length;
maxp = i;
for (i++; i < stop; i++)
if (data[i] > data[maxp])
maxp = i;
if (maxp == mid) {
if (overThreshold(data, maxp, width, threshold, isRelative, av)){
peaks.add(new Integer(maxp));
}
}
mid++;
}
return peaks;
}
private static boolean overThreshold(float[] data, int index, int width,
float threshold, boolean isRelative,
float av) {
int pre = 3;
int post = 1;
if (data[index] < av)
return false;
if (isRelative) {
int iStart = index - pre * width;
if (iStart < 0)
iStart = 0;
int iStop = index + post * width;
if (iStop > data.length)
iStop = data.length;
float sum = 0;
int count = iStop - iStart;
while (iStart < iStop)
sum += data[iStart++];
return (data[index] > sum / count + threshold);
} else
return (data[index] > threshold);
}
@SuppressWarnings("serial")
public static class ShowPNG extends JFrame
{
JLabel mLabel;
ImageIcon mIcon;
public ShowPNG(){
mLabel = new JLabel();
this.setLayout(new GridLayout(1,1));
this.setSize(800, 640);
this.add(mLabel);
this.setVisible(true);
}
public void set(BufferedImage bi){
mIcon = new ImageIcon(bi);
mLabel.setVisible(false);
mLabel.setIcon(mIcon);
mLabel.setVisible(true);
}
}
}

31
java/src/main/java/bpmEstimation/AccelerometerData.java Normal file
View File

@@ -0,0 +1,31 @@
package bpmEstimation;
/**
* Created by toni on 15/12/17.
*/
public class AccelerometerData {
public double x,y,z;
public long ts;
public AccelerometerData(long ts, double x, double y, double z){
this.ts = ts;
this.x = x;
this.y = y;
this.z = z;
}
@Override
public boolean equals(Object other){
if (this == other)
return true;
if (!(other instanceof AccelerometerData)) {
return false;
}
AccelerometerData ad = (AccelerometerData) other;
return ts == ad.ts;
}
}

114
java/src/main/java/bpmEstimation/AccelerometerInterpolator.java Normal file
View File

@@ -0,0 +1,114 @@
package bpmEstimation;
import utilities.Utils;
import java.util.Arrays;
/**
* Created by toni on 16/12/17.
*/
public class AccelerometerInterpolator {
private double[] mX;
private double[] mY;
private double[] mZ;
private long[] mTsInterp;
public AccelerometerInterpolator(AccelerometerWindowBuffer ab, double sampleRate_ms){
if(sampleRate_ms <= 0){
sampleRate_ms = Utils.mean(Utils.diff(ab.getTs()));
}
long size = (ab.getYongest().ts - (ab.getOldest().ts - (long) sampleRate_ms)) / (long) sampleRate_ms;
mTsInterp = new long[(int)size];
int j = 0;
for(long i = ab.getOldest().ts; i <= ab.getYongest().ts; i += sampleRate_ms){
mTsInterp[j++] = i;
}
mX = interpLinear(ab.getTs(), ab.getX(), mTsInterp);
mY = interpLinear(ab.getTs(), ab.getY(), mTsInterp);
mZ = interpLinear(ab.getTs(), ab.getZ(), mTsInterp);
}
public long[] getTs(){
return mTsInterp;
}
public double[] getX(){
return mX;
}
public double[] getY(){
return mY;
}
public double[] getZ(){
return mZ;
}
public static double[] interpLinear(double[] x, double[] y, double[] xi) throws IllegalArgumentException {
if (x.length != y.length) {
throw new IllegalArgumentException("X and Y must be the same length");
}
if (x.length == 1) {
throw new IllegalArgumentException("X must contain more than one value");
}
double[] dx = new double[x.length - 1];
double[] dy = new double[x.length - 1];
double[] slope = new double[x.length - 1];
double[] intercept = new double[x.length - 1];
// Calculate the line equation (i.e. slope and intercept) between each point
for (int i = 0; i < x.length - 1; i++) {
dx[i] = x[i + 1] - x[i];
if (dx[i] == 0) {
throw new IllegalArgumentException("X must be montotonic. A duplicate " + "x-value was found");
}
if (dx[i] < 0) {
throw new IllegalArgumentException("X must be sorted");
}
dy[i] = y[i + 1] - y[i];
slope[i] = dy[i] / dx[i];
intercept[i] = y[i] - x[i] * slope[i];
}
// Perform the interpolation here
double[] yi = new double[xi.length];
for (int i = 0; i < xi.length; i++) {
if ((xi[i] > x[x.length - 1]) || (xi[i] < x[0])) {
yi[i] = Double.NaN;
}
else {
int loc = Arrays.binarySearch(x, xi[i]);
if (loc < -1) {
loc = -loc - 2;
yi[i] = slope[loc] * xi[i] + intercept[loc];
}
else {
yi[i] = y[loc];
}
}
}
return yi;
}
public static double[] interpLinear(long[] x, double[] y, long[] xi) throws IllegalArgumentException {
double[] xd = new double[x.length];
for (int i = 0; i < x.length; i++) {
xd[i] = x[i];
}
double[] xid = new double[xi.length];
for (int i = 0; i < xi.length; i++) {
xid[i] = xi[i];
}
return interpLinear(xd, y, xid);
}
}

89
java/src/main/java/bpmEstimation/AccelerometerWindowBuffer.java Normal file
View File

@@ -0,0 +1,89 @@
package bpmEstimation;
import utilities.Utils;
import java.util.ArrayList;
/**
* Created by toni on 15/12/17.
*/
public class AccelerometerWindowBuffer extends ArrayList<AccelerometerData> {
private int mWindowSize;
private int mOverlapSize;
private int mOverlapCounter;
private double[] mX;
private double[] mY;
private double[] mZ;
private long[] mTs;
public AccelerometerWindowBuffer(int windowSize, int overlap){
this.mWindowSize = windowSize;
this.mOverlapSize = overlap;
mOverlapCounter = 1;
mX = new double[this.mWindowSize];
mY = new double[this.mWindowSize];
mZ = new double[this.mWindowSize];
mTs = new long[this.mWindowSize];
}
//TODO: add exception handling. falseArgument if ad has no numeric x,y,z
public boolean add(AccelerometerData ad){
//do not add duplicates!
if(!isEmpty() && getYongest().equals(ad)){
return false;
}
boolean r = super.add(ad);
if (size() > mWindowSize){
removeRange(0, size() - mWindowSize);
}
//TODO: bullshit... jedes mal wieder alle neue anlegen ist ja nur langsam
//update the double arrays.
for (int i = 0; i < size(); ++i) {
mX[i] = get(i).x;
mY[i] = get(i).y;
mZ[i] = get(i).z;
mTs[i] = get(i).ts;
}
++mOverlapCounter;
return r;
}
public boolean isNextWindowReady(){
if(size() == mWindowSize && mOverlapCounter > mOverlapSize){
mOverlapCounter = 1;
return true;
}
return false;
}
public AccelerometerData getYongest() {
return get(size() - 1);
}
public AccelerometerData getOldest() {
return get(0);
}
public double[] getX(){
return mX;
}
public double[] getY(){
return mY;
}
public double[] getZ(){
return mZ;
}
public long[] getTs(){
return mTs;
}
}

79
java/src/main/java/bpmEstimation/AutoCorrelation.java Normal file
View File

@@ -0,0 +1,79 @@
package bpmEstimation;
import org.jtransforms.fft.DoubleFFT_1D;
import utilities.Utils;
import java.util.Arrays;
/**
* Created by toni on 15/12/17.
*/
public class AutoCorrelation {
int mMaxLag;
double[] mCorr;
public AutoCorrelation(double[] data, int maxLag){
mMaxLag = maxLag;
mCorr = fft(data);
}
public double[] getCorr(){
return mCorr;
}
private double[] fft(double[] data) {
if(mMaxLag < 1){
throw new RuntimeException("maxlag has to be greater 1");
}
int n = data.length;
double[] x = Arrays.copyOf(data, n);
int mxl = Math.min(mMaxLag, n - 1);
int ceilLog2 = Utils.nextPow2(2*n -1);
int n2 = (int) Math.pow(2,ceilLog2);
// x - mean(x) (pointwise)
double x_mean = Utils.mean(x);
for(int i = 0; i < x.length; ++i){
x[i] -= x_mean;
}
// double the size of x and fill up with zeros. if x is not even, add additional 0
double[] x2 = new double[n2 * 2]; //need double the size for fft.realForwardFull (look into method description)
Arrays.fill(x2, 0);
System.arraycopy(x,0, x2, 0, x.length);
// x_fft calculate fft 1D
DoubleFFT_1D fft = new DoubleFFT_1D(n2);
fft.realForwardFull(x2);
// Cr = abs(x_fft).^2 (absolute with complex numbers is (r^2) + (i^2)
double[] Cr = new double[n2 * 2];
int j = 0;
for(int i = 0; i < x2.length; ++i){
Cr[j++] = Utils.sqr(x2[i]) + Utils.sqr(x2[i+1]);
++i; //skip the complex part
}
// ifft(Cr,[],1)
DoubleFFT_1D ifft = new DoubleFFT_1D(n2);
ifft.realInverseFull(Cr, true);
// remove complex part and scale/normalize
double[] c1 = new double[n2];
j = 0;
for(int i = 0; i < Cr.length; ++i){
c1[j++] = Cr[i] / Cr[0];
++i; //skip the complex part
}
// Keep only the lags we want and move negative lags before positive lags.
double[] c = new double[(mxl * 2) + 1];
System.arraycopy(c1, 0, c, mxl, mxl + 1); // +1 to place the 1.0 in the middle of correlation
System.arraycopy(c1, n2 - mxl, c, 0, mxl);
return c;
}
}

207
java/src/main/java/bpmEstimation/Peaks.java Normal file
View File

@@ -0,0 +1,207 @@
package bpmEstimation;
import utilities.Utils;
import java.util.LinkedList;
/**
* Created by toni on 15/12/17.
*/
public class Peaks {
private LinkedList<Integer> mPeaksIdx; //provide the idx within the given data array
private LinkedList<Integer> mPeaksPos; //the real position within the data-rang e.g. lag -1024 to 1024
private LinkedList<Double> mPeaksValue; //the value at mPeaksPos
private double[] mData;
/**
* Simple method for finding local maxima in an array
* @param data input
* @param width minimum distance between peaks
* @param threshold minimum value of peaks
* @param decayRate how quickly previous peaks are forgotten
* @param isRelative minimum value of peaks is relative to local average
* @return array of peaks
*/
public Peaks(double[] data, int width, double threshold, double decayRate, boolean isRelative){
this.mData = data;
this.mPeaksIdx = new LinkedList<>();
this.mPeaksPos = new LinkedList<>();
this.mPeaksValue = new LinkedList<>();
//create the peaks
simplePeakFinder(data, width, threshold, decayRate, isRelative);
updateLists();
}
public LinkedList<Integer> getPeaksIdx() {
return mPeaksIdx;
}
public int[] getPeaksIdxAsArray() {
return mPeaksIdx.stream().mapToInt(i->i).toArray();
}
public LinkedList<Integer> getPeaksPos() {
return mPeaksPos;
}
public double[] getPeaksPosAsArray() {
return mPeaksPos.stream().mapToDouble(i -> i).toArray();
}
public LinkedList<Double> getPeaksValue() {
return mPeaksValue;
}
public double[] getPeaksValueAsArray() {
return mPeaksValue.stream().mapToDouble(i -> i).toArray();
}
public double[] getData(){
return mData;
}
//TODO: implement findFalseDetectedPeaks
public void improveResults(double[] data){
updateLists();
}
public double[] getPeaksValueWithoutZeroIdx() {
double[] values = new double[mPeaksIdx.size() - 1];
int i = 0;
int mid = (mData.length / 2);
for(Integer idx : mPeaksIdx) {
if(!idx.equals(mid)){
values[i] = mPeaksValue.get(idx);
++i;
}
}
return values;
}
public double[] getPeaksValueWithoutNegativeValues() {
double[] values = new double[mPeaksIdx.size() - 1];
int i = 0;
for(Integer idx : mPeaksIdx) {
double curVal = mPeaksValue.get(idx);
if(curVal > 0){
values[i] = curVal;
++i;
}
}
return values;
}
public double[] getPeaksValueWithoutZeroIdxAndNegativeValues(){
double[] values = new double[mPeaksIdx.size() - 1];
int i = 0;
int mid = (mData.length / 2);
for(Integer idx : mPeaksIdx) {
double curVal = mPeaksValue.get(idx);
if(!idx.equals(mid) && curVal> 0){
values[i] = curVal;
++i;
}
}
return values;
}
public boolean hasPeaks() {
return mPeaksIdx.size() > 1;
}
/**
* Provides an estimation of beats per minute given a samplerate in milliseconds
* @param sampleRate_ms
* @return bpm if peaks found and conducting activity recognized, else -1
*/
public double getBPM(double sampleRate_ms){
//todo: rückweisungsklasse kann auch hier mit rein.
if(hasPeaks()){
//todo diff and mean method for linkedlists for speed
return 60000 / Utils.mean(Utils.diff(mPeaksPos.stream().mapToDouble(i -> i * sampleRate_ms).toArray()));
}
return -1;
}
/**
* updates the position and values of the found peaks.
* call this if peaks are somewhat changed.
*/
private void updateLists(){
//fill the position and the value lists
for(Integer idx : mPeaksIdx){
int mid = (mData.length / 2);
mPeaksPos.add(idx - mid);
mPeaksValue.add(mData[idx]);
}
}
//TODO: findPeaks method identical to Matlab... with PeakProminence
private void simplePeakFinder(double[] data, int width, double threshold, double decayRate, boolean isRelative) {
int maxp = 0;
int mid = 0;
int end = data.length;
double av = data[0];
while (mid < end) {
av = decayRate * av + (1 - decayRate) * data[mid];
if (av < data[mid])
av = data[mid];
int i = mid - width;
if (i < 0)
i = 0;
int stop = mid + width + 1;
if (stop > data.length)
stop = data.length;
maxp = i;
for (i++; i < stop; i++)
if (data[i] > data[maxp])
maxp = i;
if (maxp == mid) {
if (overThreshold(data, maxp, width, threshold, isRelative, av)){
this.mPeaksIdx.add(new Integer(maxp));
}
}
mid++;
}
}
private boolean overThreshold(double[] data, int index, int width,
double threshold, boolean isRelative,
double av) {
int pre = 3;
int post = 1;
if (data[index] < av)
return false;
if (isRelative) {
int iStart = index - pre * width;
if (iStart < 0)
iStart = 0;
int iStop = index + post * width;
if (iStop > data.length)
iStop = data.length;
double sum = 0;
int count = iStop - iStart;
while (iStart < iStop)
sum += data[iStart++];
return (data[index] > sum / count + threshold);
} else
return (data[index] > threshold);
}
}

232
java/src/main/java/utilities/Utils.java Normal file
View File

@@ -0,0 +1,232 @@
package utilities;
import bpmEstimation.Peaks;
import javax.swing.*;
import java.awt.*;
import java.awt.image.BufferedImage;
import java.util.Arrays;
//TODO: change from double to generic type
public class Utils {
public static double getDistance(double x1, double y1, double x2, double y2) {
return (double) Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
}
public static double sqr(double x) {
return x * x;
}
public static int nextPow2(int a){
return a == 0 ? 0 : 32 - Integer.numberOfLeadingZeros(a - 1);
}
public static double sum(double[] data){
double sum = 0;
for (int i = 0; i < data.length; i++) {
sum += data[i];
}
return sum;
}
public static long sum(long[] data){
long sum = 0;
for (int i = 0; i < data.length; i++) {
sum += data[i];
}
return sum;
}
//TODO: Could be slow.. faster method?
public static double rms(double[] nums) {
double sum = 0.0f;
for (double num : nums)
sum += num * num;
return (double) Math.sqrt(sum / nums.length);
}
public static double[] diff(double[] data){
double[] diff = new double[data.length - 1];
int i=0;
for(int j = 1; j < data.length; ++j){
diff[i] = data[j] - data[i];
++i;
}
return diff;
}
public static long[] diff(long[] data){
long[] diff = new long[data.length - 1];
int i=0;
for(int j = 1; j < data.length; ++j){
diff[i] = data[j] - data[i];
++i;
}
return diff;
}
public static double mean(double[] data){
return sum(data) / data.length;
}
public static long mean(long[] data){
return sum(data) / data.length;
}
//TODO: Could be slow.. faster method?
public static double geometricMean(double[] data) {
double sum = data[0];
for (int i = 1; i < data.length; i++) {
sum *= data[i];
}
return (double) Math.pow(sum, 1.0 / data.length);
}
public static int intersectionNumber(double[] signal, double border){
int cnt = 0;
boolean isSmallerValue = false;
boolean isBiggerValue = false;
for(double value : signal){
if(value < border){
if(isBiggerValue){
cnt++;
}
isSmallerValue = true;
isBiggerValue = false;
}
else {
if(isSmallerValue){
cnt++;
}
isSmallerValue = false;
isBiggerValue = true;
}
}
return cnt;
}
public static double getBestBpmEstimation(Peaks peaksX, Peaks peaksY, Peaks peaksZ) throws IllegalArgumentException {
int cntNumAxis = 0;
double sumCorr = 1; //to prevent division by zero
double sumRms = 1;
double sumNumInter = 1;
double corrMeanX = 0, corrRmsX = 0;
int corrNumInterX = 0;
if(peaksX.hasPeaks()){
corrMeanX = geometricMean(peaksX.getPeaksValueWithoutZeroIdxAndNegativeValues());
corrRmsX = rms(peaksX.getPeaksValueWithoutZeroIdx());
corrNumInterX = intersectionNumber(peaksX.getData(), 0.2f);
++cntNumAxis;
sumCorr += corrMeanX;
sumRms += corrRmsX;
sumNumInter += corrNumInterX;
}
double corrMeanY = 0, corrRmsY = 0;
int corrNumInterY = 0;
if(peaksY.hasPeaks()){
corrMeanY = geometricMean(peaksY.getPeaksValueWithoutZeroIdxAndNegativeValues());
corrRmsY = rms(peaksY.getPeaksValueWithoutZeroIdx());
corrNumInterY = intersectionNumber(peaksY.getData(), 0.2f);
++cntNumAxis;
sumCorr += corrMeanY;
sumRms += corrRmsY;
sumNumInter += corrNumInterY;
}
double corrMeanZ = 0, corrRmsZ = 0;
int corrNumInterZ = 0;
if(peaksZ.hasPeaks()){
corrMeanZ = geometricMean(peaksZ.getPeaksValueWithoutZeroIdxAndNegativeValues());
corrRmsZ = rms(peaksZ.getPeaksValueWithoutZeroIdx());
corrNumInterZ = intersectionNumber(peaksZ.getData(), 0.2f);
++cntNumAxis;
sumCorr += corrMeanZ;
sumRms += corrRmsZ;
sumNumInter += corrNumInterZ;
}
if(cntNumAxis == 0){
throw new IllegalArgumentException("All Peaks are empty! -> Reject Estimation");
}
double quantityX = ((corrMeanX / sumCorr) + (corrRmsX / sumRms) + (corrNumInterX / sumNumInter)) / cntNumAxis;
double quantityY = ((corrMeanY / sumCorr) + (corrRmsY / sumRms) + (corrNumInterY / sumNumInter)) / cntNumAxis;
double quantityZ = ((corrMeanZ / sumCorr) + (corrRmsZ / sumRms) + (corrNumInterZ / sumNumInter)) / cntNumAxis;
//get best axis by quantity and estimate bpm
if(quantityX > quantityY && quantityX > quantityZ){
}
else if(quantityY > quantityZ){
}
else {
}
//check if values are empty.
return 0.0f;
}
public static double[] removeZero(double[] array){
int j = 0;
for( int i=0; i<array.length; i++ )
{
if (array[i] != 0){
array[j++] = array[i];
}
}
double[] newArray = new double[j];
System.arraycopy( array, 0, newArray, 0, j );
return newArray;
}
public static double[] greaterZero(double[] array){
int j = 0;
for( int i=0; i < array.length; i++ )
{
if (array[i] > 0){
array[j++] = array[i];
}
}
double[] newArray = new double[j];
System.arraycopy( array, 0, newArray, 0, j );
return newArray;
}
@SuppressWarnings("serial")
public static class ShowPNG extends JFrame
{
JLabel mLabel;
ImageIcon mIcon;
public ShowPNG(){
mLabel = new JLabel();
this.setLayout(new GridLayout(1,1));
this.setSize(800, 640);
this.add(mLabel);
this.setVisible(true);
}
public void set(BufferedImage bi){
mIcon = new ImageIcon(bi);
mLabel.setVisible(false);
mLabel.setIcon(mIcon);
mLabel.setVisible(true);
}
}
}