diff --git a/ASICPhaseCalculator/Source/HTransformers.cpp b/ASICPhaseCalculator/Source/HTransformers.cpp deleted file mode 100644 index c480ce1..0000000 --- a/ASICPhaseCalculator/Source/HTransformers.cpp +++ /dev/null @@ -1,255 +0,0 @@ -/* ------------------------------------------------------------------- - -This file is part of a plugin for the Open Ephys GUI -Copyright (C) 2018 Translational NeuroEngineering Lab - ------------------------------------------------------------------- - -This program is free software: you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -This program is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with this program. If not, see . - -*/ - -#include "HTransformers.h" - -namespace PhaseCalculator -{ - namespace - { - /** Helper to allow assigning to each array one band at a time (within the constructor) **/ - struct HilbertInfo - { - const String delta{ L"\u03b4" }; - const String gamma{ L"\u03b3" }; - const String beta{ L"\u03b2" }; - const String Theta{ L"\u03b8" }; - const String alpha{ L"\u03b1" }; - - String bandName[NUM_BANDS]; - Array validBand[NUM_BANDS]; - Array defaultBand[NUM_BANDS]; - Array extrema[NUM_BANDS]; - int delay[NUM_BANDS]; - Array transformer[NUM_BANDS]; - - static String validBandToString(const Array& band) - { - jassert(band.size() == 2); - return " (" + String(band[0]) + "-" + String(band[1]) + " Hz)"; - } - - HilbertInfo() - { - // Added by sumedh and modified - int del = 7; - Array htcoeff = Array({ -0.433593750000000,0.0,-0.136718750000000,0.0,-0.216796875000000,0.0,-0.638671875000000 }); - - /*delta band*/ - validBand[DELTA] = Array({ 1, 4 }); - bandName[DELTA] = delta + validBandToString(validBand[DELTA]); - defaultBand[DELTA] = Array({ 1, 4 }); - extrema[DELTA] = Array(/* none */); - delay[DELTA] = del; - transformer[DELTA] = htcoeff; - - /*theta band*/ - validBand[THETA] = Array({ 4, 8 }); - bandName[THETA] = Theta + validBandToString(validBand[THETA]); - defaultBand[THETA] = Array({ 4, 8 }); - extrema[THETA] = Array(/* none */); - delay[THETA] = del; - transformer[THETA] = htcoeff; - - /*alpha band*/ - validBand[ALPHA] = Array({ 8, 13 }); - bandName[ALPHA] = alpha + validBandToString(validBand[ALPHA]); - defaultBand[ALPHA] = Array({ 8, 13 }); - extrema[ALPHA] = Array(/* none */); - delay[ALPHA] = del; - transformer[ALPHA] = htcoeff; - - /*beta band*/ - validBand[BETA] = Array({ 12, 30 }); - bandName[BETA] = beta + validBandToString(validBand[BETA]); - defaultBand[BETA] = Array({ 12, 30 }); - extrema[BETA] = Array(/* none */); - delay[BETA] = del; - transformer[BETA] = htcoeff; - - /*lower gamma band*/ - validBand[LOW_GAM] = Array({ 30, 55 }); - bandName[LOW_GAM] = "Lo " + gamma + validBandToString(validBand[LOW_GAM]); - defaultBand[LOW_GAM] = Array({ 30, 55 }); - extrema[LOW_GAM] = Array(/* none */); - delay[LOW_GAM] = del; - transformer[LOW_GAM] = htcoeff; - - /*mid gamma band*/ - validBand[MID_GAM] = Array({ 40, 90 }); - bandName[MID_GAM] = "Mid " + gamma + validBandToString(validBand[MID_GAM]); - //defaultBand[MID_GAM] = Array({ 40, 90 }); - extrema[MID_GAM] = Array(/* none */); - extrema[MID_GAM] = Array({ 64.4559 }); - delay[MID_GAM] = del; - transformer[MID_GAM] = htcoeff; - - /*high gamma band*/ - validBand[HIGH_GAM] = Array({ 60, 200 }); - bandName[HIGH_GAM] = "Hi " + gamma + validBandToString(validBand[HIGH_GAM]); - defaultBand[HIGH_GAM] = Array({ 60, 200 }); - extrema[HIGH_GAM] = Array(/* none */); - //extrema[HIGH_GAM] = Array({ 81.6443, 123.1104, 169.3574 }); - delay[HIGH_GAM] = del; - transformer[HIGH_GAM] = htcoeff; - - /*Ripple band*/ - validBand[RIPPLE] = Array({ 150, 250 }); - bandName[RIPPLE] = "RIP " + validBandToString(validBand[RIPPLE]); - defaultBand[RIPPLE] = Array({ 150, 250 }); - extrema[RIPPLE] = Array(/* none */); - delay[RIPPLE] = del; - transformer[RIPPLE] = htcoeff; - } - }; - - static const HilbertInfo hilbertInfo; // instantiates all the data through the constructor - - struct BandpassfiltInfo - { - String bandName[NUM_BANDS]; - int delay[NUM_BANDS]; - Array transformer[NUM_BANDS]; - BandpassfiltInfo() - { - int del = 42; - /*delta band*/ - delay[DELTA] = del; - transformer[DELTA] = Array({ - 0.0761718750000000,0.0800781250000000,0.0957031250000000,0.123046875000000,0.158203125000000,0.201171875000000,0.253906250000000,0.312500000000000,0.376953125000000,0.443359375000000,0.513671875000000,0.585937500000000,0.654296875000000,0.722656250000000,0.785156250000000,0.841796875000000,0.892578125000000,0.933593750000000,0.964843750000000,0.986328125000000,0.998046875000000,0.998046875000000,0.986328125000000,0.964843750000000,0.933593750000000,0.892578125000000,0.841796875000000,0.785156250000000,0.722656250000000,0.654296875000000,0.585937500000000,0.513671875000000,0.443359375000000,0.376953125000000,0.312500000000000,0.253906250000000,0.201171875000000,0.158203125000000,0.123046875000000,0.0957031250000000,0.0800781250000000,0.0761718750000000 }); - - /*theta band*/ - delay[THETA] = del; - transformer[THETA] = Array({ - 0.0566406250000000,0.0625000000000000,0.0761718750000000,0.0996093750000000,0.132812500000000,0.171875000000000,0.222656250000000,0.277343750000000,0.341796875000000,0.408203125000000,0.480468750000000,0.552734375000000,0.626953125000000,0.697265625000000,0.765625000000000,0.826171875000000,0.880859375000000,0.925781250000000,0.960937500000000,0.986328125000000,0.998046875000000,0.998046875000000,0.986328125000000,0.960937500000000,0.925781250000000,0.880859375000000,0.826171875000000,0.765625000000000,0.697265625000000,0.626953125000000,0.552734375000000,0.480468750000000,0.408203125000000,0.341796875000000,0.277343750000000,0.222656250000000,0.171875000000000,0.132812500000000,0.0996093750000000,0.0761718750000000,0.0625000000000000,0.0566406250000000 }); - - /*alpha band*/ - delay[ALPHA] = del; - transformer[ALPHA] = Array({ - 0.0175781250000000,0.0234375000000000,0.0351562500000000,0.0507812500000000,0.0742187500000000,0.107421875000000,0.148437500000000,0.199218750000000,0.259765625000000,0.326171875000000,0.400390625000000,0.478515625000000,0.558593750000000,0.638671875000000,0.716796875000000,0.789062500000000,0.855468750000000,0.910156250000000,0.953125000000000,0.982421875000000,0.998046875000000,0.998046875000000,0.982421875000000,0.953125000000000,0.910156250000000,0.855468750000000,0.789062500000000,0.716796875000000,0.638671875000000,0.558593750000000,0.478515625000000,0.400390625000000,0.326171875000000,0.259765625000000,0.199218750000000,0.148437500000000,0.107421875000000,0.0742187500000000,0.0507812500000000,0.0351562500000000,0.0234375000000000,0.0175781250000000 }); - - /*beta band*/ - delay[BETA] = del; - transformer[BETA] = Array({ - -0.0566406250000000,-0.0585937500000000,-0.0644531250000000,-0.0722656250000000,-0.0800781250000000,-0.0839843750000000,-0.0781250000000000,-0.0605468750000000,-0.0273437500000000,0.0234375000000000,0.0917968750000000,0.175781250000000,0.277343750000000,0.388671875000000,0.505859375000000,0.623046875000000,0.734375000000000,0.833984375000000,0.912109375000000,0.968750000000000,0.998046875000000,0.998046875000000,0.968750000000000,0.912109375000000,0.833984375000000,0.734375000000000,0.623046875000000,0.505859375000000,0.388671875000000,0.277343750000000,0.175781250000000,0.0917968750000000,0.0234375000000000,-0.0273437500000000,-0.0605468750000000,-0.0781250000000000,-0.0839843750000000,-0.0800781250000000,-0.0722656250000000,-0.0644531250000000,-0.0585937500000000,-0.0566406250000000 }); - - /*low gamma band*/ - delay[LOW_GAM] = del; - transformer[LOW_GAM] = Array({ - 0.0351562500000000,0.0273437500000000,0.0156250000000000,-0.00390625000000000,-0.0371093750000000,-0.0878906250000000,-0.156250000000000,-0.238281250000000,-0.322265625000000,-0.396484375000000,-0.443359375000000,-0.447265625000000,-0.398437500000000,-0.289062500000000,-0.125000000000000,0.0839843750000000,0.318359375000000,0.552734375000000,0.761718750000000,0.916015625000000,0.998046875000000,0.998046875000000,0.916015625000000,0.761718750000000,0.552734375000000,0.318359375000000,0.0839843750000000,-0.125000000000000,-0.289062500000000,-0.398437500000000,-0.447265625000000,-0.443359375000000,-0.396484375000000,-0.322265625000000,-0.238281250000000,-0.156250000000000,-0.0878906250000000,-0.0371093750000000,-0.00390625000000000,0.0156250000000000,0.0273437500000000,0.0351562500000000 }); - - /*mid gamma band*/ - delay[MID_GAM] = del; - transformer[MID_GAM] = Array({ - 0.00195312500000000,0.0,0.00195312500000000,0.0117187500000000,0.0312500000000000,0.0566406250000000,0.0839843750000000,0.0937500000000000,0.0703125000000000,-0.00390625000000000,-0.132812500000000,-0.298828125000000,-0.466796875000000,-0.582031250000000,-0.595703125000000,-0.476562500000000,-0.220703125000000,0.128906250000000,0.501953125000000,0.818359375000000,0.998046875000000,0.998046875000000,0.818359375000000,0.501953125000000,0.128906250000000,-0.220703125000000,-0.476562500000000,-0.595703125000000,-0.582031250000000,-0.466796875000000,-0.298828125000000,-0.132812500000000,-0.00390625000000000,0.0703125000000000,0.0937500000000000,0.0839843750000000,0.0566406250000000,0.0312500000000000,0.0117187500000000,0.00195312500000000,0.0,0.00195312500000000 }); - - /*High gamma band*/ - delay[HIGH_GAM] = del; - transformer[HIGH_GAM] = Array({ - -0.00195312500000000,-0.00781250000000000,-0.0117187500000000,-0.00195312500000000,0.0117187500000000,0.0175781250000000,0.00390625000000000,0.0,0.0390625000000000,0.0917968750000000,0.0820312500000000,-0.0117187500000000,-0.0859375000000000,-0.0371093750000000,0.0468750000000000,-0.0546875000000000,-0.392578125000000,-0.640625000000000,-0.390625000000000,0.341796875000000,0.998046875000000,0.998046875000000,0.341796875000000,-0.390625000000000,-0.640625000000000,-0.392578125000000,-0.0546875000000000,0.0468750000000000,-0.0371093750000000,-0.0859375000000000,-0.0117187500000000,0.0820312500000000,0.0917968750000000,0.0390625000000000,0.0,0.00390625000000000,0.0175781250000000,0.0117187500000000,-0.00195312500000000,-0.0117187500000000,-0.00781250000000000,-0.00195312500000000 }); - - /*Ripple band*/ - delay[RIPPLE] = del; - transformer[RIPPLE] = Array({ - 0.00195312500000000,-0.00195312500000000,0.00195312500000000,0.0195312500000000,0.00976562500000000,-0.0390625000000000,-0.0527343750000000,0.0234375000000000,0.0800781250000000,0.0195312500000000,-0.0234375000000000,0.0292968750000000,-0.0410156250000000,-0.251953125000000,-0.121093750000000,0.449218750000000,0.580078125000000,-0.269531250000000,-0.998046875000000,-0.335937500000000,0.921875000000000,0.921875000000000,-0.335937500000000,-0.998046875000000,-0.269531250000000,0.580078125000000,0.449218750000000,-0.121093750000000,-0.251953125000000,-0.0410156250000000,0.0292968750000000,-0.0234375000000000,0.0195312500000000,0.0800781250000000,0.0234375000000000,-0.0527343750000000,-0.0390625000000000,0.00976562500000000,0.0195312500000000,0.00195312500000000,-0.00195312500000000,0.00195312500000000 }); - - } - }; - - static const BandpassfiltInfo bandpassfiltInfo; // instantiates all the data through the constructor - - struct LowpassfiltInfo - { - String bandName[NUM_BANDS]; - int delay[NUM_BANDS]; - Array transformer[NUM_BANDS]; - - LowpassfiltInfo() - { - int del = 25; - Array lpfcoeff = Array({ 0.0,0.00195312500000000,0.00585937500000000,0.00781250000000000,0.0,-0.0332031250000000,-0.0800781250000000,-0.0917968750000000,0.0,0.238281250000000,0.574218750000000,0.875000000000000,0.998046875000000,0.875000000000000,0.574218750000000,0.238281250000000,0.0,-0.0917968750000000,-0.0800781250000000,-0.0332031250000000,0.0,0.00781250000000000,0.00585937500000000,0.00195312500000000,0.0 }); - /*delta band*/ - delay[DELTA] = del; - transformer[DELTA] = lpfcoeff; - /*theta band*/ - delay[THETA] = del; - transformer[THETA] = lpfcoeff; - /*alpha band*/ - delay[ALPHA] = del; - transformer[ALPHA] = lpfcoeff; - /*beta band*/ - delay[BETA] = del; - transformer[BETA] = lpfcoeff; - /*low gamma band*/ - delay[LOW_GAM] = del; - transformer[LOW_GAM] = lpfcoeff; - /*mid gamma band*/ - delay[MID_GAM] = del; - transformer[MID_GAM] = lpfcoeff; - /*High gamma band*/ - delay[HIGH_GAM] = del; - transformer[HIGH_GAM] = lpfcoeff; - /*Ripple band*/ - delay[RIPPLE] = del; - transformer[RIPPLE] = lpfcoeff; - } - }; - - static const LowpassfiltInfo lowpassfiltInfo; // instantiates all the data through the constructor - } - - namespace Hilbert - { - // exported constants - extern const String* const bandName = hilbertInfo.bandName; - - extern const Array* const validBand = hilbertInfo.validBand; - - extern const Array* const defaultBand = hilbertInfo.defaultBand; - - extern const Array* const extrema = hilbertInfo.extrema; - - extern const int* const delay = hilbertInfo.delay; - - extern const Array* const transformer = hilbertInfo.transformer; - } - - namespace bandpassfilt - { - // exported constants - extern const String* const bandName = bandpassfiltInfo.bandName; - - extern const int* const delay = bandpassfiltInfo.delay; - - extern const Array* const transformer = bandpassfiltInfo.transformer; - } - namespace lowpassfilt - { - // exported constants - extern const String* const bandName = lowpassfiltInfo.bandName; - - extern const int* const delay = lowpassfiltInfo.delay; - - extern const Array* const transformer = lowpassfiltInfo.transformer; - } -} \ No newline at end of file diff --git a/ASICPhaseCalculator/backupSource.zip b/ASICPhaseCalculator/backupSource.zip deleted file mode 100644 index 3a2e273..0000000 Binary files a/ASICPhaseCalculator/backupSource.zip and /dev/null differ diff --git a/ASICPhaseCalculator/Build/.gitignore b/PhaseCalculator/Build/.gitignore similarity index 100% rename from ASICPhaseCalculator/Build/.gitignore rename to PhaseCalculator/Build/.gitignore diff --git a/ASICPhaseCalculator/CMAKE_README.txt b/PhaseCalculator/CMAKE_README.txt similarity index 100% rename from ASICPhaseCalculator/CMAKE_README.txt rename to PhaseCalculator/CMAKE_README.txt diff --git a/ASICPhaseCalculator/CMakeLists.txt b/PhaseCalculator/CMakeLists.txt similarity index 100% rename from ASICPhaseCalculator/CMakeLists.txt rename to PhaseCalculator/CMakeLists.txt diff --git a/ASICPhaseCalculator/Source/ARModeler.h b/PhaseCalculator/Source/ARModeler.h similarity index 100% rename from ASICPhaseCalculator/Source/ARModeler.h rename to PhaseCalculator/Source/ARModeler.h diff --git a/PhaseCalculator/Source/HTransformers.cpp b/PhaseCalculator/Source/HTransformers.cpp new file mode 100644 index 0000000..8d70cce --- /dev/null +++ b/PhaseCalculator/Source/HTransformers.cpp @@ -0,0 +1,146 @@ +/* +------------------------------------------------------------------ + +This file is part of a plugin for the Open Ephys GUI +Copyright (C) 2018 Translational NeuroEngineering Lab + +------------------------------------------------------------------ + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . + +*/ + +#include "HTransformers.h" + +namespace PhaseCalculator +{ + namespace + { + /** Helper to allow assigning to each array one band at a time (within the constructor) **/ + struct HilbertInfo + { + const String gamma{ L"\u03b3" }; + const String beta{ L"\u03b2" }; + const String alphaTheta{ L"\u03b1/\u03b8" }; + + String bandName[NUM_BANDS]; + Array validBand[NUM_BANDS]; + Array defaultBand[NUM_BANDS]; + Array extrema[NUM_BANDS]; + int delay[NUM_BANDS]; + Array transformer[NUM_BANDS]; + + static String validBandToString(const Array& band) + { + jassert(band.size() == 2); + return " (" + String(band[0]) + "-" + String(band[1]) + " Hz)"; + } + + HilbertInfo() + { + validBand[ALPHA_THETA] = Array({ 4, 18 }); + bandName[ALPHA_THETA] = alphaTheta + validBandToString(validBand[ALPHA_THETA]); + defaultBand[ALPHA_THETA] = Array({ 4, 8 }); + extrema[ALPHA_THETA] = Array(/* none */); + delay[ALPHA_THETA] = 9; + // from Matlab: firpm(18, [4 246]/250, [1 1], 'hilbert') + transformer[ALPHA_THETA] = Array({ + -0.28757250783614413, + 0.000027647225074994485, + -0.094611325643268351, + -0.00025887439499763831, + -0.129436276914844, + -0.0001608427426424053, + -0.21315096860055227, + -0.00055322197399797961, + -0.63685698210351149 + }); + + + validBand[BETA] = Array({ 10, 40 }); + bandName[BETA] = beta + validBandToString(validBand[BETA]); + defaultBand[BETA] = Array({ 12, 30 }); + extrema[BETA] = Array({ 21.5848 }); + delay[BETA] = 9; + // from Matlab: firpm(18, [12 30 40 240]/250, [1 1 0.7 0.7], 'hilbert') + transformer[BETA] = Array({ + -0.099949575596234311, + -0.020761484963254036, + -0.080803573080958854, + -0.027365064225587619, + -0.11114477443975329, + -0.025834076852645271, + -0.16664116044989324, + -0.015661948619847599, + -0.45268524264113719 + }); + + + validBand[LOW_GAM] = Array({ 30, 55 }); + bandName[LOW_GAM] = "Lo " + gamma + validBandToString(validBand[LOW_GAM]); + defaultBand[LOW_GAM] = Array({ 30, 55 }); + extrema[LOW_GAM] = Array({ 43.3609 }); + delay[LOW_GAM] = 2; + // from Matlab: firls(4, [30 55]/250, [1 1], 'hilbert') + transformer[LOW_GAM] = Array({ + -1.5933788446351915, + 1.7241339075391682 + }); + + + validBand[MID_GAM] = Array({ 40, 90 }); + bandName[MID_GAM] = "Mid " + gamma + validBandToString(validBand[MID_GAM]); + defaultBand[MID_GAM] = Array({ 40, 90 }); + extrema[MID_GAM] = Array({ 64.4559 }); + delay[MID_GAM] = 2; + // from Matlab: firls(4, [35 90]/250, [1 1], 'hilbert') + transformer[MID_GAM] = Array({ + -0.487176162115735, + -0.069437334858668653 + }); + + + validBand[HIGH_GAM] = Array({ 60, 200 }); + bandName[HIGH_GAM] = "Hi " + gamma + validBandToString(validBand[HIGH_GAM]); + defaultBand[HIGH_GAM] = Array({ 70, 150 }); + extrema[HIGH_GAM] = Array({ 81.6443, 123.1104, 169.3574 }); + delay[HIGH_GAM] = 3; + // from Matlab: firls(6, [60 200]/250, [1 1], 'hilbert') + transformer[HIGH_GAM] = Array({ + -0.10383410506573287, + 0.0040553935691102303, + -0.59258484603659545 + }); + } + }; + + static const HilbertInfo hilbertInfo; // instantiates all the data through the constructor + } + + namespace Hilbert + { + // exported constants + extern const String* const bandName = hilbertInfo.bandName; + + extern const Array* const validBand = hilbertInfo.validBand; + + extern const Array* const defaultBand = hilbertInfo.defaultBand; + + extern const Array* const extrema = hilbertInfo.extrema; + + extern const int* const delay = hilbertInfo.delay; + + extern const Array* const transformer = hilbertInfo.transformer; + } +} \ No newline at end of file diff --git a/ASICPhaseCalculator/Source/HTransformers.h b/PhaseCalculator/Source/HTransformers.h similarity index 68% rename from ASICPhaseCalculator/Source/HTransformers.h rename to PhaseCalculator/Source/HTransformers.h index 4b51751..4525ddd 100644 --- a/ASICPhaseCalculator/Source/HTransformers.h +++ b/PhaseCalculator/Source/HTransformers.h @@ -46,14 +46,11 @@ namespace PhaseCalculator { enum Band { - DELTA = 0, - THETA, - ALPHA, + ALPHA_THETA = 0, BETA, LOW_GAM, MID_GAM, HIGH_GAM, - RIPPLE, NUM_BANDS }; @@ -79,52 +76,6 @@ namespace PhaseCalculator // contain the first delay[band] coefficients; the rest are redundant and can be inferred extern const Array* const transformer; } - - namespace bandpassfilt - { - const int fs = 500; - - // Each pointer below points to an array of length NUM_BANDS. - - extern const String* const bandName; - - // each is a pair (lower limit, upper limit) - extern const Array* const validBand; - - // each is a pair (low cut, high cut) - extern const Array* const defaultBand; - - extern const Array* const extrema; - - // samples of group delay (= order of filter / 2) - extern const int* const delay; - - // contain the first delay[band] coefficients; the rest are redundant and can be inferred - extern const Array* const transformer; - } - - namespace lowpassfilt - { - const int fs = 500; - - // Each pointer below points to an array of length NUM_BANDS. - - extern const String* const bandName; - - // each is a pair (lower limit, upper limit) - extern const Array* const validBand; - - // each is a pair (low cut, high cut) - extern const Array* const defaultBand; - - extern const Array* const extrema; - - // samples of group delay (= order of filter / 2) - extern const int* const delay; - - // contain the first delay[band] coefficients; the rest are redundant and can be inferred - extern const Array* const transformer; - } } #endif // H_TRANSFORMERS_H_INCLUDED \ No newline at end of file diff --git a/ASICPhaseCalculator/Source/OpenEphysLib.cpp b/PhaseCalculator/Source/OpenEphysLib.cpp similarity index 97% rename from ASICPhaseCalculator/Source/OpenEphysLib.cpp rename to PhaseCalculator/Source/OpenEphysLib.cpp index a275cc0..3c02236 100644 --- a/ASICPhaseCalculator/Source/OpenEphysLib.cpp +++ b/PhaseCalculator/Source/OpenEphysLib.cpp @@ -48,7 +48,7 @@ extern "C" EXPORT int getPluginInfo(int index, Plugin::PluginInfo* info) { case 0: info->type = Plugin::PLUGIN_TYPE_PROCESSOR; - info->processor.name = "ASIC"; + info->processor.name = "Phase Calculator"; info->processor.type = Plugin::FilterProcessor; info->processor.creator = &(Plugin::createProcessor); break; diff --git a/ASICPhaseCalculator/Source/PhaseCalculator.cpp b/PhaseCalculator/Source/PhaseCalculator.cpp similarity index 72% rename from ASICPhaseCalculator/Source/PhaseCalculator.cpp rename to PhaseCalculator/Source/PhaseCalculator.cpp index 79c570a..c4fcca7 100644 --- a/ASICPhaseCalculator/Source/PhaseCalculator.cpp +++ b/PhaseCalculator/Source/PhaseCalculator.cpp @@ -37,6 +37,9 @@ namespace PhaseCalculator // (meant to be larger than actual minimum floating-point eps) static const float passbandEps = 0.01F; + // priority of the AR model calculating thread (0 = lowest, 10 = highest) + static const int arPriority = 3; + // "glitch limit" (how long of a segment is allowed to be unwrapped or smoothed, in samples) static const int glitchLimit = 200; @@ -46,6 +49,7 @@ namespace PhaseCalculator static const int visMinDelayMs = 675; static const int visMaxDelayMs = 1000; + /*** ReverseStack ***/ ReverseStack::ReverseStack(int size) : freeSpace (size) @@ -134,34 +138,36 @@ namespace PhaseCalculator void ActiveChannelInfo::update() { const Node& p = chanInfo.owner; - + int arOrder = p.getAROrder(); float highCut = p.getHighCut(); float lowCut = p.getLowCut(); Band band = p.getBand(); - int bpforder = p.getBpfOrder(); + + // update length of history based on sample rate + // the history buffer should have enough samples to calculate phases for the viusalizer + // with the proper Hilbert transform length AND train an AR model of the requested order, + // using at least 1 second of data + int newHistorySize = chanInfo.dsFactor * jmax( + visHilbertLengthMs * Hilbert::fs / 1000, + arOrder + 1, + 1 * Hilbert::fs); + + history.resetAndResize(newHistorySize); + // set filter parameters - // set filter parameters - for (auto filt : { &filter, &reverseFilter }) - { - filt->setup( - bpforder, // order - chanInfo.sampleRate, // sample rate - (highCut + lowCut) / 2, // center frequency - highCut - lowCut); // bandwidth - } - // added by sumedh - int lpforder = p.getLpfOrder(); - - - // low pass at 250 Hz - for (auto lpffilt1 : { &filterlpf }) - { - lpffilt1->setup(lpforder, (double)chanInfo.sampleRate, (double)(.8*(chanInfo.sampleRate/2)/ (chanInfo.sampleRate / 4000)),0.05); - } + for (auto filt : { &filter, &reverseFilter }) + { + filt->setup( + 2, // order + chanInfo.sampleRate, // sample rate + (highCut + lowCut) / 2, // center frequency + highCut - lowCut); // bandwidth + } + + arModeler.setParams(arOrder, newHistorySize, chanInfo.dsFactor); htState.resize(Hilbert::delay[band] * 2 + 1); - bpfState.resize(bandpassfilt::delay[band]); - lpfState.resize(lowpassfilt::delay[band]); + // visualization stuff hilbertLengthMultiplier = Hilbert::fs * chanInfo.dsFactor / 1000; visHilbertBuffer.resize(visHilbertLengthMs * hilbertLengthMultiplier); @@ -173,11 +179,8 @@ namespace PhaseCalculator { history.reset(); filter.reset(); - //added by sumedh - filterlpf.reset(); + arModeler.reset(); FloatVectorOperations::clear(htState.begin(), htState.size()); - FloatVectorOperations::clear(bpfState.begin(), bpfState.size()); - FloatVectorOperations::clear(lpfState.begin(), lpfState.size()); interpCountdown = 0; lastComputedPhase = 0; lastComputedMag = 0; @@ -247,17 +250,16 @@ namespace PhaseCalculator /**** phase calculator node ****/ Node::Node() - : GenericProcessor("ASIC Phase Calculator") + : GenericProcessor("Phase Calculator") , Thread("AR Modeler") - // added by sumedh - , lpforder(8) - , bpforder(2) + , calcInterval(50) + , arOrder(20) , outputMode(PH) , visEventChannel(-1) , visContinuousChannel(-1) { setProcessorType(PROCESSOR_TYPE_FILTER); - setBand(THETA, true); + setBand(ALPHA_THETA, true); } Node::~Node() {} @@ -309,10 +311,27 @@ namespace PhaseCalculator void Node::setParameter(int parameterIndex, float newValue) { switch (parameterIndex) { + case RECALC_INTERVAL: + calcInterval = int(newValue); + break; + + case AR_ORDER: + arOrder = int(newValue); + updateActiveChannels(); + break; case BAND: setBand(Band(int(newValue))); break; + + case LOWCUT: + setLowCut(newValue); + break; + + case HIGHCUT: + setHighCut(newValue); + break; + case OUTPUT_MODE: { OutputMode oldMode = outputMode; @@ -369,147 +388,168 @@ namespace PhaseCalculator { continue; } - int downsampleA = 4000; - int downsampleB = 1000; - /*Get a pointer to the data wpInput*/ - float* const wpInput = buffer.getWritePointer(chan); - /*active channel info -> acInfo*/ - /*Step I: IIR Low pass filter over entire data */ - acInfo->filterlpf.process(nSamples, &wpInput); - std::vector dslpf; - int preDSample = int(floor(float(acInfo->chanInfo.sampleRate) / downsampleA)); - for (int i = 0; i < nSamples; i = i + preDSample) { - dslpf.push_back(wpInput[i]); - } - /*Step II: FIR low pass filter from MATLAB*/ - /*lpf again with the coeff from MATLAB*/ - std::vector lpfData; - for (int lpfIndex = 0; lpfIndex < dslpf.size(); ++lpfIndex) { - lpfData.push_back(lpfFilterSamp(dslpf[lpfIndex], band, acInfo->lpfState)); - } - /*Step III: Downsampling from 40000K to 1K*/ - - int new_sr_size = downsampleA / downsampleB;// int(floor(float(acInfo->chanInfo.sampleRate) / float(new_sampling_rate))); - std::vector dsLpfData; - for (int i = 0; i < dslpf.size(); i=i+new_sr_size) { - dsLpfData.push_back(lpfData[i]); - } - /*Step IV: FIR Band pass filtering on the downsampled data*/ - int newdatasize = dsLpfData.size(); - std::vector dsBPFdata; - for (int bpfIndex = 0; bpfIndex < dsLpfData.size(); ++bpfIndex) { - dsBPFdata.push_back(bpfFilterSamp(dsLpfData[bpfIndex], band, acInfo->bpfState)); - } - /*Step V: calculate Hilbert transform*/ - int htOutputSamps = dsBPFdata.size(); - if (htOutput.size() < htOutputSamps) - { - htOutput.resize(htOutputSamps); - } - htOutput.resize(htOutputSamps); - /*Changes Both real and imaginary part are calculated from the same*/ - //float* wpOut = buffer.getWritePointer(chan); - for (int hilbertIndex = 0; hilbertIndex < htOutputSamps; ++hilbertIndex) - { - double samp = htFilterSamp(dsBPFdata[hilbertIndex], band, acInfo->htState); - double rc = dsBPFdata[hilbertIndex]; - double ic = htScaleFactor * samp; - htOutput.set(hilbertIndex, std::complex(rc, ic)); - //wpOut[hilbertIndex] = (float)rc; - } - /*Step VI: Interpolation: No new changes proposed*/ - int stride = int(floor(float(acInfo->chanInfo.sampleRate) / float(downsampleB))) + 2; - float* wpOut = buffer.getWritePointer(chan); - float* wpOut2; - if (outputMode == PH_AND_MAG) - { - // second output channel - int outChan2 = getNumInputs() + ac; - jassert(outChan2 < buffer.getNumChannels()); - wpOut2 = buffer.getWritePointer(outChan2); - } - - double nextComputedPhase, phaseStep; - double nextComputedMag, magStep; - bool needPhase = outputMode != MAG; - bool needMag = outputMode != PH; - - if (needPhase) - { - nextComputedPhase = LAA(htOutput[0]);//std::arg(htOutput[0]); - phaseStep = circDist(nextComputedPhase, acInfo->lastComputedPhase, Dsp::doublePi) / stride; - } - if (needMag) - { - nextComputedMag = std::abs(htOutput[0]); - magStep = (nextComputedMag - acInfo->lastComputedMag) / stride; - } - - for (int i = 0, frame = 0; i < nSamples; ++i, --acInfo->interpCountdown) - { - if (acInfo->interpCountdown == 0) - { - // update interpolation frame - ++frame; - acInfo->interpCountdown = stride; - - if (needPhase) - { - acInfo->lastComputedPhase = nextComputedPhase; - nextComputedPhase = LAA(htOutput[frame]);//std::arg(htOutput[frame]); - phaseStep = circDist(nextComputedPhase, acInfo->lastComputedPhase, Dsp::doublePi) / stride; - } - if (needMag) - { - acInfo->lastComputedMag = nextComputedMag; - nextComputedMag = std::abs(htOutput[frame]); - magStep = (nextComputedMag - acInfo->lastComputedMag) / stride; - } - } - - double thisPhase, thisMag; - if (needPhase) - { - thisPhase = circDist(nextComputedPhase, phaseStep * acInfo->interpCountdown, Dsp::doublePi); - } - if (needMag) - { - thisMag = nextComputedMag - magStep * acInfo->interpCountdown; - } - - switch (outputMode) - { - case MAG: - wpOut[i] = float(thisMag); - break; - - case PH_AND_MAG: - wpOut2[i] = float(thisMag); - // fall through - case PH: - // output in degrees - wpOut[i] = float(thisPhase * (180.0 / Dsp::doublePi)); - break; - - case IM: - wpOut[i] = float(thisMag * std::sin(thisPhase)); - break; - } - } - // unwrapping / smoothing - if (outputMode == PH || outputMode == PH_AND_MAG) - { - unwrapBuffer(wpOut, nSamples, acInfo->lastPhase); - smoothBuffer(wpOut, nSamples, acInfo->lastPhase); - acInfo->lastPhase = wpOut[nSamples - 1]; - } - // if this is the monitored channel for events, check whether we can add a new phase - /*if (hasCanvas && chan == visContinuousChannel && acInfo->history.isFull()) - { - calcVisPhases(acInfo, getTimestamp(chan) + getNumSamples(chan)); - }*/ + // filter the data + float* const wpIn = buffer.getWritePointer(chan); + acInfo->filter.process(nSamples, &wpIn); + + // enqueue as much new data as can fit into history + acInfo->history.enqueue(wpIn, nSamples); + + // calc phase and write out (only if AR model has been calculated) + if (acInfo->history.isFull() && acInfo->arModeler.hasBeenFit()) + { + // read current AR parameters safely (uses lock internally) + acInfo->arModeler.getModel(localARParams); + + // use AR model to fill predSamps (which is downsampled) based on past data. + int htDelay = Hilbert::delay[band]; + int stride = acInfo->chanInfo.dsFactor; + + double* pPredSamps = predSamps.getRawDataPointer(); + const double* pLocalParam = localARParams.getRawDataPointer(); + arPredict(acInfo->history, acInfo->interpCountdown, pPredSamps, pLocalParam, + htDelay + 1, stride, arOrder); + + // identify indices of current buffer to execute HT + htInds.clearQuick(); + for (int i = acInfo->interpCountdown; i < nSamples; i += stride) + { + htInds.add(i); + } + + int htOutputSamps = htInds.size() + 1; + if (htOutput.size() < htOutputSamps) + { + htOutput.resize(htOutputSamps); + } + + // execute tranformer on current buffer + int kOut = -htDelay; + for (int kIn = 0; kIn < htInds.size(); ++kIn, ++kOut) + { + double samp = htFilterSamp(wpIn[htInds[kIn]], band, acInfo->htState); + if (kOut >= 0) + { + double rc = wpIn[htInds[kOut]]; + double ic = htScaleFactor * samp; + htOutput.set(kOut, std::complex(rc, ic)); + } + } + + // copy state to transform prediction without changing the end-of-buffer state + htTempState = acInfo->htState; + + // execute transformer on prediction + for (int i = 0; i <= htDelay; ++i, ++kOut) + { + double samp = htFilterSamp(predSamps[i], band, htTempState); + if (kOut >= 0) + { + double rc = i == htDelay ? predSamps[0] : wpIn[htInds[kOut]]; + double ic = htScaleFactor * samp; + htOutput.set(kOut, std::complex(rc, ic)); + } + } + + // output with upsampling (interpolation) + float* wpOut = buffer.getWritePointer(chan); + float* wpOut2; + if (outputMode == PH_AND_MAG) + { + // second output channel + int outChan2 = getNumInputs() + ac; + jassert(outChan2 < buffer.getNumChannels()); + wpOut2 = buffer.getWritePointer(outChan2); + } + + double nextComputedPhase, phaseStep; + double nextComputedMag, magStep; + bool needPhase = outputMode != MAG; + bool needMag = outputMode != PH; + + if (needPhase) + { + nextComputedPhase = std::arg(htOutput[0]); + phaseStep = circDist(nextComputedPhase, acInfo->lastComputedPhase, Dsp::doublePi) / stride; + } + if (needMag) + { + nextComputedMag = std::abs(htOutput[0]); + magStep = (nextComputedMag - acInfo->lastComputedMag) / stride; + } + + for (int i = 0, frame = 0; i < nSamples; ++i, --acInfo->interpCountdown) + { + if (acInfo->interpCountdown == 0) + { + // update interpolation frame + ++frame; + acInfo->interpCountdown = stride; + + if (needPhase) + { + acInfo->lastComputedPhase = nextComputedPhase; + nextComputedPhase = std::arg(htOutput[frame]); + phaseStep = circDist(nextComputedPhase, acInfo->lastComputedPhase, Dsp::doublePi) / stride; + } + if (needMag) + { + acInfo->lastComputedMag = nextComputedMag; + nextComputedMag = std::abs(htOutput[frame]); + magStep = (nextComputedMag - acInfo->lastComputedMag) / stride; + } + } + + double thisPhase, thisMag; + if (needPhase) + { + thisPhase = circDist(nextComputedPhase, phaseStep * acInfo->interpCountdown, Dsp::doublePi); + } + if (needMag) + { + thisMag = nextComputedMag - magStep * acInfo->interpCountdown; + } + + switch (outputMode) + { + case MAG: + wpOut[i] = float(thisMag); + break; + + case PH_AND_MAG: + wpOut2[i] = float(thisMag); + // fall through + case PH: + // output in degrees + wpOut[i] = float(thisPhase * (180.0 / Dsp::doublePi)); + break; + + case IM: + wpOut[i] = float(thisMag * std::sin(thisPhase)); + break; + } + } + + // unwrapping / smoothing + if (outputMode == PH || outputMode == PH_AND_MAG) + { + unwrapBuffer(wpOut, nSamples, acInfo->lastPhase); + smoothBuffer(wpOut, nSamples, acInfo->lastPhase); + acInfo->lastPhase = wpOut[nSamples - 1]; + } + } + else // fifo not full or AR model not ready + { + // just output zeros + buffer.clear(chan, 0, nSamples); + } + // if this is the monitored channel for events, check whether we can add a new phase + if (hasCanvas && chan == visContinuousChannel && acInfo->history.isFull()) + { + calcVisPhases(acInfo, getTimestamp(chan) + getNumSamples(chan)); + } } } @@ -518,7 +558,7 @@ namespace PhaseCalculator { if (isEnabled) { - //startThread(arPriority); + startThread(arPriority); // have to manually enable editor, I guess... Editor* editor = static_cast(getEditor()); @@ -576,6 +616,34 @@ namespace PhaseCalculator Array reverseData; reverseData.resize(maxHistoryLength); + + uint32 startTime, endTime; + while (!threadShouldExit()) + { + startTime = Time::getMillisecondCounter(); + + for (auto acInfo : activeChans) + { + if (!acInfo->history.isFull()) + { + continue; + } + + // unwrap reversed history and add to temporary data array + double* dataPtr = reverseData.getRawDataPointer(); + acInfo->history.unwrapAndCopy(dataPtr, true); + + // calculate parameters + acInfo->arModeler.fitModel(reverseData); + } + + endTime = Time::getMillisecondCounter(); + int remainingInterval = calcInterval - (endTime - startTime); + if (remainingInterval >= 10) // avoid WaitForSingleObject + { + sleep(remainingInterval); + } + } } void Node::updateSettings() @@ -657,16 +725,10 @@ namespace PhaseCalculator return int(getProcessorFullId(sourceNodeId, subProcessorIdx)); } - // Added by Sumedh to access private variables - int Node::getLpfOrder() const - { - return lpforder; - } - - int Node::getBpfOrder() const - { - return bpforder; - } + int Node::getAROrder() const + { + return arOrder; + } float Node::getHighCut() const { @@ -779,10 +841,71 @@ namespace PhaseCalculator const Array& defaultBand = Hilbert::defaultBand[band]; lowCut = defaultBand[0]; highCut = defaultBand[1]; + + auto editor = static_cast(getEditor()); + if (editor) + { + editor->refreshLowCut(); + editor->refreshHighCut(); + } + + updateScaleFactor(); + updateActiveChannels(); + } + + void Node::setLowCut(float newLowCut) + { + if (newLowCut == lowCut) { return; } + + auto editor = static_cast(getEditor()); + const Array& validBand = Hilbert::validBand[band]; + + if (newLowCut < validBand[0] || newLowCut >= validBand[1]) + { + // invalid; don't set parameter and reset editor + editor->refreshLowCut(); + CoreServices::sendStatusMessage("Low cut outside valid band of selected filter."); + return; + } + + lowCut = newLowCut; + if (lowCut >= highCut) + { + // push highCut up + highCut = jmin(lowCut + passbandEps, validBand[1]); + editor->refreshHighCut(); + } + updateScaleFactor(); updateActiveChannels(); } + void Node::setHighCut(float newHighCut) + { + if (newHighCut == highCut) { return; } + + auto editor = static_cast(getEditor()); + const Array& validBand = Hilbert::validBand[band]; + + if (newHighCut <= validBand[0] || newHighCut > validBand[1]) + { + // invalid; don't set parameter and reset editor + editor->refreshHighCut(); + CoreServices::sendStatusMessage("High cut outside valid band of selected filter."); + return; + } + + highCut = newHighCut; + if (highCut <= lowCut) + { + // push lowCut down + lowCut = jmax(highCut - passbandEps, validBand[0]); + editor->refreshLowCut(); + } + + updateScaleFactor(); + updateActiveChannels(); + } void Node::setVisContChan(int newChan) { @@ -815,7 +938,7 @@ namespace PhaseCalculator void Node::updateScaleFactor() { - htScaleFactor = getScaleFactor(THETA, 4.0, 8.0); + htScaleFactor = getScaleFactor(band, lowCut, highCut); } void Node::unwrapBuffer(float* wp, int nSamples, float lastPhase) @@ -1064,7 +1187,7 @@ namespace PhaseCalculator acInfo->reverseFilter.reset(); acInfo->reverseFilter.process(hilbertLength, &wpHilbert); - //un-reverse values + // un-reverse values acInfo->visHilbertBuffer.reverseReal(hilbertLength); // Hilbert transform! @@ -1143,6 +1266,29 @@ namespace PhaseCalculator channelInfo.getUnchecked(chan)->deactivate(); } + void Node::arPredict(const ReverseStack& history, int interpCountdown, double* prediction, + const double* params, int samps, int stride, int order) + { + const double* rpHistory = history.begin(); + int histSize = history.size(); + int histStart = history.getHeadOffset() + stride - interpCountdown; + + // s = index to write output + for (int s = 0; s < samps; ++s) + { + prediction[s] = 0; + + // p = which AR param we are on + for (int p = 0; p < order; ++p) + { + double pastSamp = p < s + ? prediction[s - 1 - p] + : rpHistory[(histStart + (p - s) * stride) % histSize]; + + prediction[s] -= params[p] * pastSamp; + } + } + } double Node::getScaleFactor(Band band, double lowCut, double highCut) { @@ -1212,117 +1358,6 @@ namespace PhaseCalculator std::memmove(state_p, state_p + 1, order * sizeof(double)); return sampOut; } - //Added by Sumedh - /* - input: gets the current value to be filtered - Band: The band that needs to be filtered, actually this is identified with Matlab function - state: the bandpass filter states needs to keep track of the data - Steps - 1) move the memory and add the input at the 0 - 2) std::multiplies state * coefficient - 3) sum the entire states - 4) return the sum - */ - double Node::lpfFilterSamp(double input, Band band, Array& state) - { - double* state_p = state.getRawDataPointer(); - int nCoefs = lowpassfilt::delay[band]; - /*Considering the entire coefficient*/ - int order = nCoefs; // considering the current even order - //to the right - //double temp = state_p[order - 1]; //remember last element - for (int i = order - 1; i >= 0; i--) - { - state_p[i + 1] = state_p[i]; //move all element to the right except last one - } - state_p[0] = 0.0; //assign remembered value to first element - //std::memmove(state_p + 1, state_p, order * sizeof(double)); - state_p[0] = input; - double retValue = 0.0; - const double* transf = lowpassfilt::transformer[band].begin(); - for (int kCoef = 0; kCoef < nCoefs; ++kCoef) - { - retValue = retValue + state_p[kCoef] * transf[kCoef]; - } - - return retValue; - } - /* - input: gets the current value to be filtered - Band: The band that needs to be filtered, actually this is identified with Matlab function - state: the bandpass filter states needs to keep track of the data - Steps - 1) move the memory and add the input at the 0 - 2) std::multiplies state * coefficient - 3) sum the entire states - 4) return the sum - */ - double Node::bpfFilterSamp(double input, Band band, Array& state) - { - double* state_p = state.getRawDataPointer(); - int nCoefs = bandpassfilt::delay[band]; - double retValue = 0.0; - int order = nCoefs; // considering the current even order - state_p[order - 1] = 0.0; - for (int i = order - 1; i >= 0; i--) - { - state_p[i + 1] = state_p[i]; //move all element to the right except last one - } - state_p[0] = 0.0; //assign remembered value to first element - //std::memmove(state_p + 1, state_p, order-1 * sizeof(double)); - state_p[0] = input; - const double* transf = bandpassfilt::transformer[band].begin(); - for (int kCoef = 0; kCoef < nCoefs; ++kCoef) - { - retValue = retValue + state_p[kCoef] * transf[kCoef]; - } - return retValue; - } - - double Node::LAA(std::complex c) - { - // Hold phase before quantization - double q; - - // Determine phase based on octant. See LAA alg details. - if (std::abs(c.real()) >= std::abs(c.imag())) - { - if (c.real() >= 0) - { - q = (1. / 8.) * (c.imag() / c.real()); // octant 1 and 8 - - } - else - { - if (c.imag() >= 0) - { - q = .5 + (1. / 8.) * (c.imag() / c.real()); // octant 4 - } - else - { - q = -.5 + (1. / 8.) * (c.imag() / c.real()); // octant 5 - } - } - } - else - { - if (c.imag() >= 0) - { - q = 0.25 - (1. / 8.) * (c.real() / c.imag()); // octant 2 and 3 - } - else - { - q = -.25 - (1. / 8.) * (c.real() / c.imag()); // octant 6 and 7 - } - } - - // Do quantization on phase (based on bit percision). We are testing 8 bits. - double lsb = 1. / pow(2., 8.); - double ans = floor(q / lsb) * lsb * 3.14 / 0.5; - if (std::isnan(ans)) - ans = 0.0; - return ans; - } } diff --git a/ASICPhaseCalculator/Source/PhaseCalculator.h b/PhaseCalculator/Source/PhaseCalculator.h similarity index 90% rename from ASICPhaseCalculator/Source/PhaseCalculator.h rename to PhaseCalculator/Source/PhaseCalculator.h index 1445e10..24e17fc 100644 --- a/ASICPhaseCalculator/Source/PhaseCalculator.h +++ b/PhaseCalculator/Source/PhaseCalculator.h @@ -62,13 +62,11 @@ namespace PhaseCalculator enum Param { RECALC_INTERVAL, + AR_ORDER, BAND, LOWCUT, HIGHCUT, OUTPUT_MODE, - // Added by Sumedh - LPF_ORDER, - BPF_ORDER, VIS_E_CHAN, VIS_C_CHAN }; @@ -110,18 +108,10 @@ namespace PhaseCalculator // reset to perform after end of acquisition or update void reset(); - // Added by Sumedh to create low pass filter - // filter design copied from FilterNode - using LowPassFilter = Dsp::SimpleFilter - , // order - 1, // number of channels - Dsp::DirectFormI>; // realization - LowPassFilter filterlpf; // filter design copied from FilterNode using BandpassFilter = Dsp::SimpleFilter - , // order 1, // number of channels Dsp::DirectFormII>; // realization @@ -130,11 +120,9 @@ namespace PhaseCalculator BandpassFilter filter; + ARModeler arModeler; Array htState; - Array bpfState; - Array lpfState; - Array lpfStateMain; // number of samples until a new non-interpolated output. e.g. if this // equals 1 after a buffer is processed, then there is one interpolated @@ -232,13 +220,11 @@ namespace PhaseCalculator int getFullSourceId(int chan); // getters + int getAROrder() const; float getHighCut() const; float getLowCut() const; Band getBand() const; - - // Added by Sumedh to access the private variables - int getLpfOrder() const; - int getBpfOrder() const; + // reads from the visPhaseBuffer if it can acquire a TryLock. returns true if successful. bool tryToReadVisPhases(std::queue& other); @@ -267,6 +253,11 @@ namespace PhaseCalculator // Resets lowCut and highCut to defaults for the current band void resetCutsToDefaults(); + // Sets lowCut (which in turn influences highCut) + void setLowCut(float newLowCut); + + // Sets highCut (which in turn influences lowCut) + void setHighCut(float newHighCut); // Sets visContinuousChannel and updates the visualization filter void setVisContChan(int newChan); @@ -324,6 +315,22 @@ namespace PhaseCalculator void deactivateInputChannel(int chan); + // ---- static utility methods ---- + + /* + * arPredict: use autoregressive model of order to predict future data. + * + * lastSample points to the most recent sample of past data that will be used to + * compute phase, and there must be at least stride * (order - 1) samples + * preceding it in order to do the AR prediction. + * + * Input params is an array of coefficients of an AR model of length 'order'. + * + * Writes samps future data values to prediction. + */ + static void arPredict(const ReverseStack& history, int interpCountdown, double* prediction, + const double* params, int samps, int stride, int order); + // Get the htScaleFactor for the given band's Hilbert transformer, // over the range from lowCut and highCut. This is the reciprocal of the geometric // mean (i.e. mean in decibels) of the maximum and minimum magnitude responses over the range. @@ -331,17 +338,14 @@ namespace PhaseCalculator // Execute the hilbert transformer on one sample and update the state. static double htFilterSamp(double input, Band band, Array& state); - //Added by Sumedh - static double bpfFilterSamp(double input, Band band, Array& state); - static double lpfFilterSamp(double input, Band band, Array& state); - static double lpfFilter(double input, Band band, Array& state); - // ---- customizable parameters ------ - // Added by Sumedh - // Use LAA to return phase angle in degrees - double LAA(std::complex); - int lpforder; - int bpforder; + // ---- customizable parameters ------ + + // time to wait between AR model recalculations in ms + int calcInterval; + + // order of the AR model + int arOrder; OutputMode outputMode; diff --git a/ASICPhaseCalculator/Source/PhaseCalculatorCanvas.cpp b/PhaseCalculator/Source/PhaseCalculatorCanvas.cpp similarity index 100% rename from ASICPhaseCalculator/Source/PhaseCalculatorCanvas.cpp rename to PhaseCalculator/Source/PhaseCalculatorCanvas.cpp diff --git a/ASICPhaseCalculator/Source/PhaseCalculatorCanvas.h b/PhaseCalculator/Source/PhaseCalculatorCanvas.h similarity index 100% rename from ASICPhaseCalculator/Source/PhaseCalculatorCanvas.h rename to PhaseCalculator/Source/PhaseCalculatorCanvas.h diff --git a/ASICPhaseCalculator/Source/PhaseCalculatorEditor.cpp b/PhaseCalculator/Source/PhaseCalculatorEditor.cpp similarity index 64% rename from ASICPhaseCalculator/Source/PhaseCalculatorEditor.cpp rename to PhaseCalculator/Source/PhaseCalculatorEditor.cpp index 7c47462..7a1ecc5 100644 --- a/ASICPhaseCalculator/Source/PhaseCalculatorEditor.cpp +++ b/PhaseCalculator/Source/PhaseCalculatorEditor.cpp @@ -28,14 +28,14 @@ along with this program. If not, see . #include // abs namespace PhaseCalculator -{ //added by sumedh +{ Editor::Editor(Node* parentNode, bool useDefaultParameterEditors) - : VisualizerEditor(parentNode, 250, useDefaultParameterEditors) + : VisualizerEditor(parentNode, 220, useDefaultParameterEditors) , extraChanManager(parentNode) , prevExtraChans(0) { tabText = "Event Phase Plot"; - int filterWidth = 5; + int filterWidth = 120; // make the canvas now, so that restoring its parameters always works. canvas = new Canvas(parentNode); @@ -57,8 +57,88 @@ namespace PhaseCalculator bandBox->addListener(this); addAndMakeVisible(bandBox); + lowCutLabel = new Label("lowCutL", "Low:"); + lowCutLabel->setBounds(5, 70, 50, 20); + lowCutLabel->setFont({ "Small Text", 12, Font::plain }); + lowCutLabel->setColour(Label::textColourId, Colours::darkgrey); + addAndMakeVisible(lowCutLabel); + + lowCutEditable = new Label("lowCutE"); + lowCutEditable->setEditable(true); + lowCutEditable->addListener(this); + lowCutEditable->setBounds(50, 70, 35, 18); + lowCutEditable->setText(String(parentNode->lowCut), dontSendNotification); + lowCutEditable->setColour(Label::backgroundColourId, Colours::grey); + lowCutEditable->setColour(Label::textColourId, Colours::white); + addAndMakeVisible(lowCutEditable); + + lowCutUnit = new Label("lowCutU", "Hz"); + lowCutUnit->setBounds(85, 70, 25, 18); + lowCutUnit->setFont({ "Small Text", 12, Font::plain }); + lowCutUnit->setColour(Label::textColourId, Colours::darkgrey); + addAndMakeVisible(lowCutUnit); + + highCutLabel = new Label("highCutL", "High:"); + highCutLabel->setBounds(5, 100, 50, 20); + highCutLabel->setFont({ "Small Text", 12, Font::plain }); + highCutLabel->setColour(Label::textColourId, Colours::darkgrey); + addAndMakeVisible(highCutLabel); + + highCutEditable = new Label("highCutE"); + highCutEditable->setEditable(true); + highCutEditable->addListener(this); + highCutEditable->setBounds(50, 100, 35, 18); + highCutEditable->setText(String(parentNode->highCut), dontSendNotification); + highCutEditable->setColour(Label::backgroundColourId, Colours::grey); + highCutEditable->setColour(Label::textColourId, Colours::white); + addAndMakeVisible(highCutEditable); + + highCutUnit = new Label("highCutU", "Hz"); + highCutUnit->setBounds(85, 100, 25, 18); + highCutUnit->setFont({ "Small Text", 12, Font::plain }); + highCutUnit->setColour(Label::textColourId, Colours::darkgrey); + addAndMakeVisible(highCutUnit); + + recalcIntervalLabel = new Label("recalcL", "AR Refresh:"); + recalcIntervalLabel->setBounds(filterWidth, 25, 100, 20); + recalcIntervalLabel->setFont({ "Small Text", 12, Font::plain }); + recalcIntervalLabel->setColour(Label::textColourId, Colours::darkgrey); + addAndMakeVisible(recalcIntervalLabel); + + recalcIntervalEditable = new Label("recalcE"); + recalcIntervalEditable->setEditable(true); + recalcIntervalEditable->addListener(this); + recalcIntervalEditable->setBounds(filterWidth + 5, 44, 55, 18); + recalcIntervalEditable->setColour(Label::backgroundColourId, Colours::grey); + recalcIntervalEditable->setColour(Label::textColourId, Colours::white); + recalcIntervalEditable->setText(String(parentNode->calcInterval), dontSendNotification); + recalcIntervalEditable->setTooltip(recalcIntervalTooltip); + addAndMakeVisible(recalcIntervalEditable); + + recalcIntervalUnit = new Label("recalcU", "ms"); + recalcIntervalUnit->setBounds(filterWidth + 60, 47, 25, 15); + recalcIntervalUnit->setFont({ "Small Text", 12, Font::plain }); + recalcIntervalUnit->setColour(Label::textColourId, Colours::darkgrey); + addAndMakeVisible(recalcIntervalUnit); + + arOrderLabel = new Label("arOrderL", "Order:"); + arOrderLabel->setBounds(filterWidth, 65, 60, 20); + arOrderLabel->setFont({ "Small Text", 12, Font::plain }); + arOrderLabel->setColour(Label::textColourId, Colours::darkgrey); + addAndMakeVisible(arOrderLabel); + + arOrderEditable = new Label("arOrderE"); + arOrderEditable->setEditable(true); + arOrderEditable->addListener(this); + arOrderEditable->setBounds(filterWidth + 55, 66, 25, 18); + arOrderEditable->setColour(Label::backgroundColourId, Colours::grey); + arOrderEditable->setColour(Label::textColourId, Colours::white); + arOrderEditable->setText(String(parentNode->arOrder), sendNotificationAsync); + arOrderEditable->setTooltip(arOrderTooltip); + addAndMakeVisible(arOrderEditable); + outputModeLabel = new Label("outputModeL", "Output:"); - outputModeLabel->setBounds(filterWidth, 72, 70, 20); + outputModeLabel->setBounds(filterWidth, 87, 70, 20); outputModeLabel->setFont({ "Small Text", 12, Font::plain }); outputModeLabel->setColour(Label::textColourId, Colours::darkgrey); addAndMakeVisible(outputModeLabel); @@ -70,7 +150,7 @@ namespace PhaseCalculator outputModeBox->addItem("IMAG", IM); outputModeBox->setSelectedId(parentNode->outputMode); outputModeBox->setTooltip(outputModeTooltip); - outputModeBox->setBounds(filterWidth + 5, 92, 76, 19); + outputModeBox->setBounds(filterWidth + 5, 105, 76, 19); outputModeBox->addListener(this); addAndMakeVisible(outputModeBox); @@ -97,6 +177,47 @@ namespace PhaseCalculator void Editor::labelTextChanged(Label* labelThatHasChanged) { Node* processor = static_cast(getProcessor()); + + if (labelThatHasChanged == recalcIntervalEditable) + { + int intInput; + bool valid = updateControl(labelThatHasChanged, 0, INT_MAX, processor->calcInterval, intInput); + + if (valid) + { + processor->setParameter(RECALC_INTERVAL, static_cast(intInput)); + } + } + else if (labelThatHasChanged == arOrderEditable) + { + int intInput; + bool valid = updateControl(labelThatHasChanged, 1, INT_MAX, processor->arOrder, intInput); + + if (valid) + { + processor->setParameter(AR_ORDER, static_cast(intInput)); + } + } + else if (labelThatHasChanged == lowCutEditable) + { + float floatInput; + bool valid = updateControl(labelThatHasChanged, 0.0f, FLT_MAX, processor->lowCut, floatInput); + + if (valid) + { + processor->setParameter(LOWCUT, floatInput); + } + } + else if (labelThatHasChanged == highCutEditable) + { + float floatInput; + bool valid = updateControl(labelThatHasChanged, 0.0f, FLT_MAX, processor->highCut, floatInput); + + if (valid) + { + processor->setParameter(HIGHCUT, floatInput); + } + } } void Editor::channelChanged(int chan, bool newState) @@ -144,7 +265,10 @@ namespace PhaseCalculator void Editor::startAcquisition() { - bandBox->setEnabled(true); + bandBox->setEnabled(false); + lowCutEditable->setEnabled(false); + highCutEditable->setEnabled(false); + arOrderEditable->setEnabled(false); outputModeBox->setEnabled(false); channelSelector->inactivateButtons(); } @@ -152,6 +276,9 @@ namespace PhaseCalculator void Editor::stopAcquisition() { bandBox->setEnabled(true); + lowCutEditable->setEnabled(true); + highCutEditable->setEnabled(true); + arOrderEditable->setEnabled(true); outputModeBox->setEnabled(true); channelSelector->activateButtons(); } @@ -223,6 +350,8 @@ namespace PhaseCalculator Node* processor = (Node*)(getProcessor()); XmlElement* paramValues = xml->createNewChildElement("VALUES"); + paramValues->setAttribute("calcInterval", processor->calcInterval); + paramValues->setAttribute("arOrder", processor->arOrder); paramValues->setAttribute("lowCut", processor->lowCut); paramValues->setAttribute("highCut", processor->highCut); paramValues->setAttribute("outputMode", processor->outputMode); @@ -239,11 +368,26 @@ namespace PhaseCalculator forEachXmlChildElementWithTagName(*xml, xmlNode, "VALUES") { // some parameters have two fallbacks for backwards compatability + recalcIntervalEditable->setText(xmlNode->getStringAttribute("calcInterval", recalcIntervalEditable->getText()), sendNotificationSync); + arOrderEditable->setText(xmlNode->getStringAttribute("arOrder", arOrderEditable->getText()), sendNotificationSync); bandBox->setSelectedId(selectBandFromSavedParams(xmlNode) + 1, sendNotificationSync); + lowCutEditable->setText(xmlNode->getStringAttribute("lowCut", lowCutEditable->getText()), sendNotificationSync); + highCutEditable->setText(xmlNode->getStringAttribute("highCut", highCutEditable->getText()), sendNotificationSync); outputModeBox->setSelectedId(xmlNode->getIntAttribute("outputMode", outputModeBox->getSelectedId()), sendNotificationSync); } } + void Editor::refreshLowCut() + { + auto p = static_cast(getProcessor()); + lowCutEditable->setText(String(p->lowCut), dontSendNotification); + } + + void Editor::refreshHighCut() + { + auto p = static_cast(getProcessor()); + highCutEditable->setText(String(p->highCut), dontSendNotification); + } void Editor::refreshVisContinuousChan() { diff --git a/ASICPhaseCalculator/Source/PhaseCalculatorEditor.h b/PhaseCalculator/Source/PhaseCalculatorEditor.h similarity index 88% rename from ASICPhaseCalculator/Source/PhaseCalculatorEditor.h rename to PhaseCalculator/Source/PhaseCalculatorEditor.h index a9864e5..174f0e1 100644 --- a/ASICPhaseCalculator/Source/PhaseCalculatorEditor.h +++ b/PhaseCalculator/Source/PhaseCalculatorEditor.h @@ -62,6 +62,8 @@ namespace PhaseCalculator void loadCustomParameters(XmlElement* xml) override; // display updaters - do not trigger listeners. + void refreshLowCut(); + void refreshHighCut(); void refreshVisContinuousChan(); private: @@ -139,15 +141,30 @@ namespace PhaseCalculator ScopedPointer