david/hsh-prevent
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

initial

Browse Source
This commit is contained in:
David Madl
2026-03-02 06:19:34 +01:00
commit d83a7506d6
273 changed files with 364230 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
signal/CMakeLists.txt Normal file
View File

@@ -0,0 +1,44 @@
# For more information about using CMake with Android Studio, read the
# documentation: https://d.android.com/studio/projects/add-native-code.html
# Sets the minimum version of CMake required to build the native library.
cmake_minimum_required(VERSION 3.4.1)
# Creates and names a library, sets it as either STATIC
# or SHARED, and provides the relative paths to its source code.
# You can define multiple libraries, and CMake builds them for you.
# Gradle automatically packages shared libraries with your APK.
add_library( # Sets the name of the library.
native-lib
# Sets the library as a shared library.
SHARED
# Provides a relative path to your source file(s).
src/main/cpp/native-lib.cpp )
# Searches for a specified prebuilt library and stores the path as a
# variable. Because CMake includes system libraries in the search path by
# default, you only need to specify the name of the public NDK library
# you want to add. CMake verifies that the library exists before
# completing its build.
find_library( # Sets the name of the path variable.
log-lib
# Specifies the name of the NDK library that
# you want CMake to locate.
log )
# Specifies libraries CMake should link to your target library. You
# can link multiple libraries, such as libraries you define in this
# build script, prebuilt third-party libraries, or system libraries.
target_link_libraries( # Specifies the target library.
native-lib
# Links the target library to the log library
# included in the NDK.
${log-lib} )

55
signal/build.gradle Normal file
View File

@@ -0,0 +1,55 @@
apply plugin: 'com.android.library'
android {
compileSdkVersion 25
buildToolsVersion "25.0.2"
defaultConfig {
minSdkVersion 23
targetSdkVersion 25
versionCode 1
versionName "1.0"
testInstrumentationRunner "android.support.test.runner.AndroidJUnitRunner"
/*externalNativeBuild {
cmake {
cppFlags "-std=c++11"
}
}*/
}
buildTypes {
release {
minifyEnabled false
proguardFiles getDefaultProguardFile('proguard-android.txt'), 'proguard-rules.pro'
}
}
/*externalNativeBuild {
cmake {
path "CMakeLists.txt"
}
}*/
}
allprojects {
repositories {
maven { url "https://jitpack.io" } // for iirj
}
}
dependencies {
compile fileTree(dir: 'libs', include: ['*.jar'])
androidTestCompile('com.android.support.test.espresso:espresso-core:2.2.2', {
exclude group: 'com.android.support', module: 'support-annotations'
})
compile 'com.android.support:appcompat-v7:25.3.0'
compile group: 'org.apache.commons', name: 'commons-math3', version: '3.6.1'
testCompile 'junit:junit:4.12'
testCompile 'org.jfree:jfreechart:1.0.19'
//compile group: 'uk.me.berndporr', name:'iirj', version: '1.0'
compile 'com.github.berndporr:iirj:v1.0.4'
compile group: 'org.apache.commons', name: 'commons-math3', version: '3.6.1'
testCompile 'org.jfree:jcommon:1.0.24'
testCompile 'org.json:json:20160810'
testCompile files('/usr/lib/jvm/java-8-openjdk/jre/lib/rt.jar') // for JFreeChart on PC, need javax.swing etc.
}

25
signal/proguard-rules.pro vendored Normal file
View File

@@ -0,0 +1,25 @@
# Add project specific ProGuard rules here.
# By default, the flags in this file are appended to flags specified
# in /opt/android-sdk/tools/proguard/proguard-android.txt
# You can edit the include path and order by changing the proguardFiles
# directive in build.gradle.
#
# For more details, see
# http://developer.android.com/guide/developing/tools/proguard.html
# Add any project specific keep options here:
# If your project uses WebView with JS, uncomment the following
# and specify the fully qualified class name to the JavaScript interface
# class:
#-keepclassmembers class fqcn.of.javascript.interface.for.webview {
# public *;
#}
# Uncomment this to preserve the line number information for
# debugging stack traces.
#-keepattributes SourceFile,LineNumberTable
# If you keep the line number information, uncomment this to
# hide the original source file name.
#-renamesourcefileattribute SourceFile

128
signal/signal.iml Normal file
View File

@@ -0,0 +1,128 @@
<?xml version="1.0" encoding="UTF-8"?>
<module external.linked.project.id=":signal" external.linked.project.path="$MODULE_DIR$" external.root.project.path="$MODULE_DIR$/.." external.system.id="GRADLE" type="JAVA_MODULE" version="4">
<component name="FacetManager">
<facet type="android-gradle" name="Android-Gradle">
<configuration>
<option name="GRADLE_PROJECT_PATH" value=":signal" />
</configuration>
</facet>
<facet type="android" name="Android">
<configuration>
<option name="SELECTED_BUILD_VARIANT" value="debug" />
<option name="ASSEMBLE_TASK_NAME" value="assembleDebug" />
<option name="COMPILE_JAVA_TASK_NAME" value="compileDebugSources" />
<afterSyncTasks>
<task>generateDebugSources</task>
</afterSyncTasks>
<option name="ALLOW_USER_CONFIGURATION" value="false" />
<option name="MANIFEST_FILE_RELATIVE_PATH" value="/src/main/AndroidManifest.xml" />
<option name="RES_FOLDER_RELATIVE_PATH" value="/src/main/res" />
<option name="RES_FOLDERS_RELATIVE_PATH" value="file://$MODULE_DIR$/src/main/res" />
<option name="ASSETS_FOLDER_RELATIVE_PATH" value="/src/main/assets" />
<option name="PROJECT_TYPE" value="1" />
</configuration>
</facet>
</component>
<component name="NewModuleRootManager" LANGUAGE_LEVEL="JDK_1_7" inherit-compiler-output="false">
<output url="file://$MODULE_DIR$/build/intermediates/classes/debug" />
<output-test url="file://$MODULE_DIR$/build/intermediates/classes/test/debug" />
<exclude-output />
<content url="file://$MODULE_DIR$">
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/r/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/aidl/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/buildConfig/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/rs/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/apt/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/rs/debug" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/resValues/debug" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/r/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/aidl/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/buildConfig/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/rs/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/apt/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/rs/androidTest/debug" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/resValues/androidTest/debug" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/res" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/resources" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/assets" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/aidl" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/java" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/rs" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/shaders" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/res" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/resources" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/assets" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/aidl" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/java" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/rs" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/shaders" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/main/res" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/main/resources" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/main/assets" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/main/aidl" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/main/java" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/main/rs" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/main/shaders" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/androidTest/res" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/androidTest/resources" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/androidTest/assets" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/androidTest/aidl" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/androidTest/java" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/androidTest/rs" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/androidTest/shaders" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/test/res" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/test/resources" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/test/assets" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/test/aidl" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/test/java" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/test/rs" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/test/shaders" isTestSource="true" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/annotations" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/blame" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/bundles" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/classes" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/incremental" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/incremental-safeguard" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/jniLibs" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/manifests" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/res" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/rs" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/shaders" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/symbols" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/transforms" />
<excludeFolder url="file://$MODULE_DIR$/build/outputs" />
<excludeFolder url="file://$MODULE_DIR$/build/tmp" />
</content>
<orderEntry type="jdk" jdkName="Android API 25 Platform" jdkType="Android SDK" />
<orderEntry type="sourceFolder" forTests="false" />
<orderEntry type="library" exported="" scope="TEST" name="rt" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="runner-0.5" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="espresso-idling-resource-2.2.2" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="hamcrest-library-1.3" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="hamcrest-integration-1.3" level="project" />
<orderEntry type="library" exported="" name="support-core-ui-25.3.0" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="jcommon-1.0.24" level="project" />
<orderEntry type="library" exported="" name="support-core-utils-25.3.0" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="jsr305-2.0.1" level="project" />
<orderEntry type="library" exported="" name="support-fragment-25.3.0" level="project" />
<orderEntry type="library" exported="" name="commons-math3-3.6.1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="espresso-core-2.2.2" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="exposed-instrumentation-api-publish-0.5" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="rules-0.5" level="project" />
<orderEntry type="library" exported="" name="iirj-v1.0.4" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="javax.annotation-api-1.2" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="jfreechart-1.0.19" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="javax.inject-1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="json-20160810" level="project" />
<orderEntry type="library" exported="" name="support-compat-25.3.0" level="project" />
<orderEntry type="library" exported="" name="support-v4-25.3.0" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="javawriter-2.1.1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="hamcrest-core-1.3" level="project" />
<orderEntry type="library" exported="" name="support-media-compat-25.3.0" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="junit-4.12" level="project" />
<orderEntry type="library" exported="" name="appcompat-v7-25.3.0" level="project" />
<orderEntry type="library" exported="" name="animated-vector-drawable-25.3.0" level="project" />
<orderEntry type="library" exported="" name="support-annotations-25.3.0" level="project" />
<orderEntry type="library" exported="" name="support-vector-drawable-25.3.0" level="project" />
</component>
</module>

77
signal/src/androidTest/java/net/heartshield/filter/FirPerfTest.java Normal file
View File

@@ -0,0 +1,77 @@
package net.heartshield.filter;
import junit.framework.TestCase;
import net.heartshield.filter.FFirFilterRunning;
import net.heartshield.filter.FirFilterRunning;
import net.heartshield.signal.Signal;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-15
*/
public class FirPerfTest extends TestCase {
public void testPerformanceDouble() {
double[] taps = Signal.triangle(200);
double[] sig = new double[20000];
double FPS = 200.0;
double F = 0.1;
for(int i = 0; i < sig.length; i++)
sig[i] = Math.sin(2*Math.PI*F/FPS*i);
FirFilterRunning fir = new FirFilterRunning(taps);
long t1 = System.nanoTime();
double[] y = fir.batch(sig);
long t2 = System.nanoTime();
throw new RuntimeException("time elapsed: " + (t2-t1) + " ns");
/*
* Nexus 5x: ~ 40 ms
* 28033752 ns
* 40981567 ns
* 46468442 ns
* 48293911 ns
*/
}
private static float[] d2f(double[] x) {
float[] f = new float[x.length];
for(int i = 0; i < x.length; i++)
f[i] = (float) x[i];
return f;
}
public void testPerformanceFloat() {
double[] taps = Signal.triangle(200);
float[] sig = new float[20000];
float FPS = 200.0f;
float F = 0.1f;
for(int i = 0; i < sig.length; i++)
sig[i] = (float) Math.sin(2*Math.PI*F/FPS*i);
FFirFilterRunning fir = new FFirFilterRunning(d2f(taps));
long t1 = System.nanoTime();
float[] y = fir.batch(sig);
long t2 = System.nanoTime();
throw new RuntimeException("time elapsed: " + (t2-t1) + " ns");
/*
* Nexus 5x:
* 31485888 ns
* 58385840 ns
* 54760839 ns
* 61453548 ns
*
* not faster - why?
*
* - CPU ALU is the same speed
* - debugging overhead prevails
*/
}
}

18
signal/src/androidTest/java/net/heartshield/signal/AudioDecoderTests.java Normal file
View File

@@ -0,0 +1,18 @@
package net.heartshield.signal;
import android.content.Context;
import android.support.test.InstrumentationRegistry;
import junit.framework.TestCase;
import java.io.IOException;
public class AudioDecoderTests extends TestCase {
public void testReadWav() throws IOException {
Context appContext = InstrumentationRegistry.getTargetContext();
AudioDecoder dec = new AudioDecoder(appContext);
short[] arr = dec.decodeToMemory(net.heartshield.signal.test.R.raw.test100hz);
if(arr.length == 0)
throw new IllegalStateException("failed");
}
}

26
signal/src/androidTest/java/net/heartshield/signal/ExampleInstrumentedTest.java Normal file
View File

@@ -0,0 +1,26 @@
package net.heartshield.signal;
import android.content.Context;
import android.support.test.InstrumentationRegistry;
import android.support.test.runner.AndroidJUnit4;
import org.junit.Test;
import org.junit.runner.RunWith;
import static org.junit.Assert.*;
/**
* Instrumentation test, which will execute on an Android device.
*
* @see <a href="http://d.android.com/tools/testing">Testing documentation</a>
*/
@RunWith(AndroidJUnit4.class)
public class ExampleInstrumentedTest {
@Test
public void useAppContext() throws Exception {
// Context of the app under test.
Context appContext = InstrumentationRegistry.getTargetContext();
assertEquals("net.heartshield.firdes.test", appContext.getPackageName());
}
}

BIN
signal/src/androidTest/res/raw/test100hz.wav Normal file
View File

Binary file not shown.

10
signal/src/main/AndroidManifest.xml Normal file
View File

@@ -0,0 +1,10 @@
<manifest xmlns:android="http://schemas.android.com/apk/res/android"
package="net.heartshield.signal">
<application android:allowBackup="true" android:label="@string/app_name"
android:supportsRtl="true">
</application>
</manifest>

11
signal/src/main/cpp/native-lib.cpp Normal file
View File

@@ -0,0 +1,11 @@
#include <jni.h>
#include <string>
extern "C"
JNIEXPORT jstring JNICALL
Java_net_heartshield_signal_FirFilters_stringFromJNI(
JNIEnv* env,
jobject /* this */) {
std::string hello = "Hello from C++ in firdes";
return env->NewStringUTF(hello.c_str());
}

57
signal/src/main/java/net/heartshield/filter/AlivecorFilter.java Normal file
View File

@@ -0,0 +1,57 @@
package net.heartshield.filter;
/**
* Realtime FM demodulator for AliveCor Kardia hand-held ECG.
* Introduces a delay and outputs some initial invalid samples (initial settling of the filter).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-29
*/
public class AlivecorFilter extends RFilterBlock {
private BandPassFilter mPrefilter;
private Hilbert mHilbert;
private PllFilter mPll;
private LowPassFilter mLowpass;
private BandRejectFilter mBandreject;
private int mDelay;
public AlivecorFilter(double fps) {
mPrefilter = new BandPassFilter(fps, 17000, 21000, 500);
mHilbert = new Hilbert();
mPll = new PllFilter(fps, 1500, 17000, 21000);
mLowpass = new LowPassFilter(fps, 200, 50);
//mBandreject = new BandRejectFilter(fps, 40, 60, 3);
mPrefilter.connect(mHilbert);
mHilbert.connect(mPll);
mPll.connect(mLowpass);
//mLowpass.connect(mBandreject);
mDelay = mPrefilter.getDelay() + mHilbert.getDelay() + mLowpass.getDelay(); // + mBandreject.getDelay();
}
public double getFps() {
return mPrefilter.getFps();
}
/** filter delay in number of samples */
public int getDelay() {
return mDelay;
}
@Override
protected double[] batch(double[] x) {
throw new IllegalStateException();
}
@Override
public void connect(IRFilterBlock consumer) {
//mBandreject.connect(consumer);
mLowpass.connect(consumer);
}
@Override
public void put(double[] x) {
mPrefilter.put(x);
}
}

16
signal/src/main/java/net/heartshield/filter/BandPassFilter.java Normal file
View File

@@ -0,0 +1,16 @@
package net.heartshield.filter;
import net.heartshield.signal.FirFilters;
/**
* Realtime band-pass filter.
* Introduces a delay and outputs some initial invalid samples (initial settling of the filter).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-29
*/
public class BandPassFilter extends FirFilter {
public BandPassFilter(double fps, double lowCutoffFreq, double highCutoffFreq, double transitionWidth) {
super(fps, FirFilters.bandPass(fps, lowCutoffFreq, highCutoffFreq, transitionWidth));
}
}

16
signal/src/main/java/net/heartshield/filter/BandRejectFilter.java Normal file
View File

@@ -0,0 +1,16 @@
package net.heartshield.filter;
import net.heartshield.signal.FirFilters;
/**
* Realtime band-reject filter.
* Introduces a delay and outputs some initial invalid samples (initial settling of the filter).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-29
*/
public class BandRejectFilter extends FirFilter {
public BandRejectFilter(double fps, double lowCutoffFreq, double highCutoffFreq, double transitionWidth) {
super(fps, FirFilters.bandReject(fps, lowCutoffFreq, highCutoffFreq, transitionWidth));
}
}

26
signal/src/main/java/net/heartshield/filter/ButterBandPassFilter.java Normal file
View File

@@ -0,0 +1,26 @@
package net.heartshield.filter;
import uk.me.berndporr.iirj.Butterworth;
/**
* Butterworth band-pass filter.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-15
*/
public class ButterBandPassFilter extends RFilterBlock {
private Butterworth mButterworth;
public ButterBandPassFilter(int order, double fps, double lowCutoffFreq, double highCutoffFreq) {
mButterworth = new Butterworth();
mButterworth.bandPass(order, fps, (lowCutoffFreq + highCutoffFreq) / 2.0, highCutoffFreq - lowCutoffFreq);
}
@Override
protected double[] batch(double[] x) {
double[] y = new double[x.length];
for(int i = 0; i < x.length; i++)
y[i] = mButterworth.filter(x[i]);
return y;
}
}

32
signal/src/main/java/net/heartshield/filter/CDataSink.java Normal file
View File

@@ -0,0 +1,32 @@
package net.heartshield.filter;
import net.heartshield.signal.CVec;
import org.apache.commons.math3.complex.Complex;
/**
* Filter block accepting data and storing it into an array.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public class CDataSink implements CFilterBlock {
private Complex[] mData;
CDataSink() {
mData = new Complex[0];
}
public Complex[] getData() {
return mData;
}
@Override
public void put(Complex[] x) {
mData = CVec.stack(mData, x);
}
protected Complex[] batch(Complex[] x) {
throw new IllegalStateException("should not call batch() here");
}
}

15
signal/src/main/java/net/heartshield/filter/CFilterBlock.java Normal file
View File

@@ -0,0 +1,15 @@
package net.heartshield.filter;
import org.apache.commons.math3.complex.Complex;
/**
* Realtime batch-processing filter block interface.
* Complex numbers.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public interface CFilterBlock {
/** batch-add an array and push results through consumers */
void put(Complex[] x);
}

27
signal/src/main/java/net/heartshield/filter/CRFilterBlock.java Normal file
View File

@@ -0,0 +1,27 @@
package net.heartshield.filter;
import org.apache.commons.math3.complex.Complex;
/**
* Realtime batch-processing filter block interface.
* Complex input, real output.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-29
*/
public abstract class CRFilterBlock implements CFilterBlock {
protected IRFilterBlock mConsumer;
public void connect(IRFilterBlock consumer) {
mConsumer = consumer;
}
/** batch-add an array and push results through consumers */
@Override
public void put(Complex[] x) {
mConsumer.put(this.batch(x));
}
/** batch-process an array and return array of output values */
protected abstract double[] batch(Complex[] x);
}

32
signal/src/main/java/net/heartshield/filter/DataSink.java Normal file
View File

@@ -0,0 +1,32 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
/**
* Filter block accepting data and storing it into an array.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public class DataSink extends SinkBlock {
private double[] mData;
public DataSink() {
mData = new double[0];
}
public double[] getData() {
return mData;
}
@Override
public void put(double[] x) {
mData = DVec.stack(mData, x);
}
@Override
protected double[] batch(double[] x) {
throw new IllegalStateException("should not call batch() here");
}
}

30
signal/src/main/java/net/heartshield/filter/Delay.java Normal file
View File

@@ -0,0 +1,30 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
/**
* Simple time delay filter block. Initialized with zeros.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public class Delay extends RFilterBlock {
private int mDelay;
private double[] mBuf;
/**
* @param delay delay in number of samples
*/
public Delay(int delay) {
mDelay = delay;
mBuf = DVec.zeros(delay);
}
@Override
protected double[] batch(double[] x) {
double[] buf = DVec.stack(mBuf, x);
double[] y = DVec.get(buf, 0, buf.length - mDelay);
mBuf = DVec.get(buf, buf.length - mDelay, buf.length);
return y;
}
}

22
signal/src/main/java/net/heartshield/filter/Diff.java Normal file
View File

@@ -0,0 +1,22 @@
package net.heartshield.filter;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-16
*/
public class Diff extends RFilterBlock {
double prev = 0;
// could also just have been FIR({-1, 1})
@Override
protected double[] batch(double[] x) {
double[] y = new double[x.length];
if(x.length == 0)
return y;
y[0] = x[0] - prev;
for(int i = 1; i < x.length; i++)
y[i] = x[i] - x[i-1];
prev = x[x.length - 1];
return y;
}
}

62
signal/src/main/java/net/heartshield/filter/Downsampler.java Normal file
View File

@@ -0,0 +1,62 @@
package net.heartshield.filter;
/**
* Downsamples by an integer ratio.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-30
*/
public class Downsampler extends RFilterBlock {
private int N; // left samples
private int M; // modulus
private double mFps;
public Downsampler(double fpsIn, int ratio) {
if(ratio <= 0)
throw new IllegalArgumentException("ratio must be >= 1");
M = ratio;
N = M;
mFps = fpsIn / ratio;
this.mDownsamplingRatio = ratio;
}
/** downsampled output FPS */
public double getFps() { return mFps; }
@Override
protected double[] batch(double[] x) {
/*
* N left samples 1 <= N <= M
* L input length 0 < L
* M modulus 1 <= M
*
* N || L
* +---+---+---++---+---+---+---+---+
* | | | || | | | | |
* +---+---+---++---+---+---+---+---+
*/
int L = x.length;
// too short to give even a single sample
if(N + L <= M) {
N += L;
return new double[0];
}
int start = M - N;
int outLen = (N + L - 1) / M;
double[] ret = new double[outLen];
int j = 0;
for(int i = start; i < L; i += M)
ret[j++] = x[i];
N = (N + L - 1) % M + 1;
return ret;
}
}

83
signal/src/main/java/net/heartshield/filter/EveningResampler.java Normal file
View File

@@ -0,0 +1,83 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
import java.util.ArrayList;
import java.util.List;
/**
* Evenly resamples signal values (t,x) that may be unevenly sampled.
* Linearly interpolates intermediate values.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-19
*/
public class EveningResampler {
public EveningResampler(double fps) {
mFps = fps;
}
/**
* Add uneven signal sample value `frame` at `time`.
* `time` must be strictly increasing and positive.
*/
public synchronized void addUneven(double time, double frame) {
if(time <= mPrevTime)
throw new IllegalArgumentException("time must be strictly increasing");
final int curFrame = time2frame(time), prevFrame = time2frame(mPrevTime);
final double curK = (frame - mPrevFrame) / (time - mPrevTime);
double k = mPrevK; // first entry must be added with previous slope
//System.out.println("addUneven(t=" + time + ", f=" + frame + ") prevFrame=" + prevFrame + " curFrame=" + curFrame);
//System.out.println("dt=" + (time - mPrevTime) + " dy=" + (frame - mPrevFrame));
for(int i = prevFrame; i < curFrame; i++) {
double t = (double) i / mFps;
// floor rounding to curFrame -> t < time
if(t >= time)
throw new IllegalStateException("t >= time. this should never happen");
double f = mPrevFrame + (t - mPrevTime) * k;
//System.out.println(" i=" + i + " " + "t=" + t + " delt=" + (t - mPrevTime) + " f=" + f);
mEvenTimes.add(t);
mEvenFrames.add(f);
add(t, f);
k = curK;
}
mPrevTime = time;
mPrevFrame = frame;
mPrevK = curK;
}
private int time2frame(double t) {
return (int) (mFps * t);
}
// note: these would make more sense in a combined call... since they may need to be consistent.
public synchronized double[] getEvenTimes() {
return DVec.toArray(mEvenTimes);
}
public synchronized double[] getEvenFrames() {
return DVec.toArray(mEvenFrames);
}
/** override this to receive even updates */
void add(double time, double frame) {}
public synchronized void clear() {
mPrevTime = -1e-3;
mPrevFrame = 0.0;
mPrevK = 0.0;
}
protected double mFps;
private double mPrevTime = -1e-3;
private double mPrevFrame = 0.0;
private double mPrevK = 0.0;
private List<Double> mEvenTimes = new ArrayList<>();
private List<Double> mEvenFrames = new ArrayList<>();
}

24
signal/src/main/java/net/heartshield/filter/FFilterBlock.java Normal file
View File

@@ -0,0 +1,24 @@
package net.heartshield.filter;
/**
* Realtime batch-processing filter block interface.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public abstract class FFilterBlock implements IFFilterBlock {
protected IFFilterBlock mConsumer;
public void connect(IFFilterBlock consumer) {
mConsumer = consumer;
}
/** batch-add an array and push results through consumers */
@Override
public void put(float[] x) {
mConsumer.put(this.batch(x));
}
/** batch-process an array and return array of output values */
protected abstract float[] batch(float[] x);
}

49
signal/src/main/java/net/heartshield/filter/FFirFilterRunning.java Normal file
View File

@@ -0,0 +1,49 @@
package net.heartshield.filter;
/**
* Realtime FIR filter using direct convolution.
*
* Introduces a delay that depends on the filter parameters (taps.length/2).
* Also, the first 2*delay output samples are invalid (initial settling of the filter).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-15
*/
public class FFirFilterRunning extends FFilterBlock {
private float[] taps;
private float[] xbuf;
public FFirFilterRunning(float[] taps) {
if(taps.length == 0)
throw new IllegalStateException("taps.length must be > 0");
this.taps = taps;
this.xbuf = new float[taps.length];
}
/** roll in-place by 1 position to left. in: [x_0 x_1 x_2] out: [x_1 x_2 x_0] */
private static void rollLeft(float[] x) {
if(x.length == 0)
return;
float x0 = x[0];
for(int i = 0; i < x.length-1; i++)
x[i] = x[i+1];
x[x.length-1] = x0;
}
private float convolveOne(float val) {
float sum = 0.0f;
rollLeft(xbuf);
xbuf[xbuf.length-1] = val;
for(int j = 0; j < taps.length; j++)
sum += xbuf[j] * taps[j];
return sum;
}
@Override
protected float[] batch(float[] x) {
float[] y = new float[x.length];
for(int i = 0; i < x.length; i++)
y[i] = convolveOne(x[i]);
return y;
}
}

118
signal/src/main/java/net/heartshield/filter/FirFilter.java Normal file
View File

@@ -0,0 +1,118 @@
package net.heartshield.filter;
import net.heartshield.signal.CVec;
import net.heartshield.signal.DVec;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.transform.DftNormalization;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.TransformType;
/**
* Realtime FIR filter using FFT convolution.
*
* Introduces a delay that depends on the filter parameters.
* Also, the first 2*delay output samples are invalid (initial settling of the filter).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public class FirFilter extends RFilterBlock {
private double[] mTaps;
private double mFps;
private int mDelay;
private double[] mBuf;
public FirFilter(double fps, double[] taps) {
mTaps = taps;
mFps = fps;
mDelay = taps.length / 2;
mBuf = DVec.zeros(taps.length-1);
}
/** filter delay in number of samples */
public int getDelay() {
return mDelay;
}
public double getFps() {
return mFps;
}
/*
// instead, use convolve() using FFT to speed things up
@Deprecated
@Override
protected double[] batch(double[] x) {
double[] buf = DVec.stack(mBuf, x);
int mathArraysLen = buf.length + mTaps.length - 1;
int trimLen = mathArraysLen - buf.length;
int trimLeft = trimLen / 2;
int trimRight = trimLen - trimLeft;
int trimEnd = mathArraysLen - trimRight;
// just keep trailing bit to include it into leading boundary of next batch
mBuf = DVec.get(buf, buf.length - mTaps.length+1, buf.length);
return DVec.get(MathArrays.convolve(buf, mTaps), trimLeft, trimEnd);
// to do: this has a slicing bug somewhere -- convolution using FFT yields correct results, while this convolution yields an offset between imag and real in Hilbert
// to do: (meaning there is some indexing trouble)
}
*/
@Override
protected double[] batch(double[] x) {
double[] buf = DVec.stack(mBuf, x);
// just keep trailing bit to include it into leading boundary of next batch
mBuf = DVec.get(buf, buf.length - mTaps.length+1, buf.length);
return convolveFft(buf, mTaps);
}
/**
* Applies a filter to a signal, implementing the overlap-save method (chunk up the signal)
* see https://en.wikipedia.org/wiki/Overlap%E2%80%93save_method
*
* This is equivalent to: np.convolve(x, taps, mode='valid') but more efficient for long signals
*/
private static double[] convolveFft(double[] sig, double[] taps) {
if(sig.length <= taps.length)
throw new IllegalArgumentException("signal must be longer than taps");
//double[] h = taps;
int M = taps.length;
int overlap = M - 1;
int N = nextpow2(4*overlap);
int stepSize = N - overlap;
//double[][] dataRI = new double[2][N];
//* <li>{@code dataRI[0][i]} is the real part of the {@code i}-th data point,</li>
//* <li>{@code dataRI[1][i]} is the imaginary part of the {@code i}-th data point.</li>
double[] h = new double[N];
System.arraycopy(taps, 0, h, 0, taps.length);
// rest is padded with zeros from the right
FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
Complex[] H = fft.transform(h, TransformType.FORWARD);
double[] x = DVec.stack(sig, DVec.zeros(N)); // end padding, so the last step_size batch will surely cover the end
double[] y = DVec.zeros(x.length);
for(int pos = 0; pos < x.length + 1 - N; pos += stepSize) {
Complex[] freqTransform = fft.transform(DVec.get(x, pos, pos+N), TransformType.FORWARD);
Complex[] freqFiltered = CVec.mul(freqTransform, H);
Complex[] yt = fft.transform(freqFiltered, TransformType.INVERSE);
for(int i = 0; i < stepSize; i++)
y[pos+i] = yt[M-1+i].getReal();
}
// cut back the end padding, and the overlap region where taps hang out of the signal
return DVec.get(y, 0, sig.length-taps.length+1);
}
private static int nextpow2(int n) {
double mf = Math.log(n) / Math.log(2);
double mi = Math.ceil(mf);
return (int) Math.pow(2, (int) mi);
}
}

49
signal/src/main/java/net/heartshield/filter/FirFilterRunning.java Normal file
View File

@@ -0,0 +1,49 @@
package net.heartshield.filter;
/**
* Realtime FIR filter using direct convolution.
*
* Introduces a delay that depends on the filter parameters (taps.length/2).
* Also, the first 2*delay output samples are invalid (initial settling of the filter).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-15
*/
public class FirFilterRunning extends RFilterBlock {
private double[] taps;
private double[] xbuf;
public FirFilterRunning(double[] taps) {
if(taps.length == 0)
throw new IllegalStateException("taps.length must be > 0");
this.taps = taps;
this.xbuf = new double[taps.length];
}
/** roll in-place by 1 position to left. in: [x_0 x_1 x_2] out: [x_1 x_2 x_0] */
private static void rollLeft(double[] x) {
if(x.length == 0)
return;
double x0 = x[0];
for(int i = 0; i < x.length-1; i++)
x[i] = x[i+1];
x[x.length-1] = x0;
}
private double convolveOne(double val) {
double sum = 0.0;
rollLeft(xbuf);
xbuf[xbuf.length-1] = val;
for(int j = 0; j < taps.length; j++)
sum += xbuf[j] * taps[j];
return sum;
}
@Override
protected double[] batch(double[] x) {
double[] y = new double[x.length];
for(int i = 0; i < x.length; i++)
y[i] = convolveOne(x[i]);
return y;
}
}

16
signal/src/main/java/net/heartshield/filter/HighPassFilter.java Normal file
View File

@@ -0,0 +1,16 @@
package net.heartshield.filter;
import net.heartshield.signal.FirFilters;
/**
* Realtime high-pass filter.
* Introduces a delay and outputs some initial invalid samples (initial settling of the filter).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public class HighPassFilter extends FirFilter {
public HighPassFilter(double fps, double cutoffFreq, double transitionWidth) {
super(fps, FirFilters.highPass(fps, cutoffFreq, transitionWidth));
}
}

40
signal/src/main/java/net/heartshield/filter/Hilbert.java Normal file
View File

@@ -0,0 +1,40 @@
package net.heartshield.filter;
import org.apache.commons.math3.complex.Complex;
/**
* Hilbert transform turns a series of real-valued samples into complex-valued (analytic) samples.
*
* Introduces a delay of 32 samples.
* Also, the first 32 output samples are invalid (initial settling of the filter).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public class Hilbert extends RCFilterBlock {
RFilterBlock mReal;
FirFilter mImag;
int mDelay;
public Hilbert() {
//super(FirFilters.hilbert(65), Double.NaN);
mImag = new HilbertImag();
mReal = new Delay(mImag.getDelay());
mDelay = mImag.getDelay();
}
/** filter delay in number of samples */
public int getDelay() {
return mDelay;
}
@Override
protected Complex[] batch(double[] x) {
double[] real = mReal.batch(x);
double[] imag = mImag.batch(x);
Complex[] result = new Complex[real.length];
for(int i = 0; i < real.length; i++)
result[i] = new Complex(real[i], imag[i]);
return result;
}
}

15
signal/src/main/java/net/heartshield/filter/HilbertImag.java Normal file
View File

@@ -0,0 +1,15 @@
package net.heartshield.filter;
import net.heartshield.signal.FirFilters;
/**
* Part of the Hilbert implementation, turns real part into imaginary part.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public class HilbertImag extends FirFilter {
public HilbertImag() {
super(Double.NaN, FirFilters.hilbert(65));
}
}

12
signal/src/main/java/net/heartshield/filter/IFFilterBlock.java Normal file
View File

@@ -0,0 +1,12 @@
package net.heartshield.filter;
/**
* Float-valued input, realtime batch-processing filter block interface.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public interface IFFilterBlock {
/** batch-add an array and push results through consumers */
void put(float[] x);
}

12
signal/src/main/java/net/heartshield/filter/IRFilterBlock.java Normal file
View File

@@ -0,0 +1,12 @@
package net.heartshield.filter;
/**
* Real-valued input, realtime batch-processing filter block interface.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public interface IRFilterBlock {
/** batch-add an array and push results through consumers */
void put(double[] x);
}

16
signal/src/main/java/net/heartshield/filter/LowPassFilter.java Normal file
View File

@@ -0,0 +1,16 @@
package net.heartshield.filter;
import net.heartshield.signal.FirFilters;
/**
* Realtime low-pass filter.
* Introduces a delay and outputs some initial invalid samples (initial settling of the filter).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public class LowPassFilter extends FirFilter {
public LowPassFilter(double fps, double cutoffFreq, double transitionWidth) {
super(fps, FirFilters.lowPass(fps, cutoffFreq, transitionWidth));
}
}

169
signal/src/main/java/net/heartshield/filter/PcgBeatDetector.java Normal file
View File

@@ -0,0 +1,169 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
import net.heartshield.signal.IVec;
import net.heartshield.signal.Signal;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-17
*/
public class PcgBeatDetector extends RFilterBlock {
private PcgEnvelopeFilter mPcgFilter;
private Splitter mSplitter = new Splitter();
private Splitter mSplitterPerc = new Splitter();
private Diff mDiff = new Diff(); // could also just have been FIR({-1, 1})
private FirFilterRunning mMean1;
//private FirFilterRunning mMean2;
private PercentileFilter mPercentileFilter;
private DiracUpsampler mUpsampler;
private RectUpsampler mRectUpsampler;
public PcgBeatDetector(double frontFps) {
mPcgFilter = new PcgEnvelopeFilter(frontFps);
double fps = mPcgFilter.getFps();
//mPcgFilter.connect(mSplitter);
//mPcgFilter.connect(mDiff);
double PERCENTILE_WIN_SECS = 4.0;
double PERCENTILE = 85;
// for plotting, 50.0 -- otherwise, use Double.MAX_VALUE
// why is input not normalized? (range until >30)
final double PERCENTILE_INIT = Double.MAX_VALUE; // ignore all beats until PercentileFilter has settled
mPercentileFilter = new PercentileFilter((int) (PERCENTILE_WIN_SECS * fps), PERCENTILE_INIT, PERCENTILE);
final int N1 = (int) (0.08 * fps);
//final int N2 = (int) (0.2334 * fps);
mMean1 = new FirFilterRunning(DVec.ones(N1));
//mMean2 = new FirFilterRunning(DVec.ones(N2));
//mDiff.connect(mMean1);
mPcgFilter.connect(mMean1);
//mMean1.connect(mMean2);
mMean1.connect(mSplitter);
//mPercentileFilter.connect(mPercentiles);
mPercentileFilter.connect(mSplitterPerc);
mSplitter.connect(mPercentileFilter); // must be before mBeatDetector, which uses mPercentiles.buf
mSplitter.connect(mBeatDetector);
mUpsampler = new DiracUpsampler(mPcgFilter.getDownsamplingRatio());
mRectUpsampler = new RectUpsampler(mPcgFilter.getDownsamplingRatio());
mBeatDetector.connect(mUpsampler);
mSplitterPerc.connect(mPercentiles);
mSplitterPerc.connect(mRectUpsampler);
// dummy consumer, unless someone connects
mRectUpsampler.connect(new IRFilterBlock() {
@Override
public void put(double[] x) {}
});
}
static class FilterBuf implements IRFilterBlock {
public double[] prevBuf = new double[1];
public double[] buf = new double[1];
@Override
public void put(double[] x) {
//System.out.println(" mPercentiles got x.length=" + x.length);
//if(x.length > 0) {
prevBuf = buf;
buf = x;
//}
}
}
private FilterBuf mPercentiles = new FilterBuf();
private RFilterBlock mBeatDetector = new RFilterBlock() {
private double[] prevSlice = new double[1];
// note: arrives in background thread
@Override
protected double[] batch(double[] x) {
//System.out.println(" mBeatDetector got x.length=" + x.length);
//double[] perc = mPercentileFilter.batch(x);
final double DIFF_THRESHOLD = 0.35;
double[] percentiles = DVec.stack(mPercentiles.prevBuf, mPercentiles.buf);
double[] threshold = DVec.mul(percentiles, DIFF_THRESHOLD);
double pval = mPercentiles.buf[0];
System.out.println("PcgBeatDetector: percentile = " + pval);
//Log.i(TAG, );
// would be nice if we could receive an edge-padded sig
// simulate by storing `prevSlice`
double[] sig = DVec.stack(prevSlice, x);
for(int i = 0; i < sig.length; i++)
if(sig[i] < threshold[i])
sig[i] = 0;
int[] sigBeats = Signal.heartbeat_localmax(sig);
// shave off the last entry - it may appear to be spurious max
int[] beats = IVec.get(sigBeats, sigBeats.length - x.length - 1, sigBeats.length - 1);
prevSlice = x;
return IVec.toDouble(beats);
}
};
static class DiracUpsampler extends RFilterBlock {
private int mRatio;
DiracUpsampler(int ratio) {
mRatio = ratio;
this.mUpsamplingRatio = ratio;
}
@Override
protected double[] batch(double[] x) {
double[] ret = new double[x.length * mRatio];
for(int i = 0; i < x.length; i++)
ret[i*mRatio] = x[i];
return ret;
}
}
static class RectUpsampler extends RFilterBlock {
private int mRatio;
RectUpsampler(int ratio) {
mRatio = ratio;
this.mUpsamplingRatio = ratio;
}
@Override
protected double[] batch(double[] x) {
double[] ret = new double[x.length * mRatio];
for(int i = 0; i < x.length; i++)
for(int j = i * mRatio; j < (i+1) * mRatio; j++)
ret[j] = x[i];
return ret;
}
}
@Override
protected double[] batch(double[] x) {
throw new IllegalStateException();
}
@Override
public void connect(IRFilterBlock consumer) {
mUpsampler.connect(consumer);
}
public void connectPercentiles(IRFilterBlock consumer) {
mRectUpsampler.connect(consumer);
}
@Override
public void put(double[] x) {
mPcgFilter.put(x);
}
}

60
signal/src/main/java/net/heartshield/filter/PcgEnvelopeFilter.java Normal file
View File

@@ -0,0 +1,60 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
import net.heartshield.signal.Signal;
/**
* Rough envelope detection for heart sounds (phonocardiography),
* to provide user feedback.
*
* Downsamples.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-15
*/
public class PcgEnvelopeFilter extends RFilterBlock {
private FirFilterRunning mMovingAverage;
private Downsampler mDownsampler;
private ButterBandPassFilter mBandPass;
private RectifierFilter mRectifier;
private FirFilterRunning mLowPass;
private double mBackFps;
private int mDownsamplingRatio = 24; // 48 kHz -> 2 kHz
public PcgEnvelopeFilter(double frontFps) {
// note: performance wise, we could do better than these two.
mMovingAverage = new FirFilterRunning(DVec.ones(mDownsamplingRatio)); // TODO: fps-dependent
mDownsampler = new Downsampler(frontFps, mDownsamplingRatio);
mBackFps = mDownsampler.getFps();
mBandPass = new ButterBandPassFilter(3, mBackFps, 30.0, 250.0);
mRectifier = new RectifierFilter();
mLowPass = new FirFilterRunning(Signal.triangle(200));
mMovingAverage.connect(mDownsampler);
mDownsampler.connect(mBandPass);
mBandPass.connect(mRectifier);
mRectifier.connect(mLowPass);
}
public int getDownsamplingRatio() { return mDownsamplingRatio; }
public double getFps() {
return mBackFps;
}
@Override
protected double[] batch(double[] x) {
throw new IllegalStateException();
}
@Override
public void connect(IRFilterBlock consumer) {
mLowPass.connect(consumer);
}
@Override
public void put(double[] x) {
mMovingAverage.put(x);
}
}

35
signal/src/main/java/net/heartshield/filter/PercentileFilter.java Normal file
View File

@@ -0,0 +1,35 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
import net.heartshield.signal.Signal;
/**
* Running percentile in a recent window.
*
* Includes a single-window delay and the first window's output samples are 0.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-16
*/
public class PercentileFilter extends WindowedBatchFilter {
private double q;
/**
* @param win window size in samples
* @param q percentile in range [0:100]
*/
public PercentileFilter(int win, double initial, double q) {
super(win, initial);
this.q = q;
}
@Override
protected double processWindow(double[] buf, int start) {
return Signal.percentile(DVec.get(buf, start, start + mDelay), q);
}
@Override
public double[] batch(double[] x) {
return super.batch(x);
}
}

82
signal/src/main/java/net/heartshield/filter/Pll.java Normal file
View File

@@ -0,0 +1,82 @@
package net.heartshield.filter;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.util.FastMath;
/**
* Low-level implementation of phase-locked loop.
* Works with normalized frequency internally.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
class Pll {
private double mPhase = 0;
private double mFreq = 0;
private double mAlpha, mBeta;
private double mMaxFreq, mMinFreq;
private static double TWO_PI = 2*Math.PI;
/**
* TODO convert freqs
*
* internal settings are in terms of radians per sample, not Hz.
* see https://en.wikipedia.org/wiki/Normalized_frequency_%28unit%29
*
* @param bandwidth
* @param minFreq
* @param maxFreq
*/
Pll(double bandwidth, double minFreq, double maxFreq) {
if(bandwidth < 0.0)
throw new IllegalArgumentException("bandwidth must be >= 0");
if(maxFreq < 0.0)
throw new IllegalArgumentException("maxFreq must be >= 0");
if(minFreq < 0.0)
throw new IllegalArgumentException("minFreq must be >= 0");
if(maxFreq < minFreq)
throw new IllegalArgumentException("maxFreq must be >= minFreq");
mMaxFreq = maxFreq;
mMinFreq = minFreq;
double damping = Math.sqrt(2.0) / 2.0;
double denom = (1.0 + 2.0 * damping * bandwidth + bandwidth * bandwidth);
mAlpha = (4.0 * damping * bandwidth) / denom;
mBeta = (4.0 * bandwidth * bandwidth) / denom;
}
void batch(Complex[] input, double[] output) {
double err;
for(int i = 0; i < output.length; i++) {
output[i] = mFreq;
Complex in = input[i];
double phase = FastMath.atan2(in.getImaginary(), in.getReal());
err = mod2pi(phase - mPhase);
mFreq += mBeta * err;
mPhase += mFreq + mAlpha * err;
while(mPhase > TWO_PI)
mPhase -= TWO_PI;
while(mPhase < -TWO_PI)
mPhase += TWO_PI;
if(mFreq > mMaxFreq)
mFreq = mMaxFreq;
else if(mFreq < mMinFreq)
mFreq = mMinFreq;
}
}
private static double mod2pi(double in) {
if(in > Math.PI)
return in - TWO_PI;
else if(in < -Math.PI)
return in + TWO_PI;
else
return in;
}
}

34
signal/src/main/java/net/heartshield/filter/PllFilter.java Normal file
View File

@@ -0,0 +1,34 @@
package net.heartshield.filter;
import org.apache.commons.math3.complex.Complex;
/**
* Phase-locked loop filter.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-29
*/
public class PllFilter extends CRFilterBlock {
private double mFps;
private Pll mPll;
/**
* @param bandwidth determines the lock range (Hz)
* @param minFreq minimum frequency (Hz)
* @param maxFreq maximum frequency (Hz)
*/
public PllFilter(double fps, double bandwidth, double minFreq, double maxFreq) {
mFps = fps;
mPll = new Pll(normFreq(bandwidth), normFreq(minFreq), normFreq(maxFreq));
}
/** @return normalized frequency */
private double normFreq(double f) { return f * 2.0 * Math.PI / mFps; }
@Override
protected double[] batch(Complex[] x) {
double[] ret = new double[x.length];
mPll.batch(x, ret);
return ret;
}
}

182
signal/src/main/java/net/heartshield/filter/PlotFilter.java Normal file
View File

@@ -0,0 +1,182 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
import net.heartshield.signal.FirFilters;
import net.heartshield.signal.Signal;
/**
* Processes uneven addition of samples into a curve that is plotted.
* introduces a delay of TAPS_LEN_RIGHT.
*
* Results are not precise. Rather, this filter is optimized for less delay, as intended for visuals.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-19
*/
public class PlotFilter extends EveningResampler {
private static final double[] TAPS; // FIR highpass fps=30 fc=0.5 tw=0.5
private static final int TAPS_LEN_RIGHT; // length from [center : last non-zero)
private static final double[] TAPS_SHORT; // non-centered active taps
private static final int PLOT_SMOOTHING_N;
private int lpad;
private int ltaps;
private double[] x;
private double[] y;
private double[] t;
private int iw;
private int ir;
private int iy;
private boolean initialized;
private int ypad_left;
private int measure_len;
public PlotFilter(double fps, int measure_len, int ypad_left) {
super(fps);
lpad = TAPS_SHORT.length;
ltaps = TAPS_SHORT.length;
x = DVec.stack(DVec.ones(lpad), DVec.zeros(measure_len));
y = DVec.zeros(ypad_left + TAPS_LEN_RIGHT + measure_len);
t = DVec.zeros(ypad_left + TAPS_LEN_RIGHT + measure_len);
iw = lpad-1;
ir = 0;
iy = ypad_left;
initialized = false;
this.ypad_left = ypad_left;
this.measure_len = measure_len;
}
public void clear() {
super.clear();
x = DVec.stack(DVec.ones(lpad), DVec.zeros(measure_len));
y = DVec.zeros(ypad_left + TAPS_LEN_RIGHT + measure_len);
t = DVec.zeros(ypad_left + TAPS_LEN_RIGHT + measure_len);
iw = lpad-1;
ir = 0;
iy = ypad_left;
initialized = false;
}
@Override
void add(double time, double frame) {
if(!initialized) {
initialized = true;
// pad to left with initial value
for(int i = 0; i < lpad; i++)
//x.set(i, frame);
x[i] = frame;
// retrofit some timestamps
for(int i = 0; i < ypad_left; i++)
t[i] = time + (i - ypad_left) / this.mFps;
}
// this should not happen! alas, if FPS is off, or user waits for too long without starting a measurement, it does.
if(iw >= x.length)
x = DVec.stack(x, DVec.zeros(x.length));
if(iy >= y.length) {
y = DVec.stack(y, DVec.zeros(y.length));
t = DVec.stack(t, DVec.zeros(t.length));
}
x[iw++] = frame;
y[iy] = DVec.sum(DVec.mul(DVec.get(x, ir, ir+ltaps), TAPS_SHORT));
t[iy] = time;
ir++;
iy++;
}
public int size() { return ir; }
double get() {
return y[iy-1];
}
public double[] getPlotPointsUnscaled() {
double[] buf = DVec.get(y, iy-ypad_left, iy);
buf = Signal.running_mean(Signal.running_mean(buf, PLOT_SMOOTHING_N), PLOT_SMOOTHING_N);
final int min_len = ypad_left-2*(PLOT_SMOOTHING_N-1);
if(buf.length < min_len) {
// left-pad with zeros if necessary, to come up with a plot
buf = DVec.stack(DVec.zeros(min_len), buf);
}
/*buf = DVec.mul(buf, 0.5 / DVec.max(DVec.abs(buf)));
buf = DVec.add(buf, 0.5);*/
return buf;
}
private double[] getPlotPoints() {
double[] buf = getRecentPoints(ypad_left);
buf = DVec.mul(buf, -1.0);
buf = DVec.mul(buf, 0.5 / DVec.max(DVec.abs(buf)));
buf = DVec.add(buf, 0.5);
return buf;
}
private double[] getPlotTimes() {
// to do: still off by a few samples? see PlotFilterTests
return getRecentTimes(ypad_left);
}
public synchronized Series getPlotSeries() {
return new Series(getPlotTimes(), getPlotPoints());
}
// needs synchronized access!
private double[] getRecentPoints(int length) {
double[] buf = DVec.get(y, iy-length, iy);
buf = Signal.running_mean(Signal.running_mean(buf, PLOT_SMOOTHING_N), PLOT_SMOOTHING_N);
final int min_len = length-2*(PLOT_SMOOTHING_N-1);
if(buf.length < min_len) {
// left-pad with zeros if necessary, to come up with a plot
buf = DVec.stack(DVec.zeros(min_len), buf);
}
return buf;
}
// needs synchronized access!
private double[] getRecentTimes(int length) {
// to do: still off by a few samples? see PlotFilterTests
return DVec.get(t, iy-length +1*(PLOT_SMOOTHING_N-1), iy - 1*(PLOT_SMOOTHING_N-1));
}
public static class Series {
public double[] t;
public double[] x;
public Series(double[] t, double[] x) {
this.t = t;
this.x = x;
}
}
public synchronized Series getRecentSeries(int length) {
return new Series(getRecentTimes(length), getRecentPoints(length));
}
static {
// TODO(david): fix FPS dependency (TAPS and TAPS_LEN_RIGHT are designed for 30 fps)
//TAPS = new double[]{ -2.54277053e-04, -2.76789215e-04, -2.99490261e-04, -3.22486390e-04, -3.45740002e-04, -3.69053741e-04, -3.92059737e-04, -4.14213369e-04, -4.34791815e-04, -4.52898588e-04, -4.67473350e-04, -4.77307389e-04, -4.81064926e-04, -4.77309397e-04, -4.64535347e-04, -4.41204436e-04, -4.05785773e-04, -3.56799457e-04, -2.92862998e-04, -2.12739193e-04, -1.15384959e-04, 2.91419969e-19, 1.33925787e-04, 2.86567025e-04, 4.57721006e-04, 6.46770000e-04, 8.52648926e-04, 1.07381924e-03, 1.30824978e-03, 1.55340519e-03, 1.80624390e-03, 2.06322316e-03, 2.32031452e-03, 2.57303007e-03, 2.81645381e-03, 3.04528838e-03, 3.25390534e-03, 3.43640824e-03, 3.58670158e-03, 3.69856670e-03, 3.76574416e-03, 3.78202391e-03, 3.74133443e-03, 3.63783888e-03, 3.46602965e-03, 3.22082476e-03, 2.89766071e-03, 2.49258382e-03, 2.00233469e-03, 1.42442901e-03, 7.57227419e-04, -9.38187527e-19, -8.47022515e-04, -1.78260717e-03, -2.80449004e-03, -3.90936062e-03, -5.09286067e-03, -6.34959433e-03, -7.67315552e-03, -9.05616954e-03, -1.04903458e-02, -1.19665479e-02, -1.34748695e-02, -1.50047382e-02, -1.65450070e-02, -1.80840734e-02, -1.96100064e-02, -2.11106613e-02, -2.25738306e-02, -2.39873696e-02, -2.53393278e-02, -2.66181119e-02, -2.78125945e-02, -2.89122574e-02, -2.99073178e-02, -3.07888426e-02, -3.15488540e-02, -3.21804322e-02, -3.26777995e-02, -3.30363698e-02, -3.32528502e-02, 9.66431737e-01, -3.32528502e-02, -3.30363698e-02, -3.26777995e-02, -3.21804322e-02, -3.15488540e-02, -3.07888426e-02, -2.99073178e-02, -2.89122574e-02, -2.78125945e-02, -2.66181119e-02, -2.53393278e-02, -2.39873696e-02, -2.25738306e-02, -2.11106613e-02, -1.96100064e-02, -1.80840734e-02, -1.65450070e-02, -1.50047382e-02, -1.34748695e-02, -1.19665479e-02, -1.04903458e-02, -9.05616954e-03, -7.67315552e-03, -6.34959433e-03, -5.09286067e-03, -3.90936062e-03, -2.80449004e-03, -1.78260717e-03, -8.47022515e-04, -9.38187527e-19, 7.57227419e-04, 1.42442901e-03, 2.00233469e-03, 2.49258382e-03, 2.89766071e-03, 3.22082476e-03, 3.46602965e-03, 3.63783888e-03, 3.74133443e-03, 3.78202391e-03, 3.76574416e-03, 3.69856670e-03, 3.58670158e-03, 3.43640824e-03, 3.25390534e-03, 3.04528838e-03, 2.81645381e-03, 2.57303007e-03, 2.32031452e-03, 2.06322316e-03, 1.80624390e-03, 1.55340519e-03, 1.30824978e-03, 1.07381924e-03, 8.52648926e-04, 6.46770000e-04, 4.57721006e-04, 2.86567025e-04, 1.33925787e-04, 2.91419969e-19, -1.15384959e-04, -2.12739193e-04, -2.92862998e-04, -3.56799457e-04, -4.05785773e-04, -4.41204436e-04, -4.64535347e-04, -4.77309397e-04, -4.81064926e-04, -4.77307389e-04, -4.67473350e-04, -4.52898588e-04, -4.34791815e-04, -4.14213369e-04, -3.92059737e-04, -3.69053741e-04, -3.45740002e-04, -3.22486390e-04, -2.99490261e-04, -2.76789215e-04, -2.54277053e-04 };
double fps=30.0;
double cf=0.5;
double tw=0.5;
TAPS = FirFilters.highPass(fps, cf, tw);
double tapsLenRightSecs = 0.5;
TAPS_LEN_RIGHT = (int) (tapsLenRightSecs * fps); // 15
double plotSmoothingWinSecs = 0.2;
PLOT_SMOOTHING_N = (int) (plotSmoothingWinSecs * fps); // 6
double[] sl = DVec.get(TAPS, 0, TAPS.length/2 + TAPS_LEN_RIGHT);
TAPS_SHORT = DVec.sub(sl, DVec.sum(sl) / (double) sl.length);
}
}

27
signal/src/main/java/net/heartshield/filter/RCFilterBlock.java Normal file
View File

@@ -0,0 +1,27 @@
package net.heartshield.filter;
import org.apache.commons.math3.complex.Complex;
/**
* Realtime batch-processing filter block interface.
* Real input, complex output.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public abstract class RCFilterBlock implements IRFilterBlock {
protected CFilterBlock mConsumer;
public void connect(CFilterBlock consumer) {
mConsumer = consumer;
}
/** batch-add an array and push results through consumers */
@Override
public void put(double[] x) {
mConsumer.put(this.batch(x));
}
/** batch-process an array and return array of output values */
protected abstract Complex[] batch(double[] x);
}

44
signal/src/main/java/net/heartshield/filter/RFilterBlock.java Normal file
View File

@@ -0,0 +1,44 @@
package net.heartshield.filter;
/**
* Realtime batch-processing filter block interface.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public abstract class RFilterBlock implements IRFilterBlock {
protected IRFilterBlock mConsumer;
protected int mDownsamplingRatio = 1;
protected int mUpsamplingRatio = 1;
public void connect(IRFilterBlock consumer) {
mConsumer = consumer;
}
/** batch-add an array and push results through consumers */
@Override
public void put(double[] x) {
double[] res = this.batch(x);
throwIfBadResultLength(x, res);
if(mConsumer == null)
throw new IllegalStateException("must connect a consumer to " + this.getClass().getName());
mConsumer.put(res);
}
private void throwIfBadResultLength(double[] x, double[] res) {
//if(x.length * mUpsamplingRatio != res.length * mDownsamplingRatio)
// maybe it's one more sample that's been output
double s1 = x.length * mUpsamplingRatio;
double s2 = res.length * mDownsamplingRatio;
if(mUpsamplingRatio == 1 && mDownsamplingRatio == 1 && res.length != x.length)
throw new IllegalStateException("x.length=" + x.length + " != res.length=" + res.length + " from batch(x) in " + this.getClass().getName());
if(Math.abs(s1 - s2) / Math.max(mUpsamplingRatio, mDownsamplingRatio) > 1.0)
throw new IllegalStateException("x.length=" + x.length + " should yield " + ((double) (x.length * mUpsamplingRatio)) / mDownsamplingRatio + " samples, but got " + res.length);
}
/** batch-process an array and return array of output values */
protected abstract double[] batch(double[] x);
}

14
signal/src/main/java/net/heartshield/filter/RectifierFilter.java Normal file
View File

@@ -0,0 +1,14 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-15
*/
public class RectifierFilter extends RFilterBlock {
@Override
protected double[] batch(double[] x) {
return DVec.abs(x);
}
}

19
signal/src/main/java/net/heartshield/filter/SinkBlock.java Normal file
View File

@@ -0,0 +1,19 @@
package net.heartshield.filter;
/**
* Filter block accepting data.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public abstract class SinkBlock extends RFilterBlock {
@Override
protected double[] batch(double[] x) {
throw new IllegalStateException("should not call batch() here");
}
@Override
public void put(double[] x) {
throw new IllegalStateException("override me");
}
}

14
signal/src/main/java/net/heartshield/filter/SourceBlock.java Normal file
View File

@@ -0,0 +1,14 @@
package net.heartshield.filter;
/**
* Filter block providing data, e.g. from some input hardware device.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-26
*/
public class SourceBlock extends RFilterBlock {
@Override
protected double[] batch(double[] x) {
return x;
}
}

23
signal/src/main/java/net/heartshield/filter/Splitter.java Normal file
View File

@@ -0,0 +1,23 @@
package net.heartshield.filter;
import java.util.ArrayList;
import java.util.List;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-16
*/
public class Splitter implements IRFilterBlock {
protected List<IRFilterBlock> mConsumers = new ArrayList<>();
public void connect(IRFilterBlock consumer) {
mConsumers.add(consumer);
}
/** batch-add an array and push results through consumers */
@Override
public void put(double[] x) {
for(IRFilterBlock consumer : mConsumers)
consumer.put(x);
}
}

97
signal/src/main/java/net/heartshield/filter/WindowedBatchFilter.java Normal file
View File

@@ -0,0 +1,97 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
/**
* Window-based almost-realtime processing of a signal.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-16
*/
public class WindowedBatchFilter extends RFilterBlock {
protected int mDelay;
protected double[] mBuf;
protected int N; // left samples
protected int M; // modulus
protected double mPrevVal = 0.0;
/**
* @param win window size in samples
* @param initial initialization value, pretended to be to the left of the first input batch
*/
public WindowedBatchFilter(int win, double initial) {
if(win <= 0)
throw new IllegalArgumentException("win must be >= 1");
mDelay = win;
mBuf = DVec.mul(DVec.ones(win), initial);
M = win;
N = M;
}
protected double processWindow(double[] buf, int start) {
// simple downsampler
return buf[start + mDelay];
}
@Override
protected double[] batch(double[] x) {
/*
* N left samples 1 <= N <= M
* L input length 0 < L
* M modulus 1 <= M
*
* N || L
* +---+---+---++---+---+---+---+---+
* | | | || | | | | |
* +---+---+---++---+---+---+---+---+
*/
int L = x.length;
double[] buf = DVec.stack(mBuf, x);
// too short to process even a single window
if(N + L <= M) {
N += L;
mBuf = DVec.get(buf, buf.length - mDelay, buf.length);
return DVec.mul(DVec.ones(L), mPrevVal);
}
int start = M - N;
//int outLen = (N + L - 1) / M;
//double[] ret = new double[outLen];
double[] ret = new double[x.length];
//double val = processWindow(buf, start);
for(int j = 0; j < start; j++)
ret[j] = mPrevVal;
// processWindow() result of the whole window BEFORE is used for the next window.
// i.e. [x_0 x_1 x_2 x_3] with M=2 becomes
// [0.0 0.0 y01 y01]
//
// add [x_4] -> [y23]
// add [x_5] -> [y23]
//
// where yij = processWindow([x_i x_j])
double val = mPrevVal;
for(int i = start; i < L; i += M) {
val = processWindow(buf, i);
int end = Math.min(i + M, L);
for(int k = i; k < end; k++)
ret[k] = val;
}
mPrevVal = val;
N = (N + L - 1) % M + 1;
mBuf = DVec.get(buf, buf.length - mDelay, buf.length);
return ret;
}
}

75
signal/src/main/java/net/heartshield/filter/WindowedDownsampler.java Normal file
View File

@@ -0,0 +1,75 @@
package net.heartshield.filter;
import net.heartshield.signal.DVec;
/**
* Window-based almost-realtime processing of a signal.
* Outputs one sample per window (is downsampling).
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-16
*/
public class WindowedDownsampler extends RFilterBlock {
protected int mDelay;
protected double[] mBuf;
protected int N; // left samples
protected int M; // modulus
/**
* @param win window size in samples
*/
public WindowedDownsampler(int win) {
mDelay = win;
mBuf = new double[win];
if(win <= 0)
throw new IllegalArgumentException("ratio must be >= 1");
M = win;
N = M;
}
protected double processWindow(double[] buf, int start) {
// simple downsampler
return buf[start + mDelay];
}
@Override
protected double[] batch(double[] x) {
/*
* N left samples 1 <= N <= M
* L input length 0 < L
* M modulus 1 <= M
*
* N || L
* +---+---+---++---+---+---+---+---+
* | | | || | | | | |
* +---+---+---++---+---+---+---+---+
*/
int L = x.length;
double[] buf = DVec.stack(mBuf, x);
// too short to give even a single sample
if(N + L <= M) {
N += L;
mBuf = DVec.get(buf, buf.length - mDelay, buf.length);
return new double[0];
}
int start = M - N;
int outLen = (N + L - 1) / M;
double[] ret = new double[outLen];
int j = 0;
for(int i = start; i < L; i += M)
ret[j++] = processWindow(buf, i);
N = (N + L - 1) % M + 1;
mBuf = DVec.get(buf, buf.length - mDelay, buf.length);
return ret;
}
}

139
signal/src/main/java/net/heartshield/signal/AudioDecoder.java Normal file
View File

@@ -0,0 +1,139 @@
package net.heartshield.signal;
import android.content.Context;
import android.content.res.AssetFileDescriptor;
import android.content.res.Resources;
import android.media.MediaCodec;
import android.media.MediaExtractor;
import android.media.MediaFormat;
import android.util.Log;
import java.io.IOException;
import java.nio.ByteBuffer;
import java.util.Arrays;
import static junit.framework.Assert.assertEquals;
import static junit.framework.Assert.assertTrue;
public class AudioDecoder {
private static final String TAG = "AudioDecoder";
private Resources mResources;
public AudioDecoder(Context context) {
mResources = context.getResources();
}
// https://android.googlesource.com/platform/cts/+/jb-mr2-release/tests/tests/media/src/android/media/cts/DecoderTest.java
// Apache License, Version 2.0
public short[] decodeToMemory(int resId/*, boolean reconfigure*/) throws IOException {
boolean reconfigure = false;
short [] decoded = new short[0];
int decodedIdx = 0;
AssetFileDescriptor testFd = mResources.openRawResourceFd(resId);
MediaExtractor extractor;
MediaCodec codec;
ByteBuffer[] codecInputBuffers;
ByteBuffer[] codecOutputBuffers;
extractor = new MediaExtractor();
extractor.setDataSource(testFd.getFileDescriptor(), testFd.getStartOffset(),
testFd.getLength());
testFd.close();
assertEquals("wrong number of tracks", 1, extractor.getTrackCount());
MediaFormat format = extractor.getTrackFormat(0);
String mime = format.getString(MediaFormat.KEY_MIME);
assertTrue("not an audio file", mime.startsWith("audio/"));
codec = MediaCodec.createDecoderByType(mime);
codec.configure(format, null /* surface */, null /* crypto */, 0 /* flags */);
codec.start();
codecInputBuffers = codec.getInputBuffers();
codecOutputBuffers = codec.getOutputBuffers();
if (reconfigure) {
codec.stop();
codec.configure(format, null, null, 0);
codec.start();
codecInputBuffers = codec.getInputBuffers();
codecOutputBuffers = codec.getOutputBuffers();
}
extractor.selectTrack(0);
// start decoding
final long kTimeOutUs = 5000;
MediaCodec.BufferInfo info = new MediaCodec.BufferInfo();
boolean sawInputEOS = false;
boolean sawOutputEOS = false;
int noOutputCounter = 0;
while (!sawOutputEOS && noOutputCounter < 50) {
noOutputCounter++;
if (!sawInputEOS) {
int inputBufIndex = codec.dequeueInputBuffer(kTimeOutUs);
if (inputBufIndex >= 0) {
ByteBuffer dstBuf = codecInputBuffers[inputBufIndex];
int sampleSize =
extractor.readSampleData(dstBuf, 0 /* offset */);
long presentationTimeUs = 0;
if (sampleSize < 0) {
Log.d(TAG, "saw input EOS.");
sawInputEOS = true;
sampleSize = 0;
} else {
presentationTimeUs = extractor.getSampleTime();
}
codec.queueInputBuffer(
inputBufIndex,
0 /* offset */,
sampleSize,
presentationTimeUs,
sawInputEOS ? MediaCodec.BUFFER_FLAG_END_OF_STREAM : 0);
if (!sawInputEOS) {
extractor.advance();
}
}
}
int res = codec.dequeueOutputBuffer(info, kTimeOutUs);
if (res >= 0) {
//Log.d(TAG, "got frame, size " + info.size + "/" + info.presentationTimeUs);
if (info.size > 0) {
noOutputCounter = 0;
}
if (info.size > 0 && reconfigure) {
// once we've gotten some data out of the decoder, reconfigure it again
reconfigure = false;
extractor.seekTo(0, MediaExtractor.SEEK_TO_NEXT_SYNC);
sawInputEOS = false;
codec.stop();
codec.configure(format, null /* surface */, null /* crypto */, 0 /* flags */);
codec.start();
codecInputBuffers = codec.getInputBuffers();
codecOutputBuffers = codec.getOutputBuffers();
continue;
}
int outputBufIndex = res;
ByteBuffer buf = codecOutputBuffers[outputBufIndex];
if (decodedIdx + (info.size / 2) >= decoded.length) {
decoded = Arrays.copyOf(decoded, decodedIdx + (info.size / 2));
}
for (int i = 0; i < info.size; i += 2) {
decoded[decodedIdx++] = buf.getShort(i);
}
codec.releaseOutputBuffer(outputBufIndex, false /* render */);
if ((info.flags & MediaCodec.BUFFER_FLAG_END_OF_STREAM) != 0) {
Log.d(TAG, "saw output EOS.");
sawOutputEOS = true;
}
} else if (res == MediaCodec.INFO_OUTPUT_BUFFERS_CHANGED) {
codecOutputBuffers = codec.getOutputBuffers();
Log.d(TAG, "output buffers have changed.");
} else if (res == MediaCodec.INFO_OUTPUT_FORMAT_CHANGED) {
MediaFormat oformat = codec.getOutputFormat();
Log.d(TAG, "output format has changed to " + oformat);
} else {
Log.d(TAG, "dequeueOutputBuffer returned " + res);
}
}
codec.stop();
codec.release();
return decoded;
}
}

42
signal/src/main/java/net/heartshield/signal/CVec.java Normal file
View File

@@ -0,0 +1,42 @@
package net.heartshield.signal;
import org.apache.commons.math3.complex.Complex;
/**
* Numerical utilities on Complex[] -- compare numpy functions.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-28
*/
public class CVec {
/** concatenate two or more arrays */
public static Complex[] stack(Complex[]... arrs) {
int len = 0;
for(Complex[] x : arrs)
len += x.length;
Complex[] a = new Complex[len];
int i = 0;
for(Complex[] x : arrs) {
System.arraycopy(x, 0, a, i, x.length);
i += x.length;
}
return a;
}
private static void throwIfLengthMismatch(Complex[] x, Complex[] y) {
if(x.length != y.length)
throw new IllegalArgumentException("x.length=" + x.length + " != y.length=" + y.length);
}
/** element-wise multiplication */
public static Complex[] mul(Complex[] x, Complex[] y) {
throwIfLengthMismatch(x, y);
Complex[] a = new Complex[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i].multiply(y[i]);
return a;
}
}

167
signal/src/main/java/net/heartshield/signal/CameraUtil.java Normal file
View File

@@ -0,0 +1,167 @@
package net.heartshield.signal;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-04-09
*/
public class CameraUtil {
static public double fractionRed(int[] rgb, int width, int height) {
long nred = 0;
for(int j = 0, yp = 0; j < height; j++) {
for(int i = 0; i < width; i++, yp++) {
final int r = (rgb[yp] & 0xff0000) >> 16;
final int g = (rgb[yp] & 0xff00) >> 8;
final int b = (rgb[yp] & 0xff);
if(r > g && r > b)
nred++;
}
}
return ((double) nred) / (width * height);
}
/**
* Compute column-wise mean.
*
* @param rgb input image pixel matrix
* @param width input width
* @param height input height
* @return output row (length = width)
*/
static public double[][] meanRow(int[] rgb, int width, int height) {
double[][] rgbMean = new double[3][width];
for(int j = 0, yp = 0; j < height; j++) {
for(int i = 0; i < width; i++, yp++) {
rgbMean[0][i] += (rgb[yp] & 0xff0000) >> 16;
rgbMean[1][i] += (rgb[yp] & 0xff00) >> 8;
rgbMean[2][i] += (rgb[yp] & 0xff);
}
}
for(int i = 0; i < width; i++) {
for(int k = 0; k < 3; k++)
rgbMean[k][i] /= height;
}
return rgbMean;
}
/**
* Compute row-wise mean.
*
* @param rgb input image pixel matrix
* @param width input width
* @param height input height
* @return output column (length = height)
*/
static public double[][] meanCol(int[] rgb, int width, int height) {
double[][] rgbMean = new double[3][height];
for(int j = 0, yp = 0; j < height; j++) {
for(int i = 0; i < width; i++, yp++) {
rgbMean[0][j] += (rgb[yp] & 0xff0000) >> 16;
rgbMean[1][j] += (rgb[yp] & 0xff00) >> 8;
rgbMean[2][j] += (rgb[yp] & 0xff);
}
}
for(int i = 0; i < height; i++) {
for(int k = 0; k < 3; k++)
rgbMean[k][i] /= width;
}
return rgbMean;
}
/**
* Compute column-wise standard deviation.
*
* @param rgb input image pixel matrix
* @param width input width
* @param height input height
* @return output row (length = width)
*/
static public double[][] stdRow(int[] rgb, double[][] rgbMean, int width, int height) {
double[][] rgbStd = new double[3][width];
for(int j = 0, yp = 0; j < height; j++) {
for(int i = 0; i < width; i++, yp++) {
double r = ((rgb[yp] & 0xff0000) >> 16) - rgbMean[0][i];
double g = ((rgb[yp] & 0xff00) >> 8) - rgbMean[1][i];
double b = (rgb[yp] & 0xff) - rgbMean[2][i];
rgbStd[0][i] += r * r;
rgbStd[1][i] += g * g;
rgbStd[2][i] += b * b;
}
}
for(int i = 0; i < width; i++) {
for(int k = 0; k < 3; k++)
rgbStd[k][i] = Math.sqrt(rgbStd[k][i] / height);
}
return rgbStd;
}
/**
* Compute row-wise standard deviation.
*
* @param rgb input image pixel matrix
* @param width input width
* @param height input height
* @return output row (length = height)
*/
static public double[][] stdCol(int[] rgb, double[][] rgbMean, int width, int height) {
double[][] rgbStd = new double[3][height];
for(int j = 0, yp = 0; j < height; j++) {
for(int i = 0; i < width; i++, yp++) {
double r = ((rgb[yp] & 0xff0000) >> 16) - rgbMean[0][j];
double g = ((rgb[yp] & 0xff00) >> 8) - rgbMean[1][j];
double b = (rgb[yp] & 0xff) - rgbMean[2][j];
rgbStd[0][j] += r * r;
rgbStd[1][j] += g * g;
rgbStd[2][j] += b * b;
}
}
for(int i = 0; i < height; i++) {
for(int k = 0; k < 3; k++)
rgbStd[k][i] = Math.sqrt(rgbStd[k][i] / width);
}
return rgbStd;
}
static public void decodeYUV420SP(int[] rgb, double[] mean_rgb, byte[] yuv420sp, int width, int height) {
final int frameSize = width * height;
long rr = 0, gg = 0, bb = 0;
for (int j = 0, yp = 0; j < height; j++) {
int uvp = frameSize + (j >> 1) * width, u = 0, v = 0;
for (int i = 0; i < width; i++, yp++) {
int y = (0xff & ((int) yuv420sp[yp])) - 16;
if (y < 0) y = 0;
if ((i & 1) == 0) {
u = (0xff & yuv420sp[uvp++]) - 128;
v = (0xff & yuv420sp[uvp++]) - 128;
}
int y1192 = 1192 * y;
int r = (y1192 + 1634 * v);
int g = (y1192 - 833 * v - 400 * u);
int b = (y1192 + 2066 * u);
if (r < 0) r = 0; else if (r > 262143) r = 262143;
if (g < 0) g = 0; else if (g > 262143) g = 262143;
if (b < 0) b = 0; else if (b > 262143) b = 262143;
rgb[yp] = 0xff000000 | ((r << 6) & 0xff0000) | ((g >> 2) & 0xff00) | ((b >> 10) & 0xff);
rr += (r >> 10) & 0xff;
gg += (g >> 10) & 0xff;
bb += (b >> 10) & 0xff;
}
}
mean_rgb[0] = ((double) rr) / (width * height);
mean_rgb[1] = ((double) gg) / (width * height);
mean_rgb[2] = ((double) bb) / (width * height);
}
}

388
signal/src/main/java/net/heartshield/signal/DVec.java Normal file
View File

@@ -0,0 +1,388 @@
package net.heartshield.signal;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
/**
* Numerical utilities on double[] -- compare numpy functions.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-19
*/
public class DVec {
/** return zero-array of specified length */
public static double[] zeros(int len) {
return new double[len];
}
/** return one-array of specified length */
public static double[] ones(int len) {
double[] a = new double[len];
for(int i = 0; i < len; i++)
a[i] = 1.0;
return a;
}
public static double[] arange(double start, double stop) {
return arange(start, stop, 1.0);
}
/** evenly spaced values within the interval [start, stop) */
public static double[] arange(double start, double stop, double step) {
int len = (int) ((stop - start) / step);
if(len <= 0)
return new double[0];
double[] a = new double[len];
int i = 0;
for(double d = start; i < len; d += step) {
a[i] = d;
i++;
}
return a;
}
public static double[] linspace(double start, double stop, int num) {
return linspace(start, stop, num, true);
}
/** `num` evenly spaced numbers over interval [start,stop] with endpoint=true or [start,stop) with endpoint=false */
public static double[] linspace(double start, double stop, int num, boolean endpoint) {
if(num < 0)
throw new IllegalArgumentException("num must be >= 0");
int end = endpoint ? num : (num-1);
double step = (stop - start) / (double) end;
double[] a = new double[num];
double d = start;
for(int i = 0; i < num; i++) {
a[i] = d;
d += step;
}
return a;
}
/** concatenate two or more arrays */
public static double[] stack(double[]... arrs) {
int len = 0;
for(double[] x : arrs)
len += x.length;
double[] a = new double[len];
int i = 0;
for(double[] x : arrs) {
System.arraycopy(x, 0, a, i, x.length);
i += x.length;
}
return a;
}
/** returns the sorted version of x */
public static double[] sort(double[] x) {
double[] y = Arrays.copyOf(x, x.length);
Arrays.sort(y);
return y;
}
private static void throwIfOutOfBounds(double[] x, int start, int end) {
int len = end - start;
if(x.length == 0) {
if(start != 0 || end != 0)
throw new IndexOutOfBoundsException("out of bounds on empty Vec");
} else {
// getting a zero-length slice at the end is still OK
if(start > x.length || start < 0)
throw new IndexOutOfBoundsException("start out of bounds");
if(end > x.length || end < 0)
throw new IndexOutOfBoundsException("end out of bounds");
}
if(len < 0)
throw new IndexOutOfBoundsException("end >= start was violated");
}
private static void throwIfLengthMismatch(double[] x, double[] y) {
if(x.length != y.length)
throw new IllegalArgumentException("x.length=" + x.length + " != y.length=" + y.length);
}
/** get a slice of the array[start,end) (deep copy) */
public static double[] get(double[] x, int start, int end) {
int len = end - start;
throwIfOutOfBounds(x, start, end);
double[] a = new double[len];
System.arraycopy(x, start, a, 0, len);
return a;
}
/** element-wise addition */
public static double[] add(double[] x, double[] y) {
throwIfLengthMismatch(x, y);
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] + y[i];
return a;
}
/** element-wise addition of a scalar */
public static double[] add(double[] x, double y) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] + y;
return a;
}
/** element-wise subtraction of a scalar */
public static double[] sub(double[] x, double y) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] - y;
return a;
}
/** element-wise subtraction */
public static double[] sub(double[] x, double[] y) {
throwIfLengthMismatch(x, y);
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] - y[i];
return a;
}
/** element-wise multiplication with a scalar */
public static double[] mul(double[] x, double y) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] * y;
return a;
}
/** element-wise multiplication */
public static double[] mul(double[] x, double[] y) {
throwIfLengthMismatch(x, y);
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] * y[i];
return a;
}
/** element-wise comparison, without tolerance -> boolean */
public static int[] eq(double[] x, double[] y) {
throwIfLengthMismatch(x, y);
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (x[i] == y[i]) ? 1 : 0;
return a;
}
/** element-wise comparison, without tolerance -> boolean */
public static int[] lt(double[] x, double[] y) {
throwIfLengthMismatch(x, y);
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (x[i] < y[i]) ? 1 : 0;
return a;
}
/** element-wise comparison, with tolerance -> boolean */
public static int[] eqt(double[] x, double[] y, double tol) {
throwIfLengthMismatch(x, y);
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (Math.abs(x[i] - y[i]) < tol) ? 1 : 0;
return a;
}
/** element-wise comparison to a scalar, without tolerance -> boolean */
public static int[] eq(double[] x, double y) {
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (x[i] == y) ? 1 : 0;
return a;
}
/**
* element-wise signum function;
* * zero if the argument is zero,
* * 1.0 if the argument is greater than zero,
* * -1.0 if the argument is less than zero.
*/
public static double[] sign(double[] x) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = Math.signum(x[i]);
return a;
}
/** element-wise rounding */
public static double[] round(double[] x) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = Math.round(x[i]);
return a;
}
/** element-wise abs */
public static double[] abs(double[] x) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = Math.abs(x[i]);
return a;
}
/** element-wise sin */
public static double[] sin(double[] x) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = Math.sin(x[i]);
return a;
}
/** element-wise exp */
public static double[] exp(double[] x) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = Math.exp(x[i]);
return a;
}
/** numerical differences */
public static double[] diff(double[] x) {
if(x.length == 0)
throw new IllegalArgumentException("cannot diff() empty array");
double[] a = new double[x.length - 1];
for(int i = 0; i < x.length - 1; i++)
a[i] = x[i+1] - x[i];
return a;
}
/** cumulative sum */
public static double[] cumsum(double[] x) {
double[] a = new double[x.length];
if(x.length > 0)
a[0] = x[0];
for(int i = 1; i < x.length; i++)
a[i] = a[i-1] + x[i];
return a;
}
/** sum of all elements */
public static double sum(double[] x) {
double s = 0.0;
for(int i = 0; i < x.length; i++)
s += x[i];
return s;
}
/** max of all elements */
public static double max(double[] x) {
if(x.length == 0)
throw new IllegalArgumentException("max is undefined for zero elements");
double m = Double.MIN_VALUE;
for(int i = 0; i < x.length; i++) {
if(x[i] > m)
m = x[i];
}
return m;
}
/** min of all elements */
public static double min(double[] x) {
if(x.length == 0)
throw new IllegalArgumentException("max is undefined for zero elements");
double m = Double.MAX_VALUE;
for(int i = 0; i < x.length; i++) {
if(x[i] < m)
m = x[i];
}
return m;
}
/** mean of elements */
public static double mean(double[] x) {
if(x.length == 0)
throw new IllegalArgumentException("mean is undefined for zero elements");
return sum(x) / x.length;
}
/** standard deviation of elements */
public static double std(double[] x) {
if(x.length == 0)
throw new IllegalArgumentException("std is undefined for zero elements");
double m = mean(x);
double s = 0.0;
for(int i = 0; i < x.length; i++)
s += (x[i] - m) * (x[i] - m);
return Math.sqrt(s / x.length);
}
/** indexing x by idx */
public static double[] index(double[] x, int[] idx) {
double[] a = new double[idx.length];
for(int i = 0; i < idx.length; i++)
a[i] = x[idx[i]];
return a;
}
/** integer downsampling */
public static double[] downsample(double[] x, int ratio) {
double[] a = new double[x.length / ratio];
if(x.length > 0)
a[0] = x[0];
for(int i = 1; i < a.length; i++)
a[i] = x[i*ratio];
return a;
}
public static String toString(double[] x) {
StringBuilder sb = new StringBuilder();
sb.append("{");
for(int i = 0; i < x.length; i++) {
if(i != 0)
sb.append(", ");
sb.append(Double.toString(x[i]));
}
sb.append("}");
return sb.toString();
}
/** convert to int[] */
public static int[] toInteger(double[] x) {
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (int) x[i];
return a;
}
/** convert to short[] */
public static short[] toShort(double[] x) {
short[] a = new short[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (short) x[i];
return a;
}
/** convert List<Double> into double[] */
public static double[] toArray(List<Double> lst) {
double[] a = new double[lst.size()];
int l = lst.size();
for(int i = 0; i < l; i++)
a[i] = lst.get(i);
return a;
}
/** convert double[] into List<Double> */
public static List<Double> toList(double[] x) {
List<Double> lst = new ArrayList<Double>(x.length);
for(int i = 0; i < x.length; i++)
lst.add(x[i]);
return lst;
}
}

154
signal/src/main/java/net/heartshield/signal/FirFilters.java Normal file
View File

@@ -0,0 +1,154 @@
package net.heartshield.signal;
public class FirFilters {
/**
* A native method that is implemented by the 'native-lib' native library,
* which is packaged with this application.
*/
//public native String stringFromJNI();
// Used to load the 'native-lib' library on application startup.
/*static {
System.loadLibrary("native-lib");
}*/
/** rect-window Hilbert transform impulse response */
public static double[] hilbert(int len) {
if(len % 2 == 0)
throw new IllegalArgumentException("len must be odd");
int hwin = (len-1) / 2;
double[] arr = new double[len];
double gain = 0.0;
for(int i = 1; i <= hwin; i++) {
if(i % 2 == 1) {
arr[i + hwin] = 1.0 / i;
arr[hwin - i] = -1.0 / i;
gain = arr[i + hwin] - gain;
} else {
arr[i + hwin] = 0;
arr[hwin - i] = 0;
}
}
gain = 2 * Math.abs(gain);
for(int i = 0; i < len; i++)
arr[i] /= gain;
return arr;
}
public static double[] hamming(int len) {
double[] arr = new double[len];
double M = len-1;
for(int i = 0; i < len; i++)
arr[i] = 0.54 - 0.46 * Math.cos(2 * Math.PI * i / M);
return arr;
}
/** low-pass filter impulse response design */
public static double[] lowPass(double fps, double cutoffFreq, double transitionWidth) {
int len = filterLength(fps, transitionWidth);
double[] arr = new double[len];
double[] win = hamming(len);
int M = (len-1) / 2;
double tpkm = 2 * Math.PI * cutoffFreq / fps;
for(int i = 0; i < len; i++) {
if(i == M)
arr[i] = tpkm / Math.PI;
else
arr[i] = Math.sin(tpkm * (i-M)) / (Math.PI * (i-M));
arr[i] *= win[i];
}
double denom = DVec.sum(arr);
for(int i = 0; i < len; i++)
arr[i] /= denom;
return arr;
}
/** high-pass filter impulse response design */
public static double[] highPass(double fps, double cutoffFreq, double transitionWidth) {
int len = filterLength(fps, transitionWidth);
double[] arr = new double[len];
double[] win = hamming(len);
int M = (len-1) / 2;
double tpkm = 2 * Math.PI * cutoffFreq / fps;
for(int i = 0; i < len; i++) {
if(i == M)
arr[i] = 1 - tpkm / Math.PI;
else
arr[i] = -Math.sin(tpkm * (i-M)) / (Math.PI * (i-M));
arr[i] *= win[i];
}
double denom = 0;
for(int i = 0; i < len; i++)
denom += arr[i] * Math.cos(i * Math.PI);
for(int i = 0; i < len; i++)
arr[i] /= denom;
return arr;
}
/** band-pass filter impulse response design */
public static double[] bandPass(double fps, double lowCutoffFreq, double highCutoffFreq, double transitionWidth) {
int len = filterLength(fps, transitionWidth);
double[] arr = new double[len];
double[] win = hamming(len);
int M = (len-1) / 2;
double tpkmLow = 2 * Math.PI * lowCutoffFreq / fps;
double tpkmHigh = 2 * Math.PI * highCutoffFreq / fps;
for(int i = 0; i < len; i++) {
if(i == M)
arr[i] = (tpkmHigh - tpkmLow) / Math.PI;
else
arr[i] = (Math.sin(tpkmHigh * (i-M)) - Math.sin(tpkmLow * (i-M))) / (Math.PI * (i-M));
arr[i] *= win[i];
}
double denom = 0;
for(int i = 0; i < len; i++)
denom += arr[i] * Math.cos((i-M) * (tpkmLow + tpkmHigh) / 2.0);
for(int i = 0; i < len; i++)
arr[i] /= denom;
return arr;
}
/** band-reject filter impulse response design */
public static double[] bandReject(double fps, double lowCutoffFreq, double highCutoffFreq, double transitionWidth) {
int len = filterLength(fps, transitionWidth);
double[] arr = new double[len];
double[] win = hamming(len);
int M = (len-1) / 2;
double tpkmLow = 2 * Math.PI * lowCutoffFreq / fps;
double tpkmHigh = 2 * Math.PI * highCutoffFreq / fps;
for(int i = 0; i < len; i++) {
if(i == M)
arr[i] = 1.0 - (tpkmHigh - tpkmLow) / Math.PI;
else
arr[i] = -(Math.sin(tpkmHigh * (i-M)) - Math.sin(tpkmLow * (i-M))) / (Math.PI * (i-M));
arr[i] *= win[i];
}
double denom = DVec.sum(arr);
for(int i = 0; i < len; i++)
arr[i] /= denom;
return arr;
}
private static void throwIfInvalidFreqs(double fps, double cutoffFreq, double transitionWidth) {
if(fps < 0)
throw new IllegalArgumentException("fps must be >= 0");
if(cutoffFreq < 0 || cutoffFreq > fps / 2.0)
throw new IllegalArgumentException("0 < cutoffFreq < fps/2");
if(transitionWidth < 0)
throw new IllegalArgumentException("transitionWidth must be >= 0");
}
private static int filterLength(double fps, double transitionWidth) {
int len = (int) (3.0 * fps / transitionWidth);
len += 1 - (len & 1);
return len;
}
}

249
signal/src/main/java/net/heartshield/signal/IVec.java Normal file
View File

@@ -0,0 +1,249 @@
package net.heartshield.signal;
import java.util.ArrayList;
import java.util.List;
/**
* Numerical utilities on int[] -- compare numpy functions.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-19
*/
public class IVec {
/** return zero-array of specified length */
public static int[] zeros(int len) {
return new int[len];
}
private static void throwIfLengthMismatch(int[] x, int[] y) {
if(x.length != y.length)
throw new IllegalArgumentException("x.length=" + x.length + " != y.length=" + y.length);
}
public static int[] arange(int start, int stop) {
return arange(start, stop, 1);
}
/** evenly spaced values within the interval [start, stop) */
public static int[] arange(int start, int stop, int step) {
int len = ((stop - start) / step);
if(len <= 0)
return new int[0];
int[] a = new int[len];
int i = 0;
for(int d = start; d < stop; d += step) {
a[i] = d;
i++;
}
return a;
}
/** concatenate two or more arrays */
public static int[] stack(int[]... arrs) {
int len = 0;
for(int[] x : arrs)
len += x.length;
int[] a = new int[len];
int i = 0;
for(int[] x : arrs) {
System.arraycopy(x, 0, a, i, x.length);
i += x.length;
}
return a;
}
/** element-wise addition */
public static int[] add(int[] x, int[] y) {
throwIfLengthMismatch(x, y);
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] + y[i];
return a;
}
/** element-wise addition of a scalar */
public static int[] add(int[] x, int y) {
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] + y;
return a;
}
/** element-wise subtraction */
public static int[] sub(int[] x, int[] y) {
throwIfLengthMismatch(x, y);
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] - y[i];
return a;
}
/** element-wise multiplication */
public static int[] mul(int[] x, int[] y) {
throwIfLengthMismatch(x, y);
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i] * y[i];
return a;
}
/** element-wise comparison to a scalar -> boolean */
public static int[] gt(int[] x, int y) {
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (x[i] > y) ? 1 : 0;
return a;
}
/** element-wise comparison to a scalar -> boolean */
public static int[] lt(int[] x, int y) {
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (x[i] < y) ? 1 : 0;
return a;
}
/** element-wise comparison to a scalar -> boolean */
public static int[] ge(int[] x, int y) {
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (x[i] >= y) ? 1 : 0;
return a;
}
/** element-wise comparison -> boolean */
public static int[] eq(int[] x, int[] y) {
throwIfLengthMismatch(x, y);
int[] a = new int[x.length];
for(int i = 0; i < x.length; i++)
a[i] = (x[i] == y[i]) ? 1 : 0;
return a;
}
/** numerical differences */
public static int[] diff(int[] x) {
if(x.length == 0)
throw new IllegalArgumentException("cannot diff() empty array");
int[] a = new int[x.length - 1];
for(int i = 0; i < x.length - 1; i++)
a[i] = x[i+1] - x[i];
return a;
}
/** indexing x by idx */
public static int[] index(int[] x, int[] idx) {
int[] a = new int[idx.length];
for(int i = 0; i < idx.length; i++)
a[i] = x[idx[i]];
return a;
}
/** array that has ones at `idx` locations and zeros otherwise. */
public static int[] dirac(int len, int[] idx) {
int[] a = new int[len];
for(int i = 0; i < idx.length; i++)
a[idx[i]] = 1;
return a;
}
/** true if array is purely ones, i.e. x[i] == 1 for all i */
public static boolean all(int[] x) {
for(int i = 0; i < x.length; i++)
if(x[i] != 1)
return false;
return true;
}
/** sum of all elements */
public static int sum(int[] x) {
int s = 0;
for(int i = 0; i < x.length; i++)
s += x[i];
return s;
}
public static String toString(int[] x) {
StringBuilder sb = new StringBuilder();
sb.append("{");
for(int i = 0; i < x.length; i++) {
if(i != 0)
sb.append(", ");
sb.append(Integer.toString(x[i]));
}
sb.append("}");
return sb.toString();
}
/** convert to double[] */
public static double[] toDouble(int[] x) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i];
return a;
}
/** convert List<Integer> into int[] */
public static int[] toArray(List<Integer> lst) {
int[] a = new int[lst.size()];
int l = lst.size();
for(int i = 0; i < l; i++)
a[i] = lst.get(i);
return a;
}
/** convert int[] into List<Integer> */
public static List<Integer> toList(int[] x) {
List<Integer> lst = new ArrayList<Integer>(x.length);
for(int i = 0; i < x.length; i++)
lst.add(x[i]);
return lst;
}
/** array of indices i where x[i] == 1 */
public static int[] where(int[] x) {
int nnz = 0;
int l = x.length;
for(int i = 0; i < l; i++)
if(x[i] == 1) nnz++;
int[] a = new int[nnz];
int j = 0;
for(int i = 0; i < l; i++)
if(x[i] == 1)
a[j++] = i;
return a;
}
private static void throwIfOutOfBounds(int[] x, int start, int end) {
int len = end - start;
if(x.length == 0) {
if(start != 0 || end != 0)
throw new IndexOutOfBoundsException("out of bounds on empty Vec");
} else {
if(start >= x.length || start < 0)
throw new IndexOutOfBoundsException("start out of bounds");
if(end > x.length || end < 0)
throw new IndexOutOfBoundsException("end out of bounds");
}
if(len < 0)
throw new IndexOutOfBoundsException("end >= start was violated");
}
/** get a slice of the array[start,end) (deep copy) */
public static int[] get(int[] x, int start, int end) {
int len = end - start;
throwIfOutOfBounds(x, start, end);
int[] a = new int[len];
System.arraycopy(x, start, a, 0, len);
return a;
}
}

242
signal/src/main/java/net/heartshield/signal/RunningQuality.java Normal file
View File

@@ -0,0 +1,242 @@
package net.heartshield.signal;
import android.util.Log;
import android.util.Pair;
import java.util.ArrayList;
import java.util.List;
/**
* Running signal quality indicator.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-04-23
*/
public class RunningQuality {
private static final String TAG = "RunningQuality";
/** template beat is resampled to this #samples */
private static final int BEAT_LEN = 30;
/** threshold for accepting initial beats */
private static final double BEAT_CORR_THR_1 = 0.8;
/** threshold for accepting subsequent beats */
public static final double BEAT_CORR_THR_2 = 0.6;
/** for second measurement after failed first one */
public static final double BEAT_CORR_THR_2_LOWER = 0.45;
private double mBeatCorrThr2 = BEAT_CORR_THR_2;
private Double startTime = null;
private List<Double> x = new ArrayList<>();
private List<Integer> edges = new ArrayList<>();
private double fps;
private List<double[]> beatTemplates = new ArrayList<>();
private double[] beatTemplate = null;
private List<IRange> badBeatRanges = new ArrayList<>();
private List<Pair<Integer, Double>> labels = new ArrayList<>();
private List<Double> runningCorrs = new ArrayList<>();
private Integer startIdx = null;
private int numNoisy = 0;
//private int numTotal = 0;
/** Pythonic range [start:end) */
public static class IRange {
public int start;
public int end;
public IRange(int start, int end) {
this.start = start;
this.end = end;
}
}
public RunningQuality(double fps) {
this.fps = fps;
}
public void clear() {
x = new ArrayList<>();
edges = new ArrayList<>();
beatTemplates = new ArrayList<>();
beatTemplate = null;
badBeatRanges = new ArrayList<>();
labels = new ArrayList<>();
runningCorrs = new ArrayList<>();
startIdx = null;
numNoisy = 0;
//numTotal = 0;
startTime = null;
Log.i(TAG, "clear()");
}
/*
* edge: detected downslope
* beat: timespan between two edges
*/
/**
* @param t time of sample 0 in this slice (but only used for very first call)
* @param sl recent slice to append
*/
public void append(double t, double[] sl) {
//Log.i(TAG, "append(sl.length=" + sl.length + ")");
if(startTime == null)
startTime = t;
final int beatdet_hwin = 30;
int prev_ix = Math.max(x.size() - beatdet_hwin, 0);
x.addAll(DVec.toList(sl));
// nb. CPU waste. nice-to: incremental beatdetect()
int[] idx = IVec.where(Signal.beatdetect(fps, DVec.toArray(x)));
Log.i(TAG, "append(sl.length=" + sl.length + "). all idx: " + IVec.toString(idx));
// ignore new beats unless closer within - at least 1 sec (half-window of savgol in getrr())
int[] recent = IVec.ge(idx, prev_ix);
int[] oldEnough = IVec.lt(idx, x.size() - beatdet_hwin);
int[] new_idx = IVec.index(idx, IVec.where(IVec.mul(recent, oldEnough)));
int prev_iedge = edges.size();
edges.addAll(IVec.toList(new_idx));
if(edges.size() > 1 && new_idx.length > 0)
appendBeats(prev_ix, prev_iedge);
}
public interface BeatListener {
/**
* Listener callback for first-time locking.
*
* @param startIdx index of sample where this beat starts
*/
void onLocked(int startIdx);
/**
* Listener callback for all detected beats, good and bad.
* Note that this may be dispatched several times in succession.
*
* @param startIdx index of sample where this beat starts
* @param goodBeat whether this beat is good
* @param posCorrCoeff positive-clipped correlation coefficient in range [0:1]
*/
void onBeat(int startIdx, boolean goodBeat, double posCorrCoeff);
}
private BeatListener mListener;
public void setListener(BeatListener listener) { mListener = listener; }
private void dispatchBeat(int startIdx, boolean goodBeat, double posCorrCoeff) {
if(mListener != null)
mListener.onBeat(startIdx, goodBeat, posCorrCoeff);
}
private void dispatchLocked(int startIdx) {
if(mListener != null)
mListener.onLocked(startIdx);
}
public double getTimeAt(int idx) {
return startTime + idx / fps;
}
public List<IRange> getBadBeatRanges() { return badBeatRanges; }
private void appendBeats(int prev_ix, int prev_iedge) {
int start = Math.max(prev_iedge-1, 0);
int end = edges.size()-1;
for(int i = start; i < end; i++) {
int s = edges.get(i);
int e = edges.get(i+1);
// TODO: should ignore crazy-long and very short beats here. (filter up on beat detector)
//Log.i(TAG, "appendBeat(" + s + ", " + e + ")");
double[] beat = resample(DVec.toArray(x.subList(s, e)));
double corr = Double.NaN;
double posCorr = Double.NaN;
boolean goodBeat = false;
if(beatTemplates.size() > 0) {
corr = Signal.crossCorr(beat, beatTemplate);
posCorr = Math.min(Math.max(corr, 0.0), 1.0);
double corrThreshold;
if(beatTemplates.size() > 2)
corrThreshold = mBeatCorrThr2;
else
corrThreshold = BEAT_CORR_THR_1;
goodBeat = (corr > corrThreshold);
}
boolean justLocked = false;
if(beatTemplates.size() == 0) {
// cannot correlate the first beat, no template yet
Log.i(TAG, "appendBeat(" + s + ", " + e + ") first beat");
addTemplate(beat);
} else if(beatTemplates.size() <= 2) {
Log.i(TAG, "appendBeat(" + s + ", " + e + ") init corr=" + corr);
// initialization: restart if there is no clear correlation between beats
// nb.: this may fail for low BP with double-peaks of current getrr()
if(goodBeat) {
addTemplate(beat);
if(beatTemplates.size() > 2)
justLocked = true;
} else {
replaceTemplate(beat);
badBeatRanges.clear();
}
} else {
Log.i(TAG, "appendBeat(" + s + ", " + e + ") run corr=" + corr);
// running mode: collect bad beats, but may be OK not to restart immediately
if(startIdx == null)
startIdx = s;
if(goodBeat) {
// keep updating beat template (nb. this gets O(n**2) currently)
addTemplate(beat);
} else {
badBeatRanges.add(new IRange(s, e));
numNoisy++;
}
runningCorrs.add(posCorr);
if(justLocked)
dispatchLocked(s);
dispatchBeat(s, goodBeat, posCorr);
}
if(!Double.isNaN(corr))
labels.add(new Pair<>(s, corr));
}
}
public boolean isRunning() {
return beatTemplates.size() > 2;
}
public void setBeatCorrThr2(double beatCorrThr2) {
mBeatCorrThr2 = beatCorrThr2;
}
private void addTemplate(double[] beat) {
beatTemplates.add(beat);
// nb. this gets O(n**2) currently
beatTemplate = Signal.mean(beatTemplates);
}
private void replaceTemplate(double[] beat) {
beatTemplates.clear();
beatTemplates.add(beat);
beatTemplate = Signal.mean(beatTemplates);
}
/** resample to BEAT_LEN */
private double[] resample(double[] x) {
double[] t = DVec.linspace(0, x.length, BEAT_LEN, false);
return Signal.interp(t, IVec.toDouble(IVec.arange(0, x.length)), x);
}
}

17
signal/src/main/java/net/heartshield/signal/SVec.java Normal file
View File

@@ -0,0 +1,17 @@
package net.heartshield.signal;
/**
* Numerical utilities on short[] -- compare numpy functions.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-29
*/
public class SVec {
/** convert to double[] */
public static double[] toDouble(short[] x) {
double[] a = new double[x.length];
for(int i = 0; i < x.length; i++)
a[i] = x[i];
return a;
}
}

69
signal/src/main/java/net/heartshield/signal/Series.java Normal file
View File

@@ -0,0 +1,69 @@
package net.heartshield.signal;
import org.json.JSONArray;
import org.json.JSONException;
import org.json.JSONObject;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStreamReader;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-04-21
*/
public class Series {
public double[] t;
public double[] x;
public Series(double[] t, double[] x) {
this.t = t;
this.x = x;
}
public Series(String fileName) {
try {
load(fileName);
} catch (JSONException e) {
throw new RuntimeException(e);
} catch (IOException e) {
throw new RuntimeException(e);
}
}
// note: Google Gson would be nicer, but since the code exists... why not use it.
private void load(String fileName) throws JSONException, IOException {
FileInputStream is;
is = new FileInputStream(new File(fileName));
String s = readFileContents(is);
//System.out.println(s);
JSONObject ppgSeries = new JSONObject(s);
this.t = json2array(ppgSeries.getJSONArray("t"));
this.x = json2array(ppgSeries.getJSONArray("x"));
//System.out.println("t=" + DVec.toString(t));
//System.out.println("x=" + DVec.toString(x));
}
private static double[] json2array(JSONArray array) throws JSONException {
double[] res = new double[array.length()];
for(int i = 0; i < res.length; i++)
res[i] = array.getDouble(i);
return res;
}
// to do: deduplicate (create helper class), see duplicate in AppData.java
private static String readFileContents(FileInputStream is) throws IOException {
BufferedReader buf = new BufferedReader(new InputStreamReader(is));
String line = buf.readLine();
StringBuilder sb = new StringBuilder();
while (line != null) {
sb.append(line).append("\n");
line = buf.readLine();
}
return sb.toString();
}
}

276
signal/src/main/java/net/heartshield/signal/Signal.java Normal file
View File

@@ -0,0 +1,276 @@
package net.heartshield.signal;
import android.util.Log;
import java.util.Arrays;
import java.util.List;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-12
*/
public class Signal {
private static final String TAG = "MeasureActivity"; // TODO
/**
* Calculate the local maxima of a (heartbeat) signal vector.
* @return one-hot 'boolean' vector of same length as `dv`.
*/
public static int[] heartbeat_localmax(double[] dv) {
if(dv.length < 3) // if cannot even diff twice, this input is too short
return IVec.zeros(dv.length);
// note: I ported this straight from Python, to get a feel for how numpy functions
// would look like in Java. There is no good reason not to just write a single loop.
//
// (the lack of operators in Java makes formula notation very painful;
// Scala may be preferred)
double[] minf = new double[]{Double.NEGATIVE_INFINITY};
double[] diffvec = DVec.stack(minf, dv, minf);
double[] diffScore = DVec.diff(DVec.sign(DVec.diff(diffvec)));
//System.out.println(DVec.toString(diffScore));
int[] e = new int[]{dv.length-1};
double[] a = DVec.get(dv, 0, dv.length-1);
double[] b = DVec.get(dv, 1, dv.length);
int[] ep = IVec.where(DVec.eq(a, b));
int[] runEndingPositions = IVec.stack(ep, e);
int[] m1 = new int[]{-1};
int[] runLengths = IVec.diff(IVec.stack(m1, runEndingPositions));
int[] runStarts = IVec.add(IVec.sub(runEndingPositions, runLengths), 1);
int[] greaterRuns = IVec.where(IVec.gt(runLengths, 1));
int[] realRunStarts = IVec.index(runStarts, greaterRuns);
int[] realRunStops = IVec.index(runEndingPositions, greaterRuns);
int[] realRunLengths = IVec.index(runLengths, greaterRuns);
int[] sta = DVec.eq(DVec.index(diffScore, realRunStarts), -1.0);
int[] sto = DVec.eq(DVec.index(diffScore, realRunStops), -1.0);
int[] maxRuns = IVec.mul(sta, sto);
double[] maxRunMiddles = DVec.sub(DVec.round(DVec.add(IVec.toDouble(IVec.index(realRunStarts, maxRuns)), DVec.mul(IVec.toDouble(IVec.index(realRunLengths, maxRuns)), 0.5))), 1.0);
int[] maxRunMiddlesIntDir = IVec.dirac(diffScore.length, DVec.toInteger(maxRunMiddles));
int[] maxima = DVec.eq(diffScore, -2.0);
int[] maximaPlus = IVec.add(maxima, maxRunMiddlesIntDir);
int[] imaxima = IVec.ge(maxima, 1);
return imaxima;
}
/** moving average filter of length N */
public static double[] running_mean(double[] x, int N) {
if(x.length < N)
throw new IllegalArgumentException("sequence too short to filter");
double[] p = new double[]{0.0};
double[] cs = DVec.cumsum(DVec.stack(p, x));
double[] a = DVec.get(cs, N, cs.length);
double[] b = DVec.get(cs, 0, cs.length - N);
double[] y = DVec.mul(DVec.sub(a, b), 1.0 / (double) N);
return y;
}
/**
* q-th percentile of x, interpolates linearly between values if necessary
* @param x input array
* @param q percentile in range [0:100]
*/
public static double percentile(double[] x, double q) {
if(q < 0.0 || q > 100.0)
throw new IllegalArgumentException("percentile q must be between [0:100]");
double[] y = DVec.sort(x);
double idx = q / 100.0 * (y.length - 1);
// linearly interpolate between neighboring data points y[i] and y[j]
int i = (int) idx;
double frac = idx - i;
int j = (i != y.length - 1) ? (i + 1) : i;
return y[i] + (y[j] - y[i]) * frac;
}
private static void throwIfNotAscending(double[] xp) {
for(int i = 0; i < xp.length - 1; i++)
if((xp[i+1] - xp[i]) <= 0.0)
throw new IllegalArgumentException("xp is not strictly monotonically ascending");
}
/** return np.interp(x, xp, fp) - boundary values are 0th order (constant) */
public static double[] interp(double[] x, double[] xp, double[] fp) {
if(xp.length != fp.length)
throw new IllegalArgumentException("xp.length != fp.length");
if(xp.length < 2)
throw new IllegalArgumentException("cannot interpolate without at least 2 data points");
throwIfNotAscending(xp);
double[] y = new double[x.length];
//double[] dxp = DVec.diff(xp);
//double[] dfp = DVec.diff(fp);
double[] slope = new double[xp.length-1];
double[] intercept = new double[xp.length-1];
for(int i = 0; i < xp.length-1; i++) {
double dxp = xp[i+1] - xp[i];
double dfp = fp[i+1] - fp[i];
slope[i] = dfp / dxp;
intercept[i] = fp[i] - slope[i] * xp[i];
}
for(int i = 0; i < x.length; i++) {
if(x[i] < xp[0]) {
y[i] = fp[0];
} else if(x[i] > xp[xp.length-1]) {
y[i] = fp[xp.length-1];
} else {
int idx = Arrays.binarySearch(xp, x[i]);
if (idx <= -1) {
// not found precisely. interpolate
idx = -idx - 1; // convert back to array index (insertion point)
// idx==0 or idx==x.length cannot happen, handled above.
//assert(idx > 0 && idx < xp.length);
idx = idx - 1; // one below insertion point (left point)
y[i] = slope[idx] * x[i] + intercept[idx];
} else {
// found precisely
y[i] = fp[idx];
}
}
}
return y;
}
/** return ibis[np.where((ibis>0.4)&(ibis<1.4))[0]] */
public static double[] goodIbis(double[] ibis) {
final double MIN_IBI = 0.4;
final double MAX_IBI = 1.4;
int[] largeEnough = DVec.lt(DVec.mul(DVec.ones(ibis.length), MIN_IBI), ibis);
int[] smallEnough = DVec.lt(ibis, DVec.mul(DVec.ones(ibis.length), MAX_IBI));
int[] goodIndices = IVec.where(IVec.mul(largeEnough, smallEnough));
return DVec.index(ibis, goodIndices);
}
/**
* @return one-hot 'boolean' vector of same length as `dv`.
*/
public static int[] beatdetect(double fps, double[] x) {
final double DIFF_THRESHOLD = 0.35;
final int N1 = (int) (0.2 * fps);
final int N2 = (int) (0.2334 * fps);
if(x.length < N1 + N2) {
//Log.i(TAG, "calculateBpm() signal too short");
return new int[0]; // signal too short
}
double[] ds = DVec.mul(Signal.running_mean(Signal.running_mean(DVec.diff(x), N1), N2), -1);
// make threshold adaptive to max slope
double T = DIFF_THRESHOLD * Signal.percentile(ds, 95); // outlier-robust max()
//Log.i(TAG, "calculateBpm() threshold T=" + T + " 95th=" + Signal.percentile(ds, 95) + " 50th=" + Signal.percentile(ds, 50) + " max=" + Signal.percentile(ds, 100));
for(int i = 0; i < ds.length; i++)
if(ds[i] < T)
ds[i] = 0;
return Signal.heartbeat_localmax(ds);
}
// TODO(david): testbed for calculateBpm() accuracy, on existing signals having a reference ECG
// (or even the PlotFilter variant and slope localmax() beatdetector, vs. getrr2()
/**
* detect beats and compute heart rate in a recent slice of detrended intensity signal.
* @param fps frames per second
* @param x intensity samples
* @return heart rate in BPM, or null if none could be computed
*/
public static Double calculateBpm(double fps, double[] x) {
final int NUM_RECENT_BEATS = 20;
final int MIN_NUM_BEATS = 6;
final double MIN_BPM = 30;
final double MAX_BPM = 180;
double[] t = DVec.linspace(0.0, x.length / fps, x.length, false);
// peak detection and bpm calculation (on smoothed signal; final bpm comes from server)
// ds = -1*running_mean(running_mean(np.diff(self.pf.detrended()), 6), 7)
int[] beats = beatdetect(fps, x);
int[] beatLoc = IVec.where(beats);
Log.i(TAG, "calculateBpm() beatLoc.length=" + beatLoc.length);
final int skipFront = 3;
final int skipBack = 1;
if(beatLoc.length < skipFront + skipBack + 1) {
//Log.i(TAG, "calculateBpm() not enough beats to skip");
return null; // not enough beats
}
// skip some early beats (noisy), and the most recent beat (may be spurious)
beatLoc = IVec.get(beatLoc, skipFront, beatLoc.length - skipBack);
double[] ibis = DVec.diff(DVec.index(t, beatLoc));
Log.i(TAG, "calculateBpm() ibis=" + DVec.toString(ibis));
ibis = Signal.goodIbis(ibis);
// only use recent beats to provide some visual sense of heart rate trends
ibis = DVec.get(ibis, Math.max(0, ibis.length - NUM_RECENT_BEATS), ibis.length);
if(ibis.length < MIN_NUM_BEATS) {
Log.i(TAG, "calculateBpm() not enough beats < MIN_NUM_BEATS");
return null; // not enough beats
}
double bpm = 60.0 / DVec.mean(ibis);
if(bpm > MAX_BPM || bpm < MIN_BPM) {
Log.i(TAG, "calculateBpm() computed invalid BPM = " + Double.toString(bpm));
return null; // invalid BPM
}
return bpm;
}
/** normalized cross-correlation of the two signals of same length */
public static double crossCorr(double[] x, double[] y) {
return DVec.sum(DVec.mul(x, y)) / Math.sqrt(DVec.sum(DVec.mul(x, x)) * DVec.sum(DVec.mul(y, y)));
}
/** two-dimensional mean of a collection of signals */
public static double[] mean(List<double[]> sigs) {
final int nsigs = sigs.size();
if(nsigs == 0)
throw new IllegalArgumentException("cannot compute mean of empty stack");
double[] y = new double[sigs.get(0).length];
for(int i = 0; i < nsigs; i++) {
double[] x = sigs.get(i);
if(x.length != y.length)
throw new IllegalArgumentException("stack entry " + i + " has different length from the rest");
for(int j = 0; j < y.length; j++)
y[j] += x[j];
}
for(int j = 0; j < y.length; j++)
y[j] /= nsigs;
return y;
}
public static double[] triangle(int ntaps) {
int htaps = ntaps/2;
double[] taps = new double[ntaps];
for(int i = 0; i <= htaps; i++)
taps[i] = ((double) i) / htaps;
for(int i = htaps+1; i < ntaps; i++)
taps[i] = 1.0 - ((double) (i - htaps)) / htaps;
taps = DVec.mul(taps, 1.0 / DVec.sum(taps));
return taps;
}
}

3
signal/src/main/res/values/strings.xml Normal file
View File

@@ -0,0 +1,3 @@
<resources>
<string name="app_name">Firdes</string>
</resources>

50
signal/src/test/java/net/heartshield/filter/AlivecorTests.java Normal file
View File

@@ -0,0 +1,50 @@
package net.heartshield.filter;
import junit.framework.TestCase;
import net.heartshield.signal.DVec;
import net.heartshield.signal.SVec;
public class AlivecorTests extends TestCase {
public void testAlivecor() {
// 1490155457_series.b_7F07F378
//short[] audio = AudioHelper.readAudio("signal/src/test/res/raw/test100hz.wav");
long t;
t = System.nanoTime();
AudioHelper ah = new AudioHelper("signal/src/test/res/raw/1490155457_series.b_7F07F378");
short[] audio = ah.readAudio();
System.out.println("readAudio() took " + (System.nanoTime()-t)/1000000 + " ms");
t = System.nanoTime();
double[] sig = SVec.toDouble(audio);
double[] norm = DVec.mul(sig, 1.0 / Math.pow(2, 15));
System.out.println("toDouble() and norm took " + (System.nanoTime()-t)/1000000 + " ms");
t = System.nanoTime();
AlivecorFilter af = new AlivecorFilter(ah.getFps());
DataSink out = new DataSink();
af.connect(out);
af.put(norm);
System.out.println("af.put() took " + (System.nanoTime()-t)/1000000 + " ms");
// af.put() took 5279 ms
t = System.nanoTime();
double[] res = out.getData();
double[] slow = DVec.downsample(res, 100);
double[] show = DVec.get(slow, 2000, slow.length);
show = DVec.sub(show, 2.0);
System.out.println("rest took " + (System.nanoTime()-t)/1000000 + " ms");
t = System.nanoTime();
//LineChart lc = new LineChart("hilbert", new XYSeries[]{ inSeries, outSeries });
LineChart lc = LineChart.plot(show);
lc.showChart();
System.out.println("plotting took " + (System.nanoTime()-t)/1000000 + " ms");
lc.waitChart();
}
}

95
signal/src/test/java/net/heartshield/filter/AudioHelper.java Normal file
View File

@@ -0,0 +1,95 @@
package net.heartshield.filter;
import java.io.File;
import java.io.IOException;
import java.io.InputStream;
import java.util.ArrayList;
import java.util.List;
import javax.sound.sampled.AudioFormat;
import javax.sound.sampled.AudioInputStream;
import javax.sound.sampled.AudioSystem;
import javax.sound.sampled.UnsupportedAudioFileException;
/**
* Helper util to load audio on PC.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-29
*/
public class AudioHelper {
AudioInputStream ais;
int frameRate;
int frameSize;
int channels;
int bps;
public AudioHelper(String fname) {
try {
ais = AudioSystem.getAudioInputStream(new File(fname));
init(ais);
} catch (IOException e) {
throw new RuntimeException(e);
} catch (UnsupportedAudioFileException e) {
throw new RuntimeException(e);
}
}
public AudioHelper(InputStream inputStream) {
try {
ais = AudioSystem.getAudioInputStream(inputStream);
init(ais);
} catch (IOException e) {
throw new RuntimeException(e);
} catch (UnsupportedAudioFileException e) {
throw new RuntimeException(e);
}
}
private void init(AudioInputStream ais) {
AudioFormat af = ais.getFormat();
frameRate = (int)af.getFrameRate();
frameSize = af.getFrameSize();
channels = af.getChannels();
bps = af.getSampleSizeInBits();
}
public int getFps() { return frameRate; }
public short[] readAudio() {
try {
return internalReadAudio();
} catch (IOException e) {
throw new RuntimeException(e);
} catch (UnsupportedAudioFileException e) {
throw new RuntimeException(e);
}
}
private short[] internalReadAudio() throws IOException, UnsupportedAudioFileException {
/*
System.out.println("Frame Rate: " + frameRate);
System.out.println("Frame Size: " + frameSize);
//System.out.println("Buffer Size: " + bufSize);
System.out.println("Channels: " + channels);
System.out.println("Bits per sample: " + bps);
*/
List<Short> buf = new ArrayList<>();
byte[] data = new byte[2048];
int bytesRead;
while((bytesRead = ais.read(data, 0, data.length)) != -1) {
for(int i = 0; i < bytesRead; i += 2)
buf.add((short) (data[i] + 256 * data[i+1]));
}
short[] arr = new short[buf.size()];
for(int i = 0; i < arr.length; i++)
arr[i] = buf.get(i);
return arr;
}
}

105
signal/src/test/java/net/heartshield/filter/AudioSampleReader.java Normal file
View File

@@ -0,0 +1,105 @@
package net.heartshield.filter;
import java.io.*;
import javax.sound.sampled.*;
public class AudioSampleReader {
private AudioInputStream audioInputStream;
private AudioFormat format;
public AudioSampleReader(File file)
throws UnsupportedAudioFileException, IOException {
audioInputStream = AudioSystem.getAudioInputStream(file);
format = audioInputStream.getFormat();
}
// Return audio format, and through it, most properties of
// the audio file: sample size, sample rate, etc.
public AudioFormat getFormat() {
return format;
}
// Return the number of samples of all channels
public long getSampleCount() {
long total = (audioInputStream.getFrameLength() *
format.getFrameSize() * 8) / format.getSampleSizeInBits();
return total / format.getChannels();
}
// Get the intervealed decoded samples for all channels, from sample
// index begin (included) to sample index end (excluded) and copy
// them into samples. end must not exceed getSampleCount(), and the
// number of samples must not be so large that the associated byte
// array cannot be allocated
public void getInterleavedSamples(long begin, long end,
double[] samples) throws IOException,
IllegalArgumentException {
long nbSamples = end - begin;
// nbBytes = nbSamples * sampleSizeinByte * nbChannels
long nbBytes = nbSamples * (format.getSampleSizeInBits() / 8) *
format.getChannels();
if (nbBytes > Integer.MAX_VALUE)
throw new IllegalArgumentException("too many samples");
// allocate a byte buffer
byte[] inBuffer = new byte[(int)nbBytes];
// read bytes from audio file
audioInputStream.read(inBuffer, 0, inBuffer.length);
// decode bytes into samples. Supported encodings are:
// PCM-SIGNED, PCM-UNSIGNED, A-LAW, U-LAW
decodeBytes(inBuffer, samples);
}
// Extract samples of a particular channel from interleavedSamples and
// copy them into channelSamples
public void getChannelSamples(int channel,
double[] interleavedSamples, double[] channelSamples) {
int nbChannels = format.getChannels();
for (int i = 0; i < channelSamples.length; i++) {
channelSamples[i] = interleavedSamples[nbChannels*i + channel];
}
}
// Convenience method. Extract left and right channels for common stereo
// files. leftSamples and rightSamples must be of size getSampleCount()
public void getStereoSamples(double[] leftSamples, double[] rightSamples)
throws IOException {
long sampleCount = getSampleCount();
double[] interleavedSamples = new double[(int)sampleCount*2];
getInterleavedSamples(0, sampleCount, interleavedSamples);
for (int i = 0; i < leftSamples.length; i++) {
leftSamples[i] = interleavedSamples[2*i];
rightSamples[i] = interleavedSamples[2*i+1];
}
}
// Private. Decode bytes of audioBytes into audioSamples
private void decodeBytes(byte[] audioBytes, double[] audioSamples) {
int sampleSizeInBytes = format.getSampleSizeInBits() / 8;
int[] sampleBytes = new int[sampleSizeInBytes];
int k = 0; // index in audioBytes
for (int i = 0; i < audioSamples.length; i++) {
// collect sample byte in big-endian order
if (format.isBigEndian()) {
// bytes start with MSB
for (int j = 0; j < sampleSizeInBytes; j++) {
sampleBytes[j] = audioBytes[k++];
}
} else {
// bytes start with LSB
for (int j = sampleSizeInBytes - 1; j >= 0; j--) {
sampleBytes[j] = audioBytes[k++];
}
}
// get integer value from bytes
int ival = 0;
for (int j = 0; j < sampleSizeInBytes; j++) {
ival += sampleBytes[j];
if (j < sampleSizeInBytes - 1) ival <<= 8;
}
// decode value
double ratio = Math.pow(2., format.getSampleSizeInBits() - 1);
double val = ((double) ival) / ratio;
audioSamples[i] = val;
}
}
}

44
signal/src/test/java/net/heartshield/filter/EveningResamplerTests.java Normal file
View File

@@ -0,0 +1,44 @@
package net.heartshield.filter;
import junit.framework.TestCase;
import net.heartshield.filter.EveningResampler;
import net.heartshield.signal.DVec;
import net.heartshield.signal.IVec;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-19
*/
public class EveningResamplerTests extends TestCase {
private static final double EPS = 1e-6;
public void setUp() throws Exception {
super.setUp();
}
public void tearDown() throws Exception {
super.tearDown();
}
public void testAddUneven() {
EveningResampler res = new EveningResampler(1.0);
double[] times = { 0.5, 2.5, 6.5, 8.0 };
double[] frames = { 0.0, 2.0, 3.0, 3.0 };
// the last entry is NOT expected, even if rounded onto that int.
double[] expect_times = new double[]{ 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 }; // , 8.0
double[] expect_frames = new double[]{ 0.0, 0.5, 1.5, 2.125, 2.375, 2.625, 2.875, 3.0 }; // , 3.0
for(int i = 0; i < frames.length; i++) {
res.addUneven(times[i], frames[i]);
}
double[] ts = res.getEvenTimes();
double[] fs = res.getEvenFrames();
if(!IVec.all(DVec.eqt(ts, expect_times, EPS)))
throw new IllegalStateException("test failed");
if(!IVec.all(DVec.eqt(fs, expect_frames, EPS)))
throw new IllegalStateException("test failed");
}
}

167
signal/src/test/java/net/heartshield/filter/IirTests.java Normal file
View File

@@ -0,0 +1,167 @@
package net.heartshield.filter;
import junit.framework.TestCase;
import net.heartshield.signal.DVec;
import net.heartshield.signal.Series;
import net.heartshield.signal.Signal;
import org.apache.commons.math3.util.MathArrays;
import org.jfree.data.xy.XYSeries;
import uk.me.berndporr.iirj.*;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-14
*/
public class IirTests extends TestCase {
public void testButterSweep() {
/*
fps = 2e3
T = 10.0
t = np.linspace(0, T, int(T * fps))
fs = np.logspace(0, np.log10(fps/2.0), int(T * fps))
sweep = np.cos(2*np.pi*fs*t)
with open('/mnt/hsh/protos/dha3/data/butter_sweep.json', 'w') as fo:
fo.write(json.dumps({'t': list(t), 'x': list(sweep)}))
*/
Series series = new Series("data/butter_sweep.json");
double FPS = 48000.0;
double[] arr = series.x;
//double[] filtered = out.getData();
double[] filtered = new double[arr.length];
Butterworth butterworth = new Butterworth();
int order = 3;
butterworth.lowPass(order, FPS, 4000.0);
long t1 = System.nanoTime();
for(int j = 0; j < arr.length; j++) {
//System.out.println(j + " " + filtered[j]);
filtered[j] = butterworth.filter(arr[j]);
}
long t2 = System.nanoTime();
// 7, 3, 5, 2.7, 3.7 ms -> roughly 3 ms
System.out.println("filtering " + arr.length + " points took " + (t2-t1) + " ns");
if(filtered.length != arr.length)
throw new IllegalStateException("output.length != input.length");
final XYSeries inSeries = new XYSeries("input");
final XYSeries outSeries = new XYSeries("output");
for(int j = 0; j < arr.length; j++) {
inSeries.add(j, arr[j]);
outSeries.add(j, filtered[j]);
}
LineChart lc = new LineChart("filter", new XYSeries[]{ outSeries, inSeries });
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
public void testLowPass() {
Series series = new Series("data/butter_sweep.json");
double FPS = 2000.0;
double[] arr = series.x;
LowPassFilter lowPass = new LowPassFilter(FPS, 100.0, 50.0);
DataSink out = new DataSink();
lowPass.connect(out);
long t1 = System.nanoTime();
for(int i = 0; i < arr.length; i += 200)
lowPass.put(DVec.get(arr, i, i + 200));
//lowPass.put(arr);
long t2 = System.nanoTime();
double[] filtered = out.getData();
System.out.println("filtering " + arr.length + " points took " + (t2-t1) + " ns");
if(filtered.length != arr.length)
throw new IllegalStateException("output.length=" + filtered.length + " != input.length=" + arr.length);
final XYSeries inSeries = new XYSeries("input");
final XYSeries outSeries = new XYSeries("output");
for(int j = 0; j < arr.length; j++) {
inSeries.add(j, arr[j]);
outSeries.add(j, filtered[j]);
}
// tw=10 -> 76 ms, 81 ms, 67 ms (but sharper)
// strange that raising tw does not seem to change much?! aah, chunk size...
LineChart lc = new LineChart("filter", new XYSeries[]{ outSeries, inSeries });
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
public static double[] convolve(double[] f, double[] g) {
//double[] y = new double[f.length];
double[] y = MathArrays.convolve(f, g);
return DVec.get(y, g.length/2, y.length - g.length/2 + 1);
}
public void testConvolve() {
Series series = new Series("data/butter_sweep.json");
double FPS = 2000.0;
double[] arr = series.x;
double[] tri = Signal.triangle(400);
long t1 = System.nanoTime();
double[] filtered = convolve(arr, tri);
long t2 = System.nanoTime();
System.out.println("filtering " + arr.length + " points took " + (t2-t1) + " ns");
if(filtered.length != arr.length)
throw new IllegalStateException("output.length=" + filtered.length + " != input.length=" + arr.length);
final XYSeries inSeries = new XYSeries("input");
final XYSeries outSeries = new XYSeries("output");
for(int j = 0; j < arr.length; j++) {
inSeries.add(j, arr[j]);
outSeries.add(j, filtered[j]);
}
// 50 ms-ish.
// note: can do with half the triangle length, will be faster then. -> 37, 32 ms
LineChart lc = new LineChart("filter", new XYSeries[]{ outSeries, inSeries });
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
}

164
signal/src/test/java/net/heartshield/filter/ImageProcessingTests.java Normal file
View File

@@ -0,0 +1,164 @@
package net.heartshield.filter;
import junit.framework.TestCase;
import net.heartshield.signal.CameraUtil;
import net.heartshield.signal.DVec;
import org.apache.commons.math3.analysis.function.Gaussian;
import org.apache.commons.math3.exception.TooManyIterationsException;
import org.apache.commons.math3.fitting.GaussianCurveFitter;
import org.apache.commons.math3.fitting.WeightedObservedPoints;
import org.jfree.data.xy.XYSeries;
import java.awt.FlowLayout;
import java.awt.Image;
import java.awt.image.BufferedImage;
import java.awt.image.Raster;
import java.io.File;
import java.io.IOException;
import javax.imageio.ImageIO;
import javax.swing.ImageIcon;
import javax.swing.JFrame;
import javax.swing.JLabel;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-11-21
*/
public class ImageProcessingTests extends TestCase {
public void testLoad() throws IOException, InterruptedException {
String[] files = "/home/david/heartshield/outsourcing/find_ellipse/example_ellipse1g.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse1.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse2.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse3g.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse3.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse4.png".split(" ");
BufferedImage img = ImageIO.read(new File(files[0]));
JFrame frame = new JFrame();
frame.getContentPane().setLayout(new FlowLayout());
frame.getContentPane().add(new JLabel(new ImageIcon(img)));
frame.pack();
frame.setVisible(true);
frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); // if you want the X button to close the app
while(frame.isVisible())
Thread.sleep(1000);
}
public void testMean() throws IOException {
String[] files = "/home/david/heartshield/outsourcing/find_ellipse/example_ellipse1g.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse1.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse2.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse3g.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse3.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse4.png".split(" ");
BufferedImage img = ImageIO.read(new File(files[0]));
int[] data = rgb(img);
//double[][] meanRow = CameraUtil.meanRow(data, img.getWidth(), img.getHeight());
double[][] meanRow = CameraUtil.meanRow(data, img.getWidth(), img.getHeight());
double[][] stdRow = CameraUtil.stdRow(data, meanRow, img.getWidth(), img.getHeight());
double[] redRow = meanRow[0];
double[] redStdRow = stdRow[0];
//LineChart lc = LineChart.plot(redRow);
final XYSeries inSeries = new XYSeries("in");
final XYSeries outSeries = new XYSeries("out");
WeightedObservedPoints obs = new WeightedObservedPoints();
for(int i = 0; i < redRow.length; i++) {
inSeries.add(i, redRow[i]);
obs.add(i, redRow[i]);
}
// length=3 {Normalization, Mean, Sigma} see http://commons.apache.org/proper/commons-math/javadocs/api-3.4/org/apache/commons/math3/fitting/GaussianCurveFitter.html
double[] parameters = GaussianCurveFitter.create().fit(obs.toList());
System.out.println("parameters.length=" + parameters.length);
System.out.println("parameters=" + DVec.toString(parameters));
Gaussian.Parametric p = new Gaussian.Parametric();
for(int i = 0; i < redRow.length; i++)
//outSeries.add(i, redRow[i]);
outSeries.add(i, p.value(i, parameters));
LineChart lc = new LineChart("brightness", new XYSeries[]{ inSeries, outSeries });
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
public void testAll() throws IOException {
String[] files = "/home/david/heartshield/tmp/ue/unellipse_5.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse1g.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse1.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse2.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse3g.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse3.png /home/david/heartshield/outsourcing/find_ellipse/example_ellipse4.png".split(" ");
for(int j = 0; j < files.length; j++) {
System.out.println(files[j]);
BufferedImage img = ImageIO.read(new File(files[j]));
int[] data = rgb(img);
//double[][] meanRow = CameraUtil.meanRow(data, img.getWidth(), img.getHeight());
double[][] meanRow = CameraUtil.meanRow(data, img.getWidth(), img.getHeight());
double[][] stdRow = CameraUtil.stdRow(data, meanRow, img.getWidth(), img.getHeight());
double[][] meanCol = CameraUtil.meanCol(data, img.getWidth(), img.getHeight());
double[][] stdCol = CameraUtil.stdCol(data, meanCol, img.getWidth(), img.getHeight());
String fn = files[j].split("/")[6];
showFit(meanRow[0], fn + " row");
showFit(meanCol[0], fn + " col");
}
}
private void showFit(double[] x, String title) {
final XYSeries inSeries = new XYSeries("in");
final XYSeries outSeries = new XYSeries("out");
WeightedObservedPoints obs = new WeightedObservedPoints();
for (int i = 0; i < x.length; i++) {
double v = x[i];
inSeries.add(i, v);
obs.add(i, v);
}
System.out.println("x=" + DVec.toString(x));
LineChart lc2 = LineChart.plot(DVec.add(x, -250.0));
lc2.showChart();
lc2.waitChart();
// length=3 {Normalization, Mean, Sigma} see http://commons.apache.org/proper/commons-math/javadocs/api-3.4/org/apache/commons/math3/fitting/GaussianCurveFitter.html
System.out.println("fitting...");
// LevenbergMarquardtOptimizer may get stuck underneath this.
//double[] parameters = GaussianCurveFitter.create().fit(obs.toList());
double[] parameters = null;
try {
//parameters = GaussianCurveFitter.create().withStartPoint(new double[]{norm,mean,std}).withMaxIterations(800).fit(obs.toList());
parameters = GaussianCurveFitter.create().withMaxIterations(800).fit(obs.toList());
System.out.println("done.");
} catch(TooManyIterationsException e) {
System.out.println("timed out.");
return;
}
Gaussian.Parametric p = new Gaussian.Parametric();
for (int i = 0; i < inSeries.getItemCount(); i++)
//outSeries.add(i, redRow[i]);
outSeries.add(i, p.value(i, parameters));
LineChart lc = new LineChart(title, new XYSeries[]{inSeries, outSeries});
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
private int[] rgb(BufferedImage image) {
int[] data = new int[image.getWidth()*image.getHeight()];
for(int j = 0, yp = 0; j < image.getHeight(); j++)
for(int i = 0; i < image.getWidth(); i++, yp++)
data[yp] = image.getRGB(i, j);
return data;
}
// CameraUtil
// double[][] meanCol(int[] rgb, int width, int height)
}

185
signal/src/test/java/net/heartshield/filter/LineChart.java Normal file
View File

@@ -0,0 +1,185 @@
package net.heartshield.filter;
import org.jfree.chart.ChartFactory;
import org.jfree.chart.ChartPanel;
import org.jfree.chart.JFreeChart;
import org.jfree.chart.axis.NumberAxis;
import org.jfree.chart.plot.PlotOrientation;
import org.jfree.chart.plot.XYPlot;
import org.jfree.chart.renderer.xy.XYLineAndShapeRenderer;
import org.jfree.data.xy.XYDataset;
import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;
import org.jfree.ui.RefineryUtilities;
import javax.swing.*;
import java.awt.*;
import java.awt.event.WindowAdapter;
import java.awt.event.WindowEvent;
/**
* LineChart window displaying a single curve.
*
* This is a test meant for PC.
* It will NOT work on Android, since JFreeChart needs Swing/AWT.
*
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-20
*/
public class LineChart extends JFrame {
private XYSeries[] mDataSeries;
private final Object mLock = new Object();
/**
* Create a LineChart window displaying a single curve.
*
* Making a data series:
* <pre>
* final XYSeries series1 = new XYSeries("First");
* series1.add(1.0, 1.0);
* series1.add(2.0, 4.0);
* </pre>
*
* @param title
* @param series
*/
public LineChart(String title, XYSeries[] series) {
super(title);
mDataSeries = series;
final XYDataset dataset = createDataset();
final JFreeChart chart = createChart(dataset);
final ChartPanel chartPanel = new ChartPanel(chart);
chartPanel.setPreferredSize(new java.awt.Dimension(500, 270));
setContentPane(chartPanel);
final JFrame window = this;
this.addWindowListener(new WindowAdapter() {
@Override
public void windowClosing(WindowEvent arg0) {
//System.out.println("in windowClosing()...");
window.setVisible(false);
synchronized(mLock) {
mLock.notify();
}
window.dispose();
//System.out.println("windowClosing() notified mLock.");
}
});
//this.setDefaultCloseOperation(JFrame.HIDE_ON_CLOSE);
}
public LineChart(String title, XYSeries dataSeries) {
this(title, new XYSeries[]{dataSeries});
}
public void showChart() {
this.pack();
RefineryUtilities.centerFrameOnScreen(this);
this.setVisible(true);
}
public void waitChart() {
//System.out.println("in waitChart()...");
synchronized(mLock) {
while(this.isVisible()) {
try {
//System.out.println("mLock.wait()...");
mLock.wait();
//System.out.println("mLock.wait() returned.");
} catch (InterruptedException e) {
e.printStackTrace();
}
}
}
//System.out.println("waitChart() done.");
}
private XYDataset createDataset() {
final XYSeriesCollection dataset = new XYSeriesCollection();
for(XYSeries series : mDataSeries)
dataset.addSeries(series);
return dataset;
}
private JFreeChart createChart(final XYDataset dataset) {
// create the chart...
final JFreeChart chart = ChartFactory.createXYLineChart(
getTitle(), // chart title
"X", // x axis label
"Y", // y axis label
dataset, // data
PlotOrientation.VERTICAL,
true, // include legend
true, // tooltips
false // urls
);
chart.setBackgroundPaint(Color.white);
final XYPlot plot = chart.getXYPlot();
plot.setBackgroundPaint(Color.lightGray);
// plot.setAxisOffset(new Spacer(Spacer.ABSOLUTE, 5.0, 5.0, 5.0, 5.0));
plot.setDomainGridlinePaint(Color.white);
plot.setRangeGridlinePaint(Color.white);
final XYLineAndShapeRenderer renderer = new XYLineAndShapeRenderer();
//renderer.setSeriesLinesVisible(0, false);
for(int i = 0; i < mDataSeries.length; i++)
renderer.setSeriesShapesVisible(i, false);
plot.setRenderer(renderer);
// change the auto tick unit selection to integer units only...
final NumberAxis rangeAxis = (NumberAxis) plot.getRangeAxis();
rangeAxis.setStandardTickUnits(NumberAxis.createIntegerTickUnits());
return chart;
}
public static void main(final String[] args) {
final XYSeries series2 = new XYSeries("Series title");
series2.add(1.0, 5.0);
series2.add(2.0, 7.0);
series2.add(3.0, 6.0);
series2.add(4.0, 8.0);
series2.add(5.0, 4.0);
series2.add(6.0, 4.0);
series2.add(7.0, 2.0);
series2.add(8.0, 1.0);
LineChart lc = new LineChart("Plot title", series2);
//System.out.println("showChart()");
lc.showChart();
//System.out.println("wait()...");
lc.waitChart();
System.out.println("finished");
// exit the Event Dispatch Thread, or EDT
//System.exit(0);
// not necessary if window.dispose() is called
}
public static LineChart plot(double[] x) {
final XYSeries series = new XYSeries("data");
for(int i = 0; i < x.length; i++)
series.add(i, x[i]);
return new LineChart("data", series);
}
public static LineChart plot(double[] t, double[] x) {
final XYSeries series = new XYSeries("data");
if(t.length != x.length)
throw new IllegalArgumentException("t.length=" + t.length + " != x.length=" + x.length);
for(int i = 0; i < x.length; i++)
series.add(t[i], x[i]);
return new LineChart("data", series);
}
}

293
signal/src/test/java/net/heartshield/filter/PcgFilterTests.java Normal file
View File

@@ -0,0 +1,293 @@
package net.heartshield.filter;
import junit.framework.TestCase;
import net.heartshield.signal.DVec;
import org.jfree.data.xy.XYSeries;
import java.io.File;
import java.io.IOException;
import javax.sound.sampled.UnsupportedAudioFileException;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-15
*/
public class PcgFilterTests extends TestCase {
public void testEnvelopeFilter() throws IOException, UnsupportedAudioFileException {
AudioData audio = readAudio("/mnt/hsh/data/pcg/2017-05-15/2017-05-14.wav");
double[] left = audio.left;
double frontFps = audio.fps;
System.out.println("min=" + DVec.min(left) + " max=" + DVec.max(left));
// 1e6 .. 2.3e6
PcgEnvelopeFilter pcg = new PcgEnvelopeFilter(frontFps);
DataSink sink = new DataSink();
pcg.connect(sink);
long t1 = System.nanoTime();
//pcg.put(left);
int step = 4800;
for(int i = 0; i < left.length; i += step) {
//System.out.println("put(i=" + i + ")");
pcg.put(DVec.get(left, i, Math.min(i + step, left.length)));
}
long t2 = System.nanoTime();
// 170-190 ms
System.out.println("time elapsed: " + (t2-t1)/1e6 + " ms // including DataSink.put()");
double[] filtered = sink.getData();
System.out.println("in.length=" + left.length + " out.length=" + filtered.length);
filtered = DVec.get(filtered, 50000, 60000);
System.out.println("filtered min=" + DVec.min(filtered) + " max=" + DVec.max(filtered));
LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
public void testBeatDetector() throws IOException, UnsupportedAudioFileException {
AudioData audio = readAudio("../app/src/main/res/raw/canned_pcg_short.wav");
double[] left = audio.left;
double frontFps = audio.fps;
System.out.println("min=" + DVec.min(left) + " max=" + DVec.max(left));
// 1e6 .. 2.3e6
PcgBeatDetector pcg = new PcgBeatDetector(frontFps);
//TestBeatDetector pcg = new TestBeatDetector(frontFps);
DataSink sinkPerc = new DataSink();
DataSink sinkOut = new DataSink();
//pcg.connect(sink);
//pcg.percentileFilter.connect(sinkPerc);
//pcg.output.connect(sinkOut);
//sinkPerc.put(left);
pcg.connectPercentiles(sinkPerc);
pcg.connect(sinkOut);
long t1 = System.nanoTime();
//pcg.put(left);
//int step = 48000;
int step = 4800;
// 4800/27 = 177.7777777777778 / we have len=200 samples (why?)
// so PERCENTILE_WIN_SECS = 4.0 ^= 7000 samples... of course, no percentiles.
for(int i = 0; i < left.length; i += step) {
//System.out.println("put(i=" + i + ")");
pcg.put(DVec.get(left, i, Math.min(i + step, left.length)));
}
long t2 = System.nanoTime();
// full put() -> 332.265274 ms
// with step=4800 -> 3350.72893 ms
// now better?! 2400.01039 ms
System.out.println("time elapsed: " + (t2-t1)/1e6 + " ms // including DataSink.put()");
/*
double[] filtered = sink1.getData();
filtered = DVec.get(filtered, 50000 * 24, 60000 * 24);
System.out.println("filtered min=" + DVec.min(filtered) + " max=" + DVec.max(filtered));
LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();*/
double[] perc = sinkPerc.getData();
double[] out = sinkOut.getData();
final XYSeries percSeries = new XYSeries("in");
final XYSeries outSeries = new XYSeries("output");
final double clipMax = 50.0;
//for(int j = 0; j < perc.length; j++) {
//for(int j = 40000*27; j < 85000*27; j++) {
for(int j = 0; j < perc.length; j++) {
percSeries.add(j, perc[j] > clipMax ? clipMax : perc[j]);
outSeries.add(j, out[j]);
}
double npeaks = DVec.sum(DVec.get(out, 40000*27, Math.min(85000*27, out.length)));
// 58.0
System.out.println("npeaks = " + npeaks);
// in.length=2 723 840 out.length=2 710 224 WHY? TODO
// 13616 missing
// 2723840/4800 = 567 steps
// ah, moving averages maybe?
System.out.println("in.length=" + left.length + " out.length=" + out.length);
LineChart lc = new LineChart("beatdet", new XYSeries[]{ percSeries, outSeries });
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
static class AudioData {
public double fps;
public double[] left;
public AudioData(double fps, double[] left) {
this.fps = fps;
this.left = left;
}
}
private AudioData readAudio(String fname) throws IOException, UnsupportedAudioFileException {
AudioSampleReader reader = new AudioSampleReader(new File(fname));
int len = (int) reader.getSampleCount();
double[] left = new double[len];
//double[] right = new double[len];
//reader.getStereoSamples(left, right);
double[] interleaved = new double[len*reader.getFormat().getChannels()];
reader.getInterleavedSamples(0, len, interleaved);
reader.getChannelSamples(0, interleaved, left);
return new AudioData(reader.getFormat().getFrameRate(), left);
}
private double[] getBeep() throws IOException, UnsupportedAudioFileException {
//throw new IllegalStateException(new File("").getAbsolutePath());
// from prevent-android/signal
AudioData audio = readAudio("../app/src/main/res/raw/beep.wav");
return audio.left;
}
public void testBeatDetectorPerformance() throws IOException, UnsupportedAudioFileException {
AudioData audio = readAudio("/mnt/hsh/data/pcg/2017-05-15/2017-05-14.wav");
double[] left = audio.left;
double frontFps = audio.fps;
System.out.println("min=" + DVec.min(left) + " max=" + DVec.max(left));
// 1e6 .. 2.3e6
double[] beep = getBeep();
//double[] beep = new double[(int) (0.08 * frontFps)]; // 3840 -> hopefully 8196-point FFT with the step=4800...
FirFilter beeper = new FirFilter(frontFps, beep);
//FirFilterRunning beeper = new FirFilterRunning(beep);
//mPcgDetector.connect(mBeeper);
beeper.connect(new IRFilterBlock() {
@Override
public void put(double[] x) {
// ignore
}
});
PcgBeatDetector pcg = new PcgBeatDetector(frontFps);
//TestBeatDetector pcg = new TestBeatDetector(frontFps);
DataSink sinkPerc = new DataSink();
DataSink sinkOut = new DataSink();
//pcg.connect(sink);
//pcg.percentileFilter.connect(sinkPerc);
//pcg.output.connect(sinkOut);
sinkPerc.put(left);
//pcg.connect(sinkOut);
pcg.connect(beeper);
long t1 = System.nanoTime();
//pcg.put(left);
//int step = 48000;
int step = 4800;
// 4800/27 = 177.7777777777778 / we have len=200 samples (why?)
// so PERCENTILE_WIN_SECS = 4.0 ^= 7000 samples... of course, no percentiles.
for(int i = 0; i < left.length; i += step) {
//System.out.println("put(i=" + i + ")");
pcg.put(DVec.get(left, i, Math.min(i + step, left.length)));
}
long t2 = System.nanoTime();
// 2528.27118 ms
// 2321.839134 ms
System.out.println("time elapsed: " + (t2-t1)/1e6 + " ms // including DataSink.put()");
}
static class TestBeatDetector extends RFilterBlock {
private PcgEnvelopeFilter mPcgFilter;
private Splitter mSplitter = new Splitter();
private Diff mDiff = new Diff(); // could also just have been FIR({-1, 1})
private FirFilterRunning mMean1;
private FirFilterRunning mMean2;
private PcgBeatDetector.DiracUpsampler mUpsampler;
// both are lower rate (downsampled)
public PercentileFilter percentileFilter;
public FirFilterRunning output;
public TestBeatDetector(double frontFps) {
mPcgFilter = new PcgEnvelopeFilter(frontFps);
double fps = mPcgFilter.getFps();
//mPcgFilter.connect(mSplitter);
mPcgFilter.connect(mDiff);
final double PERCENTILE_WIN_SECS = 4.0;
final double PERCENTILE = 85;
final double PERCENTILE_INIT = Double.MAX_VALUE; // ignore all beats until PercentileFilter has settled
percentileFilter = new PercentileFilter((int) (PERCENTILE_WIN_SECS * fps), PERCENTILE_INIT, PERCENTILE);
output = new FirFilterRunning(new double[]{1.0});
//final int N1 = (int) (0.08 * fps);
final int N2 = (int) (0.08 * fps);
//mMean1 = new FirFilterRunning(DVec.ones(N1));
mMean2 = new FirFilterRunning(DVec.ones(N2));
//mDiff.connect(mMean1);
mDiff.connect(mMean2);
//mPcgFilter.connect(mMean1);
//mMean1.connect(mMean2);
mMean2.connect(mSplitter);
percentileFilter.connect(mPercentiles);
mSplitter.connect(percentileFilter); // must be before mBeatDetector, which uses mPercentiles.buf
mSplitter.connect(output);
}
@Override
protected double[] batch(double[] x) {
throw new IllegalStateException();
}
static class FilterBuf implements IRFilterBlock {
public double[] buf;
@Override
public void put(double[] x) {
buf = x;
}
}
private FilterBuf mPercentiles = new FilterBuf();
@Override
public void put(double[] x) {
mPcgFilter.put(x);
}
@Override
public void connect(IRFilterBlock consumer) {
throw new IllegalStateException("connect to percentileFilter or output instead");
}
}
}

429
signal/src/test/java/net/heartshield/filter/PlotFilterTests.java Normal file
View File

@@ -0,0 +1,429 @@
package net.heartshield.filter;
import junit.framework.TestCase;
import net.heartshield.signal.DVec;
import net.heartshield.signal.FirFilters;
import net.heartshield.signal.Signal;
import org.apache.commons.math3.complex.Complex;
import org.jfree.data.xy.XYSeries;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-20
*/
public class PlotFilterTests extends TestCase {
public void testPlotSine() {
double[] x = new double[100];
final double F = 1.0, FPS = 30.0;
for(int i = 0; i < x.length; i++)
x[i] = Math.sin(2*Math.PI*F/FPS * i);
x = DVec.add(x, DVec.linspace(0.0, 3.0, x.length, false));
LineChart lc = LineChart.plot(x);
//LineChart lc = LineChart.plot(DVec.linspace(0.0, 3.0, x.length, false));
lc.showChart();
lc.waitChart();
}
public void testDetrend() {
double[] t = new double[100];
double[] x = new double[100];
final double F = 1.0, FPS = 30.0;
for(int i = 0; i < x.length; i++) {
t[i] = (double) i / FPS;
x[i] = Math.sin(2 * Math.PI * F * t[i]);
}
x = DVec.add(x, DVec.linspace(0.0, 3.0, x.length, false));
final int PLOT_SAMPLES = 100;
final double TOTAL_TIME = 15.0;
PlotFilter pf = new PlotFilter(FPS, (int) (TOTAL_TIME * FPS + 1), PLOT_SAMPLES);
for(int i = 0; i < x.length; i++)
pf.addUneven(t[i], x[i]);
// to do: why is getPlotTimes() off by a few samples?! (time starts at 0.1 sec)
PlotFilter.Series ser = pf.getPlotSeries();
double[] times = ser.t;
//double[] frames = pf.getEvenFrames();
double[] frames = ser.x;
//LineChart lc = LineChart.plot(frames);
LineChart lc = LineChart.plot(times, frames);
lc.showChart();
lc.waitChart();
}
public void testSigGen() {
double FPS = 30.0;
////
double mStartTime = ((double) System.currentTimeMillis()) / 1000;
double dt = 1.0 / FPS;
double duration = 75.0;
double f = 1.0;
double end = mStartTime + duration;
double[] arr = new double[(int) (duration*FPS+1)];
int i = 0;
for(double t = mStartTime; t < end; t += dt) {
double val = Math.sin(2 * Math.PI * f * t);
arr[i++] = val;
}
LineChart lc = LineChart.plot(arr);
lc.showChart();
lc.waitChart();
}
public void testHilbertImag() {
double[] arr = new double[200];
arr[100] = 1;
HilbertImag hi = new HilbertImag();
DataSink out = new DataSink();
hi.connect(out);
hi.put(arr);
double[] res = out.getData();
if(res.length != arr.length)
throw new IllegalStateException("output.length=" + res.length + " != input.length=" + arr.length);
for(int j = 0; j < res.length; j++) {
System.out.println(j + " " + res[j]);
}
final XYSeries inSeries = new XYSeries("impulse");
final XYSeries outSeries = new XYSeries("response");
for(int j = 0; j < arr.length; j++) {
inSeries.add(j, arr[j]);
outSeries.add(j, res[j]);
}
LineChart lc = new LineChart("hilbert", new XYSeries[]{ inSeries, outSeries });
lc.showChart();
lc.waitChart();
}
public void testHilbertImpulse() {
double[] arr = new double[200];
arr[100] = 1;
Hilbert hi = new Hilbert();
CDataSink out = new CDataSink();
hi.connect(out);
hi.put(arr);
Complex[] res = out.getData();
if(res.length != arr.length)
throw new IllegalStateException("output.length=" + res.length + " != input.length=" + arr.length);
for(int j = 0; j < res.length; j++) {
System.out.println(j + " " + res[j]);
}
final XYSeries realSeries = new XYSeries("real");
final XYSeries imagSeries = new XYSeries("imag");
for(int j = 0; j < res.length; j++) {
realSeries.add(j, res[j].getReal());
imagSeries.add(j, res[j].getImaginary());
}
LineChart lc = new LineChart("hilbert", new XYSeries[]{ realSeries, imagSeries });
lc.showChart();
lc.waitChart();
}
public void testHilbert() {
double FPS = 30.0;
////
//double mStartTime = ((double) System.currentTimeMillis()) / 1000;
double mStartTime = 0.0;
double dt = 1.0 / FPS;
double duration = 7.5;
double f = 1.0;
double end = mStartTime + duration;
double[] arr = new double[(int) (duration*FPS+1)];
int i = 0;
for(double t = mStartTime; t < end; t += dt) {
double val = Math.sin(2 * Math.PI * f * t);
arr[i++] = val;
}
Hilbert hilbert = new Hilbert();
CDataSink out = new CDataSink();
hilbert.connect(out);
hilbert.put(arr);
Complex[] analytic = out.getData();
if(analytic.length != arr.length)
throw new IllegalStateException("output.length=" + analytic.length + " != input.length=" + arr.length);
final XYSeries realSeries = new XYSeries("real");
final XYSeries imagSeries = new XYSeries("imag");
for(int j = 0; j < analytic.length; j++) {
realSeries.add(j, analytic[j].getReal());
imagSeries.add(j, analytic[j].getImaginary());
}
LineChart lc = new LineChart("hilbert", new XYSeries[]{ realSeries, imagSeries });
lc.showChart();
lc.waitChart();
}
public void testLowpass() {
double FPS = 30.0;
////
//double mStartTime = ((double) System.currentTimeMillis()) / 1000;
double mStartTime = 0.0;
double dt = 1.0 / FPS;
double duration = 75.0;
double f = 0.1;
double end = mStartTime + duration;
double[] arr = new double[(int) (duration*FPS+1)];
int i = 0;
for(double t = mStartTime; t < end; t += dt) {
double val = Math.sin(2 * Math.PI * f * t);
arr[i++] = val;
}
LowPassFilter lowPass = new LowPassFilter(FPS, 2*f, f*0.5);
DataSink out = new DataSink();
lowPass.connect(out);
lowPass.put(arr);
double[] filtered = out.getData();
if(filtered.length != arr.length)
throw new IllegalStateException("output.length != input.length");
//LineChart lc = LineChart.plot(arr);
LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
public void testHighpass() {
double FPS = 30.0;
////
//double mStartTime = ((double) System.currentTimeMillis()) / 1000;
double mStartTime = 0.0;
double dt = 1.0 / FPS;
double duration = 75.0;
double f = 0.1;
double end = mStartTime + duration;
double[] arr = new double[(int) (duration*FPS+1)];
int i = 0;
for(double t = mStartTime; t < end; t += dt) {
double val = Math.sin(2 * Math.PI * f * t);
arr[i++] = val;
}
HighPassFilter highPass = new HighPassFilter(FPS, 0.3*f, f*0.5);
System.out.println("delay=" + highPass.getDelay());
DataSink out = new DataSink();
highPass.connect(out);
highPass.put(arr);
double[] filtered = out.getData();
if(filtered.length != arr.length)
throw new IllegalStateException("output.length != input.length");
//LineChart lc = LineChart.plot(arr);
LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
public void testBandpassResponse() {
double FPS = 30.0;
double f = 0.1;
double[] imp = FirFilters.bandPass(FPS, 0.7*f, 1.3*f, f*0.2);
LineChart lc = LineChart.plot(imp);
lc.showChart();
lc.waitChart();
}
public void testBandpass() {
double FPS = 30.0;
////
//double mStartTime = ((double) System.currentTimeMillis()) / 1000;
double mStartTime = 0.0;
double dt = 1.0 / FPS;
double duration = 75.0;
double f = 0.1;
double end = mStartTime + duration;
double[] arr = new double[(int) (duration*FPS+1)];
int i = 0;
for(double t = mStartTime; t < end; t += dt) {
double val = Math.sin(2 * Math.PI * f * t) + 3 * Math.sin(2 * Math.PI * (f * 0.05) * t) + 2 * Math.sin(2 * Math.PI * (f * 3.0) * t);
//double val = Math.sin(2 * Math.PI * f * t);
arr[i++] = val;
}
BandPassFilter bandPass = new BandPassFilter(FPS, 0.2*f, 1.8*f, f*0.5);
DataSink out = new DataSink();
bandPass.connect(out);
bandPass.put(arr);
double[] filtered = out.getData();
for(int j = 0; j < filtered.length; j++) {
System.out.println(j + " " + filtered[j]);
}
if(filtered.length != arr.length)
throw new IllegalStateException("output.length != input.length");
final XYSeries inSeries = new XYSeries("input");
final XYSeries outSeries = new XYSeries("output");
for(int j = 0; j < arr.length; j++) {
inSeries.add(j, arr[j]);
outSeries.add(j, filtered[j]);
}
LineChart lc = new LineChart("bandpass", new XYSeries[]{ inSeries, outSeries });
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
public void testBandreject() {
double FPS = 30.0;
////
//double mStartTime = ((double) System.currentTimeMillis()) / 1000;
double mStartTime = 0.0;
double dt = 1.0 / FPS;
double duration = 75.0;
double f = 0.1;
double end = mStartTime + duration;
double[] arr = new double[(int) (duration*FPS+1)];
int i = 0;
for(double t = mStartTime; t < end; t += dt) {
double val = 3 * Math.sin(2 * Math.PI * f * t) + Math.sin(2 * Math.PI * (f * 0.05) * t) + 2 * Math.sin(2 * Math.PI * (f * 2.75) * t);
//double val = Math.sin(2 * Math.PI * f * t);
arr[i++] = val;
}
BandRejectFilter bandReject = new BandRejectFilter(FPS, 0.2*f, 1.8*f, f*0.5);
DataSink out = new DataSink();
bandReject.connect(out);
bandReject.put(arr);
double[] filtered = out.getData();
for(int j = 0; j < filtered.length; j++) {
System.out.println(j + " " + filtered[j]);
}
if(filtered.length != arr.length)
throw new IllegalStateException("output.length != input.length");
final XYSeries inSeries = new XYSeries("input");
final XYSeries outSeries = new XYSeries("output");
for(int j = 0; j < arr.length; j++) {
inSeries.add(j, arr[j]);
outSeries.add(j, filtered[j]);
}
LineChart lc = new LineChart("bandreject", new XYSeries[]{ inSeries, outSeries });
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
public void testImpulse() {
double[] taps = Signal.triangle(200);
double[] sig = new double[2000];
sig[1000] = 1.0;
FirFilterRunning fir = new FirFilterRunning(taps);
double[] y = fir.batch(sig);
final XYSeries inSeries = new XYSeries("input");
final XYSeries outSeries = new XYSeries("output");
for(int j = 0; j < sig.length; j++) {
inSeries.add(j, sig[j]);
outSeries.add(j, y[j] * 100);
}
LineChart lc = new LineChart("running_fir", new XYSeries[]{ inSeries, outSeries });
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
public void testSine() {
double[] taps = Signal.triangle(200);
double[] sig = new double[2000];
double FPS = 200.0;
double F = 0.1;
for(int i = 0; i < sig.length; i++)
sig[i] = Math.sin(2*Math.PI*F/FPS*i);
FirFilterRunning fir = new FirFilterRunning(taps);
double[] y = fir.batch(sig);
final XYSeries inSeries = new XYSeries("input");
final XYSeries outSeries = new XYSeries("output");
for(int j = 0; j < sig.length; j++) {
inSeries.add(j, sig[j]);
outSeries.add(j, y[j]);
}
LineChart lc = new LineChart("running_fir", new XYSeries[]{ inSeries, outSeries });
//LineChart lc = LineChart.plot(filtered);
lc.showChart();
lc.waitChart();
}
}

21
signal/src/test/java/net/heartshield/filter/PpgSigqualTests.java Normal file
View File

@@ -0,0 +1,21 @@
package net.heartshield.filter;
import junit.framework.TestCase;
import net.heartshield.signal.DVec;
import net.heartshield.signal.Series;
import org.json.JSONException;
public class PpgSigqualTests extends TestCase {
public void testLoadSeries() throws JSONException {
Series series = new Series("data/1480348240_B319CE46_meta.ppg_raw_even.json");
System.out.println("t=" + DVec.toString(series.t));
System.out.println("x=" + DVec.toString(series.x));
// test if evenly spaced
double[] dt = DVec.diff(series.t);
System.out.println("dt=" + DVec.toString(dt));
}
}

88
signal/src/test/java/net/heartshield/signal/AudioDecoderTests.java Normal file
View File

@@ -0,0 +1,88 @@
package net.heartshield.signal;
import junit.framework.TestCase;
import java.io.BufferedInputStream;
import java.io.DataInputStream;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.IOException;
import java.nio.file.Path;
import java.nio.file.Paths;
import java.util.ArrayList;
import java.util.List;
import javax.sound.sampled.AudioFormat;
import javax.sound.sampled.AudioInputStream;
import javax.sound.sampled.AudioSystem;
import javax.sound.sampled.DataLine;
import javax.sound.sampled.SourceDataLine;
import javax.sound.sampled.UnsupportedAudioFileException;
public class AudioDecoderTests extends TestCase {
public void testA() throws IOException, UnsupportedAudioFileException {
//System.out.println();
Path currentRelativePath = Paths.get("");
String s = currentRelativePath.toAbsolutePath().toString();
System.out.println("Current relative path is: " + s);
System.out.println("length=" + readAudio("signal/src/test/res/raw/test100hz.wav").length);
/*
* Frame Rate: 44100
* Frame Size: 2
* Channels: 1
* Bits per sample: 16
*/
;
}
private int[] delArray(int[] arr, int i) {
int[] newMeasurements = new int[arr.length-1];
System.arraycopy(arr, 0, newMeasurements, 0, i);
if(i < newMeasurements.length)
System.arraycopy(arr, i+1, newMeasurements, i, (newMeasurements.length - i));
return newMeasurements;
}
public void testArrayRestr() {
int[] start = new int[]{1, 2, 3, 4};
int[] res = delArray(start, 1);
int[] exp = new int[]{1, 3, 4};
if(!IVec.all(IVec.eq(res, exp)))
throw new IllegalStateException("failed delArray() test");
}
private short[] readAudio(String fname) throws IOException, UnsupportedAudioFileException {
AudioInputStream ais = AudioSystem.getAudioInputStream(new File(fname));
AudioFormat af = ais.getFormat();
//DataLine.Info info = new DataLine.Info(SourceDataLine.class, af);
int frameRate = (int)af.getFrameRate();
int frameSize = af.getFrameSize();
int channels = af.getChannels();
int bps = af.getSampleSizeInBits();
System.out.println("Frame Rate: " + frameRate);
System.out.println("Frame Size: " + frameSize);
//System.out.println("Buffer Size: " + bufSize);
System.out.println("Channels: " + channels);
System.out.println("Bits per sample: " + bps);
List<Short> buf = new ArrayList<>();
byte[] data = new byte[2048];
int bytesRead;
while((bytesRead = ais.read(data, 0, data.length)) != -1) {
for(int i = 0; i < bytesRead; i += 2)
buf.add((short) (data[i] + 256 * data[i+1]));
}
short[] arr = new short[buf.size()];
for(int i = 0; i < arr.length; i++)
arr[i] = buf.get(i);
return arr;
}
}

25
signal/src/test/java/net/heartshield/signal/FirdesTests.java Normal file
View File

@@ -0,0 +1,25 @@
package net.heartshield.signal;
import junit.framework.TestCase;
import net.heartshield.filter.LineChart;
public class FirdesTests extends TestCase {
private static final double EPS = 1e-6;
public void testHilbert() {
double[] arr = FirFilters.hilbert(65);
for(int i = 0; i < arr.length; i++)
System.out.println(i + " " + arr[i]);
double[] ref = new double[]{0.0, -0.020952552556991577, 0.0, -0.022397557273507118, 0.0, -0.024056635797023773, 0.0, -0.02598116546869278, 0.0, -0.028240399435162544, 0.0, -0.030929960310459137, 0.0, -0.0341857448220253, 0.0, -0.03820759803056717, 0.0, -0.04330194741487503, 0.0, -0.04996378347277641, 0.0, -0.05904810503125191, 0.0, -0.07216990739107132, 0.0, -0.09278988093137741, 0.0, -0.1299058347940445, 0.0, -0.21650972962379456, 0.0, -0.6495291590690613, 0.0, 0.6495291590690613, 0.0, 0.21650972962379456, 0.0, 0.1299058347940445, 0.0, 0.09278988093137741, 0.0, 0.07216990739107132, 0.0, 0.05904810503125191, 0.0, 0.04996378347277641, 0.0, 0.04330194741487503, 0.0, 0.03820759803056717, 0.0, 0.0341857448220253, 0.0, 0.030929960310459137, 0.0, 0.028240399435162544, 0.0, 0.02598116546869278, 0.0, 0.024056635797023773, 0.0, 0.022397557273507118, 0.0, 0.020952552556991577, 0.0};
if(!IVec.all(DVec.eqt(arr, ref, EPS)))
throw new IllegalStateException("test failed");
LineChart lc = LineChart.plot(arr);
lc.showChart();
lc.waitChart();
}
}

55
signal/src/test/java/net/heartshield/signal/RunningQualityTest.java Normal file
View File

@@ -0,0 +1,55 @@
package net.heartshield.signal;
import junit.framework.TestCase;
import net.heartshield.filter.LineChart;
import java.util.ArrayList;
import java.util.List;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-05-11
*/
public class RunningQualityTest extends TestCase {
private class TestResults {
int nbeats = 0;
List<Double> posCorrCoeffs = new ArrayList<>();
}
public void testRQBeatDet() {
final double FPS = 30.0;
final TestResults results = new TestResults();
RunningQuality rq = new RunningQuality(FPS);
rq.setListener(new RunningQuality.BeatListener() {
@Override
public void onLocked(int startIdx) {}
@Override
public void onBeat(int startIdx, boolean goodBeat, double posCorrCoeff) {
System.out.println("onBeat(posCorrCoeff=" + posCorrCoeff);
results.nbeats++;
results.posCorrCoeffs.add(posCorrCoeff);
}
});
Series series = new Series("data/1480348240_B319CE46_meta.ppg_raw_even.json");
// TODO: detrend
/*
LineChart lc = LineChart.plot(series.t, series.x);
lc.showChart();
lc.waitChart();
*/
// note: it works even on raw (undetrended) data.
// but in real scenario, we need to detrend, otherwise a signal with around 230 DC and only <1 AC will always have very good correlation
for(int i = 0; i < series.x.length - 30; i += 30)
rq.append(i / FPS, DVec.get(series.x, i, i + 30));
assertTrue("not enough beats detected", results.nbeats > 40);
assertTrue("too many beats detected", results.nbeats < 180);
}
}

134
signal/src/test/java/net/heartshield/signal/SignalTests.java Normal file
View File

@@ -0,0 +1,134 @@
package net.heartshield.signal;
import junit.framework.TestCase;
import net.heartshield.filter.LineChart;
import net.heartshield.filter.PlotFilter;
import net.heartshield.signal.DVec;
import net.heartshield.signal.IVec;
import net.heartshield.signal.Signal;
import java.io.FileInputStream;
import java.io.IOException;
/**
* @author David Madl (git@abanbytes.eu)
* @date 2017-03-12
*/
public class SignalTests extends TestCase {
private static final double EPS = 1e-6;
public void setUp() throws Exception {
super.setUp();
}
public void tearDown() throws Exception {
super.tearDown();
}
public void testDiff() {
double[] a = new double[]{ 1., 2., 5. };
double[] d = DVec.diff(a);
double[] expect_d = new double[]{ 1., 3. };
if(!IVec.all(DVec.eqt(d, expect_d, EPS)))
throw new IllegalStateException("test failed");
}
public void testLocalmax() {
double[] a = new double[]{ 1., 2., 5., 3., 0., 4., -1. };
int[] expect_lm = new int[]{ 0, 0, 1, 0, 0, 1, 0 };
int[] lm = Signal.heartbeat_localmax(a);
//System.out.println("heartbeat_localmax="+ IVec.toString(lm));
if(!IVec.all(IVec.eq(lm, expect_lm)))
throw new IllegalStateException("test failed");
}
public void testRunningMean() {
double[] a = new double[]{ 1., 2., 4., 3., 2., 1. };
double[] expect_m = new double[]{ 1.5, 3., 3.5, 2.5, 1.5 };
double[] m = Signal.running_mean(a, 2);
//System.out.println("heartbeat_localmax="+ IVec.toString(lm));
if(!IVec.all(DVec.eqt(m, expect_m, EPS)))
throw new IllegalStateException("test failed");
}
public void testPercentile() {
double[] a = new double[]{ 10, 7, 4, 3, 2, 1 };
double[] q = new double[]{ 0.0, 100.0, 60.0, 80.0, 50.0, 66.666666666 };
double[] expect_p = new double[]{ 1.0, 10.0, 4.0, 7.0, 3.5, 4.99999999 };
double[] p = new double[q.length];
for(int i = 0; i < q.length; i++)
p[i] = Signal.percentile(a, q[i]);
if(!IVec.all(DVec.eqt(p, expect_p, EPS)))
throw new IllegalStateException("test failed");
}
private static final double FPS = 30.0;
/** rough detrending as done by PlotFilter (more for visuals and is not precise [optimized for less filter delay]) */
private Series plotDetrend(Series series) {
final int PLOT_SAMPLES = 100;
final double TOTAL_TIME = 15.0;
PlotFilter pf = new PlotFilter(FPS, (int) (TOTAL_TIME * FPS + 1), PLOT_SAMPLES);
for(int i = 0; i < series.x.length; i++)
pf.addUneven(series.t[i], series.x[i]);
PlotFilter.Series ser = pf.getRecentSeries(pf.size());
// to do: why is getPlotTimes() off by a few samples?! (time starts at 0.1 sec)
double[] times = ser.t;
//double[] frames = pf.getEvenFrames();
double[] frames = ser.x;
return new Series(times, frames);
}
public void testDetrend1() {
Series series = new Series("data/1480348240_B319CE46_meta.ppg_raw_even.json");
Series detrended = plotDetrend(series);
detrended.t = DVec.get(detrended.t, (int) (5.0 * FPS), detrended.t.length);
detrended.x = DVec.get(detrended.x, (int) (5.0 * FPS), detrended.x.length);
//LineChart lc = LineChart.plot(series.t, series.x);
LineChart lc = LineChart.plot(detrended.t, detrended.x);
lc.showChart();
lc.waitChart();
}
public void testBeatDetector() {
Series series = new Series("data/1480348240_B319CE46_meta.ppg_raw_even.json");
Series detrended = plotDetrend(series);
//double[] sig = detrended.x;
double[] sig = DVec.get(detrended.x, 0, 450);
Double bpm = Signal.calculateBpm(FPS, sig);
assertNotNull(bpm);
System.out.println("bpm = " + bpm);
// full one: bpm = 80.49580151495677
// [:600]: bpm = 70.47058823529389
// [:450]: bpm = 69.34362934362998
}
public void testInterp() {
double[] xp = new double[]{ 1., 2., 3. };
double[] fp = new double[]{ 1., 2., 2. };
double[] x = new double[]{ 0.5, 1.5, 2.5, 3.5 };
double[] expect_y = new double[]{ 1. , 1.5, 2. , 2. };
double[] y = Signal.interp(x, xp, fp);
System.out.println("interp(x,xp,fp)="+ DVec.toString(y));
if(!IVec.all(DVec.eqt(y, expect_y, EPS)))
throw new IllegalStateException("test failed");
}
}

BIN
signal/src/test/res/raw/1490155457_series.b_7F07F378 Normal file
View File

Binary file not shown.

BIN
signal/src/test/res/raw/test100hz.wav Normal file
View File

Binary file not shown.