diff --git a/PC_front.png b/PC_front.png deleted file mode 100644 index 50f3d26..0000000 Binary files a/PC_front.png and /dev/null differ diff --git a/PhaseCalculator_flow.pdf b/PhaseCalculator_flow.pdf deleted file mode 100644 index 9ed3d50..0000000 Binary files a/PhaseCalculator_flow.pdf and /dev/null differ diff --git a/PhaseCalculator_flow.xml b/PhaseCalculator_flow.xml deleted file mode 100644 index 07c26eb..0000000 --- a/PhaseCalculator_flow.xml +++ /dev/null @@ -1 +0,0 @@ -7V1bk6M2Fv41XZV9aBdIXB/7NrupndlMZbKV7NMWbdM2OxgcwNPd+fUrAcJIOtgCBO2+dKoytsz13PSdi44u8M326e9ZsNt8SVdhfIGM1dMFvr1AyLSQQ/6hI8/ViIvtamCdRav6oMPAt+ivsB406tF9tApz7sAiTeMi2vGDyzRJwmXBjQVZlj7yhz2kMX/XXbAOpYFvyyCWR3+PVsWmGvVs4zD+jzBab9idbaP+5T5Yfl9n6T6p73eB8EP5V/28Ddi16uPzTbBKH1tD+O4C32RpWlSftk83YUxpy8hWnfep49fmubMwKVROQG7NqR9BvA/ZMzsxOfl6Ff2gT1g811Rx/tzTx7ouwqfiMoijdXKBr8gRcfhQHH4ln9b1v+VV8l2QgJd5SJPiMi9ZT69iWrsn+Sqf0oy+zybIQ/Jvui92+4J+SOJndgPyctU9+PuS4fIF2CjiHgLRlyDjm2IbkwGTPmiRpd/DmzQmt6SH4JVz79gOfdIojlvjhJtouSTjNQ1uSwLg6x9hVkREhq7q4SLdkdHHTVSE33bBkt74kagLGSsFJKQsMJpno2eHT51sNBvhIEoXptuwyAgFDHaCW8tTrW8mYgL5eBBfm2nXpiW6zYFBrTLr5uIHsSEfasnpkCJsAlIkUBwiRIv8VBxqQ4Bsmbi7OIiSy3UWhmTwep0Fq4jQhnElSZNQYFQ9xHO1HpyA4pYrExw5pkxwSwO9sTVaaZfkPcOsW23vs54KNY69jKG3UUbseZTSn8IgL2TtW9mht7Lomwb5ptQhUHk9dI8dp/mFWXHUPLrEakAgurnvnOa+6QPqpoP7CLuntY03MSeYM5iUXbzRoGBIgcLWRBRmaOVd61e3HGhQIGyfZq83lf5IzP1a44ubIF7u46AoQccqKAKKHWOCz/qih+75iPGllg8ZNWyj1YreBtRZkcsaOIFszEMHmRHYtmVOOBo4YUJ6JpA6jkq6rVqCmxMcuOHJz5u7BlXTL+s4yPP68zLdRsv6syD31u2dZ7tHsMVxBdPACdu1OE5gDzB5EITzdbCi2w+gb66C4JEJIXh2lT0b+EKAHPXnNlkYrFrGbw8YxOrWAyH8K1JC0zU51nu+DCY9ZyIlhNDEJJy/+pW6wHWw4IP/Lf47hrWwBeVvRtpCYGBZCLAWIZAtb7hah9/qr2lWbNJ1mgTx3WFUsLotRoRPUfEHHV7Y9bf/1L/8LyyK55p6wb5IydDh2p9T6jKXx3USlRj/bMnCNF41VgTZOmSHWTXZ6OMfpX0Wkqk++sGHfEYR0XvHsDEIvYel5CRgCEg6Sy+8f5gGSBK/HNCbyaBkt+3UxfHaBj8ES/46N0QPInqi8a/wUT79ar+K0m/BdheH1/uHh/LADsNKp2MtMbPTBrcjUnbU3A6BdgepNnXBZN/j5czyZHDmATgZaRAzJqotMbtL/tyHe+q0BMmKihHVSholp55Lur8n5GwLYsXRn3bVUXkpFzn9tIkeipBegEau//YhCKehmuAvYQueqAGL42oQBfskWutrKTZRfE/e/isRjc9hsiZvcOb46zzkoLH9LIKBIDkwbQC16zAJDgbkQB+AEyKCQcbI5/TAcCLVGSpkCPE/DUJUQnssKtdGe8wB0gf26lO/plGpTU/MAgvxKk/gYfWg9VkHNl5lWfDcOmxHD8glRjcPqsb77vD/YBuQF2n2zJDCUP2fQM0ck593LQtBauZCcREdatYN6QeSOsh+yVZjiDy3kTWHGtkgW9ZfDT2yYHu8EjpYDpNMBcFcT6d5JTTInv9ofznYQl3WljPhSgbWAQys+3LetCGj3ttPvyloSBbm0V/BfcxErza65Gj7+sK+paJJ6JjXtDqNSqtEfk4UJErWv5UEv7SmMXcYu0IgyHVlcwflvXQIucnCy5qlvA4DNXI+AaaQI01KIo8MJMu8aXgvJvSmXEnxc0Lkp8xPvXnxRx4fByc2AMmTPeRb6RF/BMz2Ar3rfFSL3rswi8i96Bx7S73oezL/fj2MAfRP9wW9zE1TvWYMdGv4GNzN1Z13613IKSzy26fyD5qlJ8hj+fw07QD+MVQaY2rhIQSOBR4S4u7ox+U+i5+vM8IzqvqnqiSGQCL9DMrSIqgjsZempwtZ+bzeuazCiAtpyBzD2NHBsu5yC9UMlOlBGah/VKENat6zIMkf0mw7HHe/YrPqsOxuk14ErKo9kVVtgu9Tp5aOwY1+aSdh3j+dh2I51DaKqF2GWUCDQkGYjvLL+/LyXRUSfObnKATrk98REEEjpzOkd0zTlwh7U4XQWdS9Khii4IwVKQfJmgbfm2zLIWX3U3lSVdK8DdZJVOxX5EjyGEY5FpHBKAnK1yezyI6QNSn+9rZNk2WJ/GUp7nYBElT2oiV/58mK8tZsE8KQbfLnsk1IIVD9VmwTcoGA1FS2CcvhkZs0o3l6+gCEiDuivOTjOo6K5SbM37gdscUKGrMpjOFStFP5jsASiAM78m1QHk08huQxyOjU8b/9dpfLc0RQwd0go1kCI30g/7svo/LAsW+do4jnJ7ZkfoJeiR71kmMBd8WmXMt0HROP8TFNV9XcTb2PS8O+RNbbZgjmE99IkRt6tEueRcCi7YhqDS3b3sX7dVln+kDHpUqIYtOGaPSETRY+XLQ9zE1R0NWGV/QhCWjbhclluNs854ugoN5+FCSLhMyn6NNj9D0qCyZyMhOR1/r0y93v9P/kjAt0fUfPIf/+/d8/yw4qPaY+wKgOYM5pwA5ZDK7JeDuiNq/iy2qsBRLKWW/6rRUoLE9tAUUgN1R+FU8ZGkqXMZEyJ6dOt2Nx9Y+F+Ut0pNv7X6hCxdKFBmTlTd+eRnCm8yVOuQ4WkBDEL5cQtN7zsiy+vnamalrE4sJzODQM77XY+ymKS3aR+9O6NAMlq8u0XT7RYs31viDHPlJxp6A5SFY7GpQnCKC+CMxPuUiSOJNkIl+Gw6shX2rmBUtDuDoA+uVrQAmVlCMEu+qRHFNA680lTq2I1bGExZFdYc2KqrT8XAMZHYMPKYDLyieDPo5CHDlMVle0wwb5tixx8PJCOel/mNMcH4mzGgCGhuOPdtsDgFpsbCxM8Xlz6bKi8N4wRci/eSyLqR+mEAikE6ZMHLx0fBmBeBOVePqOWHoj9sbQCBbZiwFYBupmopwHFbqoDChP/JqFq4hGrToSpR1dVwSRGt5/AJ6iAGeFB0ihubJDFwJIvuPiQNO6dcvmq74dC1hs5MoGR8caPXxyie4Llf2+n8UfWChGBnKRTWsQ7TEwF8qkTLLYbGCzpp9ptLqWKSqNe/r8p6zIXE2adKw1aSq9vlVQ/BRM6WFXxIYyUJbbBwTL1iBY3vgeYKcEa1wXMFCw6hR8uaic0D7YUl9vE/yggd17itNb6flV57T1KuVQhyUzecRpA6lL0wBmMh2L2JDh68SgM602R4Yjg9Gm9mbykKjgIpgmFpwNjT4CUnCoPwph+65XwUApP4Lggo5C2GZ94usJR/f1DhHLuHKlLVXF+Azx6Ob25+bD1St6vwRJtJvTjdPno43qHIcx76N5rMjr1NSmw0lrENlEC3NVNeqYVg7TNtOHilwrHZh8tvPEMPHgyJorIGaNsybujuno7KVbjd20Sk+pIGThOgvznJK8H9bV8YCsqVOJw+lsn2tC1q9ieat9UKpRhss1+DJ923dkw8XKE3TnUJClNX2tbKcG4YphFswF0tl1IFa7AbMxj/0cMQTUYXfkCxlICEs7pmDCOmzhABPmKkxeryORZrmifysHVCDoraWqHMhknwUwvPr1S/B0RxQp3T3PCQxnje+PNMG+oG02xlDDCbZEWDt+tBXWL47JwaqbTg7n2R1e1TyJVpe1PGFMMYVLqMJBVwx7iTzTCAdtKKzal5Onp0egTBCKFvZ2rpkR49rt2C8qBcIyRgerTYSAOAlAy1NzCrq7+SiKrc/txkE+VFccLmKQxyGIGFvsTPhdREH8a7gsqiVqJ812NWVXMbmamUCwrd7fBDLM7q3huqBxZ6igvm57VbOJNMHpZt1507NJXtRsA/1URdsyyIyrlCS9lWXofhmnrVv9LCwa5tTCQLEbBAaWLWMo28qa2Y7joNa4zdE6cA0R1O7yqd6Gn4kuN//jDss/0qSbprhsyxaZp2rU5Uu5yJ9scmf99+aUDnmiPxPpmMqJxlhAa5Y/TDbEC7mmWu6s74SPfAGhOJonfEfBQZdn23Z4rZ5ztk9runXbgm5vsdwQaLAIkqS26f9Fsr3vmvRb84Cvad5GntBhFckrwlmUfQKzr1AD+1IUvtRFYsfgUyRQpBEBHq6WYu03E2MSu+bMGWM6Uj86uDHl17qK5j12yDH4yj4bqL+CVi9oaZBj9fDiHuLwqY4YXJ8KHnTt7NTCm9h0ZFTRscqwazso8hiNKawwCtvR88Cc3vGmKof7UpEGbPgL36YzSv2f0OtZ3HFHFYcQILIwDctu/rOOX7c3LBkFaRXEkO5tuuucgOrtWGsTcNFoR49tFYU6SMeT2+xbzEHnZiYdCww9hdjrOWys16bngW2njV2zNkuiKbg0SywwG2TcwL305p/utW4Nirr8nSOtEVxAaHUAAf9k87yB5UNjVwWcDxrooyJcqn2iCZ+FtV5YJXpn5cYbIDYNCslZaK0o1OMY++OJ32zS8i5XkQ/Z3HMAg8l0ssDe4Y9vCAuYR89YoPYfEJbXs1+TQlh+eowjBKN8Ob9sQR3vtOz257/jTcqaPT4l8TfKP9DiDVo/fxC04UYRWFvGDKD2ld/++NLF9yMVSJSKcfw/wu+pthD3J1hK+tZAYw8OzgEam/LrDxUFVVSLDTaFhg4gUpkKmqBzgCb8+/umnBeAwy9aKKDgF42nQE3oo1o82TuahgKXZ/D92tUhnSUkfKirSyUhCqvH37sVDeSBFk0DefBhSY9Z0v79oPyTaHYyS2pM3yjhNfNXm4vTWyiaOe4FXBxkdK/W/JCKQS6ODis/mYtjGh8+zmkfR52Fc/g4CGxxMzsGFhotGAgCwVPF5wwoPvfmACLuYGB3rc1JxwxPVWpjGt3BsQ9L0p+js1gSQyHt3LIkau/SZV7OUWrRkfnvA++AMtzIjLqXA3Smn86t+QgAzuHW9JcCMH89nR/zkcDT68f0Z/gl1PcRas+nJzDcbyqbBBRjoWDctFjfM26HUAgU4/EkMNkbT9/6sveS+v0qSr8F210cMjRndAC6dF8oqtN5Iz11dWFlrcD6jcn2cAabck3gP/XZbw6PpiTragH1IpgMcLA6tTfcQc6tJ1NuNVtVfSQzQX8HObd7Ln/RRiH7hGrFCff0LXeQE1vBm82Ks1mayHkTrS49J83zIc0b201OXfO6Yzovqnn5Nk2LzTvWPFfUPAy02J9Q8xTgw6tf/ml5/Go3x5c7Y0DLP009sEL7ioVdtWUFodLnMFnT/Z/m2srgDLoP6pIJxxDbJwDtij1o3bUOoWBpQO2bWzCJ6PLKfiIfiPtGSV/+Su9gOLaNHZWd1+eWF/WO2KAjN25dPicdtisHvZANGA094iEHAH4rt3BdkvcMk2az5JqVyX2+q0RB/K5xu5QuiSIEFW8cxHlKTtjndCcCo6Cfl63+rdVutMSKLaNyS9tdvbVtuf0zGQ0J9KEnhj/qVy03ta1P+xHleyJef4XZorfAtgNn57dhS8d2G1zUojVwUYX7dEi6sAjVkMN94OpjBh7GifoZ5MBdfibAQAbcZLU8XNhBbE06DCAoNM4bj8GAKJleDFZzsjt5hPlmMZcyCMPT0VihzKDe4UEm9aoVqCdOHR0/0negISO4w8JQIh7bZxKKLaLxRMMGBFIEorXowOzl5+A+jL+meVST7D4tinRLDojpD9fB8vu6tJYQzToCt1KSA7DrnvHpaLu+tkluciQGoD5RmruLiFiVfBGny+/0uEG87GV0LcuG+q5OhD8x23fzJa2uuF+IAyAsqFwAswKVMRSwTpbLnQ5+uNDmTflzstxkaZLuR7SWmWOLOTHsL7PvICXqdUv8TNoEKThzNVURjKnc0aNfpUQ3HdQrZGYRaraNqt446pEefUAzXtU4ag/xm3qvKYIxOKk1h+6+cfJCGndU9bqr/scZsDq6VHpgD/tiT907YxUUwesyZ5CXNNa48QDWBJr1gG6CjkhtdzX/OG7/M3wuf56LtawD8gtOVe2eWsCaKx3cUvDopk/Ka5F5oYE8UAAIdS7RETjvLv8bJ/HEOw66t3s+S1umx3rxXUmANBPESS0NgKZAJcbCbad3m27B7d+G5Hd70J7bbKcaaqd3te+UqdzuYB7rc2JJhR7zI+wRNaP9YWEV3WKLjlQlIE5s+2ws9QbE1oRiEjos/irdU/fs/Zl80xJsvqw9UJ8+HQjIRGdhhDSu69ITMVAFWlAj0KkMXbcnqbyPeal3FS0MGke/ZA9Jf0XGAtznnDBnF4f0RElDj1b6vCENFVxKyKNkm+po9yg17LwlpjBzMpGwAbi7nnpQZ57WxyYfg26o0jf0Iy59E6+jL/JjKXTsH8u3E2tI1JkIGzfkyAnWYzzus8EKfKOx7BSvo5GdJ1fbDkQ9wYtkIjorfposoIptjYUUonRVMcdY5Qn1mGShe6enmMIQ91gZYpNtCAL3Q01DClE6UqYn9sTqBkRAH/B28nYKvAQtTZ0KMB3Z73Rk6J2uAnyNwfYpsBDEUVDxtLgr3nh/5cw1z7pQXxWsQR/9+fQR+QpQ9qNQhUPIALe7S1UEUMUWys1QqKLQ738oPFouw5wWl9KS2iyN4zLheU9fmxLzNdrfkWwWazdcecmz6QOdyvUUJCnY3w8V7uTtcc46Bq/BlgWwdrJSsw/jPB1nbbZSoNkycj7Outo7zr3rPdksYduIOfdkw2D3pA8l1aOklisqqTNdrS/5mqVUqw7xKPqqX9JVSI/4Pw== \ No newline at end of file diff --git a/README.md b/README.md index 8669e32..c1c6e48 100644 --- a/README.md +++ b/README.md @@ -1,14 +1,15 @@ -# Phase Calculator Plugin [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.2548723.svg)](https://doi.org/10.5281/zenodo.2548723) +# Phase Calculator Plugin [![DOI](https://zenodo.org/badge/134900173.svg)](https://zenodo.org/badge/latestdoi/134900173) A plugin for the [Open Ephys GUI](https://github.com/open-ephys/plugin-GUI) to estimate the phase of inputs within a specified passband in real time. Its primary purpose is to enable closed-loop stimulation, typically in combination with the [Crossing Detector](https://github.com/tne-lab/crossing-detector) and either the Pulse Pal or an external stimulation system that receives ZeroMQ events (for example, the LabVIEW implementation [here](https://github.com/tne-lab/closed-loop-stim)). It can also output the magnitude or imaginary component of the band-limited analytic signal instead of the phase. (The "PH+MAG" mode outputs both phase and magnitude, in separate channels.) Finally, the visualization tab or window can receive TTL events and display the delayed but precise phase of a specified input at the event onset samples in a rose plot. This allows real-time monitoring of stimulation accuracy. -The following article describes the algorithm used by this plugin and the closed-loop stimulation pipeline as a whole: +The following article describes the algorithm used by the previous version of this plugin and the closed-loop stimulation pipeline as a whole: Blackwood, E., Lo, M., Widge, A. S. (2018). Continuous phase estimation for phase-locked neural stimulation using an autoregressive model for signal prediction. 40th International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), Honolulu, HI, 4736-4739. If you are just using the plugin in your project, you can cite this version of the code using the DOI listed in the header of this file. - + + ## Installation @@ -16,11 +17,11 @@ If you are just using the plugin in your project, you can cite this version of t ### Step 1: Decide which version to use. -* You are on the `master` branch, which uses the Fourier-transform-based [Hilbert transform](https://en.wikipedia.org/wiki/Hilbert_transform) over a sliding window to estimate the analytic signal (from which we derive the phase). We have been testing variants of this algorithm in our lab since mid-2016, and it is now fairly polished, but is somewhat less computationally efficient than the newer version. +* You are on the `master` branch, which uses a [Hilbert transformer](https://www.intechopen.com/books/matlab-a-fundamental-tool-for-scientific-computing-and-engineering-applications-volume-1/digital-fir-hilbert-transformers-fundamentals-and-efficient-design-methods) FIR filter - actually one of several, depending on the frequency band you are filtering to. This is more efficient than the original (published) algorithm since it doesn't require as much AR model-based prediction nor calculating an FFT on each step. -* If you want, you can switch to the `hilbert-transformer` branch, which uses a [Hilbert transformer](https://www.intechopen.com/books/matlab-a-fundamental-tool-for-scientific-computing-and-engineering-applications-volume-1/digital-fir-hilbert-transformers-fundamentals-and-efficient-design-methods) FIR filter instead - actually one of several, depending on the frequency band you are filtering to. This is more efficient since it doesn't require as much AR model-based prediction nor calculating an FFT on each step. However, it's not as well-tested. +* If you want, you can switch to the `old-version` branch, which uses the Fourier-transform-based [Hilbert transform](https://en.wikipedia.org/wiki/Hilbert_transform) over a sliding window to estimate the analytic signal (from which we derive the phase). We have been testing variants of this algorithm in our lab since mid-2016, and it is now fairly polished, but is somewhat less computationally efficient than the newer version. -* Alternatively, you can use the `tnel-development` (old algorithm) or `hilbert-transformer` (new algorithm) branch of [tne-lab/plugin-GUI](https://github.com/tne-lab/plugin-GUI/tree/tnel-development), which is a fork of the Open Ephys GUI source code with this plugin and some other useful ones already built in. Then you don't even need this repository or the rest of the instructions! (If you're on Linux though, you'll still have to install FFTW as described below.) Beware though, these are development branches and will be more unstable than `open-ephys/plugin-GUI/master`. +* Alternatively, you can use the `tnel-development` branch of [tne-lab/plugin-GUI](https://github.com/tne-lab/plugin-GUI/tree/tnel-development), which is a fork of the Open Ephys GUI source code with this plugin and some other useful ones already built in. Then you don't even need this repository or the rest of the instructions! (If you're on Linux though, you'll still have to install FFTW as described below.) Beware though, these are development branches and will be more unstable than `open-ephys/plugin-GUI/master`. ### All platforms: @@ -57,9 +58,9 @@ If you are just using the plugin in your project, you can cite this version of t * ***Important!*** Since the phase estimation algorithm is somewhat processor-intensive, by default all input channels are disabled. To enable the channels you would like to estimate the phase (or other component) of, select them in the "PARAM" section of the drawer. If "PH+MAG" is selcected as the output, this will also create the additional magnitude outputs for each selected input. -* Use "Low cut" and "High cut" to select the desired frequency passband. (Inputs should be unfiltered; the Phase Calculator uses its own bandpass filter.) +* In the "freq range" dropdown menu, choose a range of frequencies that includes the band you want to filter to. This determines which of the pre-designed Hilbert transformer filters is used internally (since if we tried to use one filter for all frequencies, it would end up with terrible performance everywhere). Note that the delta band is just too low to get a reasonably accurate phase estimate, even when downsampling to 500 Hz as this plugin does (before interpolating the output). -* The "Buffer length," "Past," and "Future" settings control the size and contents of the "Hilbert buffer" that is passed through the Hilbert transform to get the analytic (complex) signal for each channel. The "Past" length is the number of received samples to include on each iteration, including the current "processing buffer" (buffer of new data, whose size can be changed in the Audio Settings). The "Future" length is the number of samples of band-limited signal predicted using an autoregressive model to include. (See Blackwood, Lo & Widge, 2018 (once published) for more explanation of the algorithm.) The defaults are reasonable for a 4-8 Hz passband with a short processing buffer size. A longer buffer length generally increases the estimation accuracy at the cost of greater processing time. Also, one should optimally center the last (current) processing buffer of "past" data within the Hilbert buffer by setting "Future" to half of "Buffer Length" minus half of the processing buffer size. This ensures the maximal accuracy of the region of the Hilbert output corresponding to the current buffer. +* Use "Low cut" and "High cut" to select the desired frequency passband. (Inputs should be unfiltered; the Phase Calculator uses its own bandpass filter internally.) Changing the frequency range will automatically set a default high and low cut, but they can be changed to filter to any band within the range. * "AR Refresh" and "Order" control the autoregressive model used to predict the "future" portion of the Hilbert buffer. AR parameters are estimated using Burg's method. The default settings seem to work well, but other settings (particularly lower orders) may also work well. diff --git a/Source/ARModeler.h b/Source/ARModeler.h index a82e520..357c3c1 100644 --- a/Source/ARModeler.h +++ b/Source/ARModeler.h @@ -24,14 +24,9 @@ class ARModeler { public: - ARModeler() : arOrder(1), inputLength(2) + ARModeler(int order = 1, int length = 2, int strideIn = 1, bool* success = nullptr) { - reallocateStorage(); - } - - ARModeler(int order, int length, bool* success = nullptr) : ARModeler() - { - bool s = setParams(order, length); + bool s = setParams(order, length, strideIn); if (success != nullptr) { *success = s; @@ -41,15 +36,18 @@ class ARModeler { ~ARModeler() { } // returns true if successful. - bool setParams(int order, int length) + bool setParams(int order, int length, int strideIn) { - if (order < 1 || order >= length) + int newStridedLength = calcStridedLength(length, strideIn); + if (order < 1 || newStridedLength < order + 1) { jassertfalse; return false; } arOrder = order; inputLength = length; + stride = strideIn; + stridedLength = newStridedLength; reallocateStorage(); return true; } @@ -76,12 +74,12 @@ class ARModeler { double sn = 0.0; double sd = 0.0; int j; - int jj = inputLength - n; + int jj = stridedLength - n; for (j = 0; j < jj; j++) { - t1 = inputseries[j + n] + pef[j]; - t2 = inputseries[j] + per[j]; + t1 = inputseries[stride * (j + n)] + pef[j]; + t2 = inputseries[stride * j] + per[j]; sn -= 2.0 * t1 * t2; sd += (t1 * t1) + (t2 * t2); } @@ -99,8 +97,8 @@ class ARModeler { for (j = 0; j < jj; j++) { - per[j] = per[j] + t1 * pef[j] + t1 * inputseries[j + n]; - pef[j] = pef[j + 1] + t1 * per[j + 1] + t1 * inputseries[j + 1]; + per[j] = per[j] + t1 * pef[j] + t1 * inputseries[stride * (j + n)]; + pef[j] = pef[j + 1] + t1 * per[j + 1] + t1 * inputseries[stride * (j + 1)]; } } } @@ -116,13 +114,21 @@ class ARModeler { void resetPredictionError() { j_per.clearQuick(); - j_per.insertMultiple(0, 0, inputLength); + j_per.insertMultiple(0, 0, stridedLength); j_pef.clearQuick(); - j_pef.insertMultiple(0, 0, inputLength); + j_pef.insertMultiple(0, 0, stridedLength); + } + + static int calcStridedLength(int inputLength, int stride) + { + jassert(stride > 0); + return (inputLength + (stride - 1)) / stride; } int arOrder; int inputLength; + int stridedLength; + int stride; Array j_per; Array j_pef; Array j_h; diff --git a/Source/HTransformers.cpp b/Source/HTransformers.cpp new file mode 100644 index 0000000..457ed7d --- /dev/null +++ b/Source/HTransformers.cpp @@ -0,0 +1,155 @@ +/* +------------------------------------------------------------------ + +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 Hilbert +{ + namespace + { + const Array HIGH_GAM_VALID({ 60, 200 }); + const Array HIGH_GAM_DEFAULT({ 70, 150 }); + const Array HIGH_GAM_EXTREMA({ 81.6443, 123.1104, 169.3574 }); + const int HIGH_GAM_DELAY = 3; + // from Matlab: firls(6, [60 200]/250, [1 1], 'hilbert') + const double HIGH_GAM_TRANSFORMER[HIGH_GAM_DELAY] = { + -0.10383410506573287, + 0.0040553935691102303, + -0.59258484603659545 + }; + + const Array MID_GAM_VALID({ 40, 90 }); + const Array MID_GAM_DEFAULT({ 40, 90 }); + const Array MID_GAM_EXTREMA({ 64.4559 }); + const int MID_GAM_DELAY = 2; + // from Matlab: firls(4, [35 90]/250, [1 1], 'hilbert') + const double MID_GAM_TRANSFORMER[MID_GAM_DELAY] = { + -0.487176162115735, + -0.069437334858668653 + }; + + const Array LOW_GAM_VALID({ 30, 55 }); + const Array LOW_GAM_DEFAULT({ 30, 55 }); + const Array LOW_GAM_EXTREMA({ 43.3609 }); + const int LOW_GAM_DELAY = 2; + // from Matlab: firls(4, [30 55]/250, [1 1], 'hilbert') + const double LOW_GAM_TRANSFORMER[LOW_GAM_DELAY] = { + -1.5933788446351915, + 1.7241339075391682 + }; + + const Array BETA_VALID({ 10, 40 }); + const Array BETA_DEFAULT({ 12, 30 }); + const Array BETA_EXTREMA({ 21.5848 }); + const int BETA_DELAY = 9; + // from Matlab: firpm(18, [12 30 40 240]/250, [1 1 0.7 0.7], 'hilbert') + const double BETA_TRANSFORMER[BETA_DELAY] = { + -0.099949575596234311, + -0.020761484963254036, + -0.080803573080958854, + -0.027365064225587619, + -0.11114477443975329, + -0.025834076852645271, + -0.16664116044989324, + -0.015661948619847599, + -0.45268524264113719 + }; + + const Array ALPHA_THETA_VALID({ 4, 18 }); + const Array ALPHA_THETA_DEFAULT({ 4, 8 }); + const Array ALPHA_THETA_EXTREMA({ /* none */ }); + const int ALPHA_THETA_DELAY = 9; + // from Matlab: firpm(18, [4 246]/250, [1 1], 'hilbert') + const double ALPHA_THETA_TRANSFORMER[ALPHA_THETA_DELAY] = { + -0.28757250783614413, + 0.000027647225074994485, + -0.094611325643268351, + -0.00025887439499763831, + -0.129436276914844, + -0.0001608427426424053, + -0.21315096860055227, + -0.00055322197399797961, + -0.63685698210351149 + }; + + const String C_GAMMA{ L"\u03b3" }; + const String C_BETA{ L"\u03b2" }; + const String C_ALPHA_THETA{ L"\u03b1/\u03b8" }; + } + + String validBandToString(const Array& band) + { + jassert(band.size() == 2); + return " (" + String(band[0]) + "-" + String(band[1]) + " Hz)"; + } + + // exported constants + + extern const std::map BAND_NAME = { + { HIGH_GAM, "Hi " + C_GAMMA + validBandToString(HIGH_GAM_VALID) }, + { MID_GAM, "Mid " + C_GAMMA + validBandToString(MID_GAM_VALID) }, + { LOW_GAM, "Lo " + C_GAMMA + validBandToString(LOW_GAM_VALID) }, + { BETA, C_BETA + validBandToString(BETA_VALID) }, + { ALPHA_THETA, C_ALPHA_THETA + validBandToString(ALPHA_THETA_VALID) } + }; + + extern const std::map> VALID_BAND = { + { HIGH_GAM, HIGH_GAM_VALID }, + { MID_GAM, MID_GAM_VALID }, + { LOW_GAM, LOW_GAM_VALID }, + { BETA, BETA_VALID }, + { ALPHA_THETA, ALPHA_THETA_VALID } + }; + + extern const std::map> DEFAULT_BAND = { + { HIGH_GAM, HIGH_GAM_DEFAULT }, + { MID_GAM, MID_GAM_DEFAULT }, + { LOW_GAM, LOW_GAM_DEFAULT }, + { BETA, BETA_DEFAULT }, + { ALPHA_THETA, ALPHA_THETA_DEFAULT } + }; + + extern const std::map> EXTREMA = { + { HIGH_GAM, HIGH_GAM_EXTREMA }, + { MID_GAM, MID_GAM_EXTREMA }, + { LOW_GAM, LOW_GAM_EXTREMA }, + { BETA, BETA_EXTREMA }, + { ALPHA_THETA, ALPHA_THETA_EXTREMA } + }; + + extern const std::map DELAY = { + { HIGH_GAM, HIGH_GAM_DELAY }, + { MID_GAM, MID_GAM_DELAY }, + { LOW_GAM, LOW_GAM_DELAY }, + { BETA, BETA_DELAY }, + { ALPHA_THETA, ALPHA_THETA_DELAY } + }; + + extern const std::map TRANSFORMER = { + { HIGH_GAM, HIGH_GAM_TRANSFORMER }, + { MID_GAM, MID_GAM_TRANSFORMER }, + { LOW_GAM, LOW_GAM_TRANSFORMER }, + { BETA, BETA_TRANSFORMER }, + { ALPHA_THETA, ALPHA_THETA_TRANSFORMER } + }; +} \ No newline at end of file diff --git a/Source/HTransformers.h b/Source/HTransformers.h new file mode 100644 index 0000000..11e84f2 --- /dev/null +++ b/Source/HTransformers.h @@ -0,0 +1,78 @@ +/* +------------------------------------------------------------------ + +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 . + +*/ + +#ifndef H_TRANSFORMERS_H_INCLUDED +#define H_TRANSFORMERS_H_INCLUDED + +#include + +#include "../../../JuceLibraryCode/JuceHeader.h" + +/* + +Defines the Hilbert transformers appropriate to use for each frequency band. +(The actual coeficcients and other values are in the corresponding cpp file.) +- BAND_NAME: display name for each frequency band. +- VALID_BAND: range of frequencies appropriate to use with each transformer +- DEFAULT_BAND: band filled in by default when selecting each transformer +- EXTREMA: locations of local extrema of the magnitude response within VALID_BAND (if any). +- used to find the maximum and minimum response given a band of interest. +- DELAY: group delay of the transformer = # of samples to predict + also, the number of unique nonzero coefficients. +- TRANSFORMER: first DELAY coefficients for each filter. + the remaining DELAY+1 are 0, followed by the leading coefficients + again, negated and in reverse order. +*/ + +enum Band +{ + ALPHA_THETA = 0, + BETA, + LOW_GAM, + MID_GAM, + HIGH_GAM, + NUM_BANDS +}; + +namespace Hilbert +{ + const int FS = 500; + + extern const std::map BAND_NAME; + + // each is a pair (lower limit, upper limit) + extern const std::map> VALID_BAND; + + // each is a pair (low cut, high cut) + extern const std::map> DEFAULT_BAND; + + extern const std::map> EXTREMA; + + // samples of group delay (= order of filter / 2) + extern const std::map DELAY; + + // contain the first delay[band] coefficients; the rest are redundant and can be inferred + extern const std::map TRANSFORMER; +} + +#endif // H_TRANSFORMERS_H_INCLUDED \ No newline at end of file diff --git a/Source/PhaseCalculator.cpp b/Source/PhaseCalculator.cpp index 37b68e8..b5f781f 100644 --- a/Source/PhaseCalculator.cpp +++ b/Source/PhaseCalculator.cpp @@ -21,6 +21,9 @@ along with this program. If not, see . */ +#include // DBL_MAX +#include // sqrt + #include "PhaseCalculator.h" #include "PhaseCalculatorEditor.h" @@ -30,12 +33,6 @@ PhaseCalculator::PhaseCalculator() : GenericProcessor ("Phase Calculator") , Thread ("AR Modeler") , calcInterval (50) - , lowCut (4.0) - , highCut (8.0) - , minNyquist (FLT_MAX) - , hilbertLength (16384) - , predictionLength (8132) // = hilbertLength / 2 - (half of 120-sample buffer) - , haveSentWarning (false) , outputMode (PH) , visEventChannel (-1) , visContinuousChannel (-1) @@ -45,6 +42,7 @@ PhaseCalculator::PhaseCalculator() { setProcessorType(PROCESSOR_TYPE_FILTER); setAROrder(20); + setBand(ALPHA_THETA, true); } PhaseCalculator::~PhaseCalculator() {} @@ -98,18 +96,6 @@ void PhaseCalculator::setParameter(int parameterIndex, float newValue) int numInputs = getNumInputs(); switch (parameterIndex) { - case HILBERT_LENGTH: - setHilbertLength(static_cast(newValue)); - break; - - case PAST_LENGTH: - setPredLength(hilbertLength - static_cast(newValue)); - break; - - case PRED_LENGTH: - setPredLength(static_cast(newValue)); - break; - case RECALC_INTERVAL: calcInterval = static_cast(newValue); break; @@ -118,6 +104,10 @@ void PhaseCalculator::setParameter(int parameterIndex, float newValue) setAROrder(static_cast(newValue)); break; + case BAND: + setBand(static_cast(static_cast(newValue))); + break; + case LOWCUT: setLowCut(newValue); break; @@ -175,38 +165,18 @@ void PhaseCalculator::process(AudioSampleBuffer& buffer) { int chan = activeInputs[activeChan]; int nSamples = getNumSamples(chan); - if (nSamples == 0) + if (nSamples == 0) // nothing to do { continue; } - // Filter the data. - float* wpIn = buffer.getWritePointer(chan); - filters[activeChan]->process(nSamples, &wpIn); + // filter the data + float* const wpIn = buffer.getWritePointer(chan); + filters[chan]->process(nSamples, &wpIn); - // If there are more samples than we have room to process, process the most recent samples and output zero - // for the rest (this is an error that should be noticed and fixed). - int hilbertPastLength = hilbertLength - predictionLength; + // shift old data and copy new data into historyBuffer (as much as can fit) int historyStartIndex = jmax(nSamples - historyLength, 0); - int outputStartIndex = jmax(nSamples - hilbertPastLength, 0); - - jassert(outputStartIndex >= historyStartIndex); // since historyLength >= hilbertPastLength - int nSamplesToEnqueue = nSamples - historyStartIndex; - int nSamplesToProcess = nSamples - outputStartIndex; - - if (outputStartIndex != 0) - { - // clear the extra samples and send a warning message - buffer.clear(chan, 0, outputStartIndex); - if (!haveSentWarning) - { - CoreServices::sendStatusMessage("WARNING: Phase Calculator buffer is shorter than the sample buffer!"); - haveSentWarning = true; - } - } - - // shift old data and copy new data into historyBuffer int nOldSamples = historyLength - nSamplesToEnqueue; const double* rpBuffer = historyBuffer.getReadPointer(activeChan, nSamplesToEnqueue); @@ -220,14 +190,14 @@ void PhaseCalculator::process(AudioSampleBuffer& buffer) // shift old data for (int i = 0; i < nOldSamples; ++i) { - *wpBuffer++ = *rpBuffer++; + *(wpBuffer++) = *(rpBuffer++); } // copy new data - wpIn += historyStartIndex; + const float* rpIn = wpIn + historyStartIndex; for (int i = 0; i < nSamplesToEnqueue; ++i) { - *wpBuffer++ = *wpIn++; + *(wpBuffer++) = *(rpIn++); } } @@ -237,22 +207,15 @@ void PhaseCalculator::process(AudioSampleBuffer& buffer) bufferFreeSpace.set(activeChan, newBufferFreeSpace); if (newBufferFreeSpace == 0) { - // now that dataToProcess for this channel is full, + // now that the historyBuffer for this channel is full, // let the thread start calculating the AR model. chanState.set(activeChan, FULL_NO_AR); } } // calc phase and write out (only if AR model has been calculated) - if (chanState[activeChan] == FULL_AR) { - - // copy data to dataToProcess - if (hilbertPastLength > 0) - { - rpBuffer = historyBuffer.getReadPointer(activeChan, historyLength - hilbertPastLength); - hilbertBuffer[activeChan]->copyFrom(rpBuffer, hilbertPastLength); - } - + if (chanState[activeChan] == FULL_AR) + { // read current AR parameters safely Array localParams; localParams.resize(arOrder); @@ -266,22 +229,57 @@ void PhaseCalculator::process(AudioSampleBuffer& buffer) pLocalParam[i] = rpParam[i]; } } - - // use AR(20) model to predict upcoming data and append to dataToProcess - // get beyond end of history buffer indirectly to avoid juce assertion failure - rpBuffer = historyBuffer.getReadPointer(activeChan, historyLength - 1) + 1; - double* wpHilbert = hilbertBuffer[activeChan]->getRealPointer(hilbertPastLength); + // use AR model to fill predSamps (which is downsampled) based on past data. + int htDelay = Hilbert::DELAY.at(band); + int stride = sampleRateMultiple[chan]; - arPredict(rpBuffer, wpHilbert, predictionLength, pLocalParam, arOrder); + rpBuffer = historyBuffer.getReadPointer(activeChan, historyLength - dsOffset[chan]); + double* pPredSamps = predSamps.getRawDataPointer(); + arPredict(rpBuffer, pPredSamps, pLocalParam, htDelay + 1, stride, arOrder); - // Hilbert-transform hilbertBuffer - forwardPlan[activeChan]->execute(); - hilbertManip(hilbertBuffer[activeChan]); - backwardPlan[activeChan]->execute(); + // identify indices of current buffer to execute HT + htInds.clearQuick(); + for (int i = stride - dsOffset[chan]; 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, *htState[activeChan]); + if (kOut >= 0) + { + double rc = wpIn[htInds[kOut]]; + double ic = htScaleFactor * samp; + htOutput.set(kOut, std::complex(rc, ic)); + } + } - // calculate phase and write out to buffer - auto rpHilbert = hilbertBuffer[activeChan]->getComplexPointer(hilbertPastLength - nSamplesToProcess); + // copy state to transform prediction without changing the end-of-buffer state + htTempState = *htState[activeChan]; + + // 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) @@ -292,39 +290,94 @@ void PhaseCalculator::process(AudioSampleBuffer& buffer) wpOut2 = buffer.getWritePointer(outChan2); } - for (int i = 0; i < nSamplesToProcess; ++i) + kOut = 0; + std::complex prevCS = lastComputedSample[activeChan]; + std::complex nextCS = htOutput[kOut]; + double prevPhase, nextPhase, phaseSpan, thisPhase; + double prevMag, nextMag, magSpan, thisMag; + bool needPhase = outputMode != MAG; + bool needMag = outputMode != PH; + + if (needPhase) { + prevPhase = std::arg(prevCS); + nextPhase = std::arg(nextCS); + phaseSpan = circDist(nextPhase, prevPhase, Dsp::doublePi); + } + if (needMag) + { + prevMag = std::abs(prevCS); + nextMag = std::abs(nextCS); + magSpan = nextMag - prevMag; + } + int subSample = dsOffset[chan] % stride; + + for (int i = 0; i < nSamples; ++i, subSample = (subSample + 1) % stride) + { + if (subSample == 0) + { + // update interpolation frame + prevCS = nextCS; + nextCS = htOutput[++kOut]; + + if (needPhase) + { + prevPhase = nextPhase; + nextPhase = std::arg(nextCS); + phaseSpan = circDist(nextPhase, prevPhase, Dsp::doublePi); + } + if (needMag) + { + prevMag = nextMag; + nextMag = std::abs(nextCS); + magSpan = nextMag - prevMag; + } + } + + if (needPhase) + { + thisPhase = prevPhase + phaseSpan * subSample / stride; + thisPhase = circDist(thisPhase, 0, Dsp::doublePi); + } + if (needMag) + { + thisMag = prevMag + magSpan * subSample / stride; + } + switch (outputMode) { case MAG: - wpOut[i + outputStartIndex] = static_cast(std::abs(rpHilbert[i])); + wpOut[i] = static_cast(thisMag); break; - + case PH_AND_MAG: - wpOut2[i + outputStartIndex] = static_cast(std::abs(rpHilbert[i])); + wpOut2[i] = static_cast(thisMag); // fall through case PH: // output in degrees - wpOut[i + outputStartIndex] = static_cast(std::arg(rpHilbert[i]) * (180.0 / Dsp::doublePi)); + wpOut[i] = static_cast(thisPhase * (180.0 / Dsp::doublePi)); break; - + case IM: - wpOut[i + outputStartIndex] = static_cast(std::imag(rpHilbert[i])); + wpOut[i] = static_cast(thisMag * std::sin(thisPhase)); break; } } + lastComputedSample.set(activeChan, prevCS); + dsOffset.set(chan, ((dsOffset[chan] + nSamples - 1) % stride) + 1); // unwrapping / smoothing if (outputMode == PH || outputMode == PH_AND_MAG) { unwrapBuffer(wpOut, nSamples, activeChan); smoothBuffer(wpOut, nSamples, activeChan); + lastPhase.set(activeChan, wpOut[nSamples - 1]); } } else // fifo not full or AR model not ready { // just output zeros - buffer.clear(chan, outputStartIndex, nSamplesToProcess); + buffer.clear(chan, 0, nSamples); } // if this is the monitored channel for events, check whether we can add a new phase @@ -332,9 +385,6 @@ void PhaseCalculator::process(AudioSampleBuffer& buffer) { calcVisPhases(getTimestamp(chan) + getNumSamples(chan)); } - - // keep track of last sample - lastSample.set(activeChan, buffer.getSample(chan, nSamples - 1)); } } @@ -360,8 +410,6 @@ bool PhaseCalculator::disable() signalThreadShouldExit(); - haveSentWarning = false; - // reset states of active inputs Array activeInputs = getActiveInputs(); int nActiveInputs = activeInputs.size(); @@ -369,8 +417,11 @@ bool PhaseCalculator::disable() { bufferFreeSpace.set(activeChan, historyLength); chanState.set(activeChan, NOT_FULL); - lastSample.set(activeChan, 0); - filters[activeChan]->reset(); + FloatVectorOperations::clear(htState[activeChan]->begin(), htState[activeChan]->size()); + lastPhase.set(activeChan, 0); + lastComputedSample.set(activeChan, 0); + dsOffset.set(activeChan, sampleRateMultiple[activeChan]); + filters[activeInputs[activeChan]]->reset(); } // clear timestamp and phase queues @@ -397,7 +448,8 @@ void PhaseCalculator::run() Array paramsTemp; paramsTemp.resize(arOrder); - int numActiveChans = getActiveInputs().size(); + Array activeInputs = getActiveInputs(); + int numActiveChans = activeInputs.size(); uint32 startTime, endTime; while (!threadShouldExit()) @@ -423,7 +475,7 @@ void PhaseCalculator::run() // end critical section // calculate parameters - arModeler.fitModel(data, paramsTemp); + arModelers[activeInputs[activeChan]]->fitModel(data, paramsTemp); // write params safely { @@ -450,10 +502,44 @@ void PhaseCalculator::run() void PhaseCalculator::updateSettings() { - // handle changed sample rates - updateMinNyquist(); + // update arrays that store one entry per input + int numInputs = getNumInputs(); + int prevNumInputs = filters.size(); + int numInputsChange = numInputs - prevNumInputs; + + if (numInputsChange > 0) + { + // (temporary, until validateSampleRate call): + sampleRateMultiple.insertMultiple(-1, 1, numInputsChange); + dsOffset.insertMultiple(-1, 0, numInputsChange); + + // add new objects at new indices + for (int i = prevNumInputs; i < numInputs; i++) + { + filters.add(new BandpassFilter()); + // (temporary, until validateSampleRate call) + arModelers.add(new ARModeler()); + } + } + else if (numInputsChange < 0) + { + // delete unneeded entries + sampleRateMultiple.removeLast(-numInputsChange); + dsOffset.removeLast(-numInputsChange); + filters.removeLast(-numInputsChange); + arModelers.removeLast(-numInputsChange); + } + + // set filter parameters (sample rates may have changed) setFilterParameters(); + // check whether active channels can be processed + Array activeInputs = getActiveInputs(); + for (int chan : activeInputs) + { + validateSampleRate(chan); + } + // create new data channels if necessary updateSubProcessorMap(); updateExtraChannels(); @@ -550,6 +636,14 @@ void PhaseCalculator::loadCustomChannelParametersFromXml(XmlElement* channelElem } } +double PhaseCalculator::circDist(double x, double ref, double cutoff) +{ + const double TWO_PI = 2 * Dsp::doublePi; + double xMod = std::fmod(x - ref, TWO_PI); + double xPos = (xMod < 0 ? xMod + TWO_PI : xMod); + return (xPos > cutoff ? xPos - TWO_PI : xPos); +} + // ------------ PRIVATE METHODS --------------- void PhaseCalculator::handleEvent(const EventChannel* eventInfo, @@ -573,66 +667,91 @@ void PhaseCalculator::handleEvent(const EventChannel* eventInfo, } } -void PhaseCalculator::setHilbertLength(int newHilbertLength) +void PhaseCalculator::setAROrder(int newOrder) { - if (newHilbertLength == hilbertLength) { return; } - - float predLengthRatio = static_cast(predictionLength) / hilbertLength; - hilbertLength = newHilbertLength; - - static_cast(getEditor())->refreshHilbertLength(); + if (newOrder == arOrder) { return; } - // update dependent variables - int newPredLength = static_cast(roundf(predLengthRatio * hilbertLength)); - setPredLength(newPredLength); + arOrder = newOrder; + updateHistoryLength(); // update dependent per-channel objects - for (int i = 0; i < numActiveChansAllocated; i++) + int numInputs = getNumInputs(); + for (int chan = 0; chan < numInputs; ++chan) { - // processing buffers - hilbertBuffer[i]->resize(hilbertLength); + bool s = arModelers[chan]->setParams(arOrder, historyLength, sampleRateMultiple[chan]); + jassert(s); + } - // FFT plans - forwardPlan.set(i, new FFTWPlan(hilbertLength, hilbertBuffer[i], FFTW_MEASURE)); - backwardPlan.set(i, new FFTWPlan(hilbertLength, hilbertBuffer[i], FFTW_BACKWARD, FFTW_MEASURE)); + for (int i = 0; i < numActiveChansAllocated; i++) + { + arParams[i]->resize(arOrder); } } -void PhaseCalculator::setPredLength(int newPredLength) +void PhaseCalculator::setBand(Band newBand, bool force) { - if (newPredLength == predictionLength) { return; } + if (!force && newBand == band) { return; } + if (newBand < 0 || newBand >= NUM_BANDS) + { + jassertfalse; + return; + } + + band = newBand; - predictionLength = newPredLength; - static_cast(getEditor())->refreshPredLength(); + // set low and high cut to the defaults for this band, making sure to notify the editor + resetCutsToDefaults(); + + // resize htState for each active channel, htTempState, and predSamps + int delay = Hilbert::DELAY.at(band); + for (int i = 0; i < numActiveChansAllocated; ++i) + { + htState[i]->resize(delay * 2 + 1); + } + htTempState.resize(delay * 2 + 1); + predSamps.resize(delay + 1); } -void PhaseCalculator::setAROrder(int newOrder) +void PhaseCalculator::resetCutsToDefaults() { - if (newOrder == arOrder) { return; } - - arOrder = newOrder; - updateHistoryLength(); - bool s = arModeler.setParams(arOrder, historyLength); - jassert(s); + auto& defaultBand = Hilbert::DEFAULT_BAND.at(band); + lowCut = defaultBand[0]; + highCut = defaultBand[1]; - // update dependent per-channel objects - for (int i = 0; i < numActiveChansAllocated; i++) + updateScaleFactor(); + setFilterParameters(); + + auto editor = static_cast(getEditor()); + if (editor) { - arParams[i]->resize(arOrder); + editor->refreshLowCut(); + editor->refreshHighCut(); } } void PhaseCalculator::setLowCut(float newLowCut) { if (newLowCut == lowCut) { return; } + + auto editor = static_cast(getEditor()); + const Array& validBand = Hilbert::VALID_BAND.at(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 = lowCut + PASSBAND_EPS; - static_cast(getEditor())->refreshHighCut(); + highCut = jmin(lowCut + PASSBAND_EPS, validBand[1]); + editor->refreshHighCut(); } + updateScaleFactor(); setFilterParameters(); } @@ -640,13 +759,25 @@ void PhaseCalculator::setHighCut(float newHighCut) { if (newHighCut == highCut) { return; } + auto editor = static_cast(getEditor()); + const Array& validBand = Hilbert::VALID_BAND.at(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 = highCut - PASSBAND_EPS; - static_cast(getEditor())->refreshLowCut(); + lowCut = jmax(highCut - PASSBAND_EPS, validBand[0]); + editor->refreshLowCut(); } + updateScaleFactor(); setFilterParameters(); } @@ -654,9 +785,8 @@ void PhaseCalculator::setVisContChan(int newChan) { if (newChan >= 0) { - Array activeInputs = getActiveInputs(); - int visActiveChan = activeInputs.indexOf(newChan); - jassert(visActiveChan >= 0 && visActiveChan < filters.size()); + jassert(newChan < filters.size()); + jassert(getActiveInputs().indexOf(newChan) != -1); // disable event receival temporarily so we can flush the buffer int tempVisEventChan = visEventChannel; @@ -669,7 +799,7 @@ void PhaseCalculator::setVisContChan(int newChan) } // update filter settings - visReverseFilter.setParams(filters[visActiveChan]->getParams()); + visReverseFilter.setParams(filters[newChan]->getParams()); visEventChannel = tempVisEventChan; } visContinuousChannel = newChan; @@ -684,13 +814,19 @@ void PhaseCalculator::setVisContChan(int newChan) void PhaseCalculator::updateHistoryLength() { - jassert(VIS_HILBERT_LENGTH >= (1 << MAX_HILB_LEN_POW)); - int newHistoryLength = juce::jmax( - VIS_HILBERT_LENGTH, // must have enough samples for current and delayed Hilbert transforms - arOrder + 1); // must be longer than arOrder to train the AR model + Array activeInputs = getActiveInputs(); + + // minimum - must have enough samples to do a Hilbert transform on past values for visualization + int newHistoryLength = VIS_HILBERT_LENGTH; + for (int chan : activeInputs) + { + newHistoryLength = jmax(newHistoryLength, + arOrder * sampleRateMultiple[chan] + 1, // minimum to train AR model + Hilbert::FS * sampleRateMultiple[chan]); // use @ least 1 second to train model + } if (newHistoryLength == historyLength) { return; } - + historyLength = newHistoryLength; // update things that depend on historyLength @@ -701,54 +837,34 @@ void PhaseCalculator::updateHistoryLength() bufferFreeSpace.set(i, historyLength); } - bool s = arModeler.setParams(arOrder, historyLength); - jassert(s); -} - -void PhaseCalculator::updateMinNyquist() -{ - float currMinNyquist = FLT_MAX; - - auto ed = static_cast(getEditor()); - int nInputs = getNumInputs(); - Array activeInputs = getActiveInputs(); for (int chan : activeInputs) { - float sampleRate = getDataChannel(chan)->getSampleRate(); - currMinNyquist = jmin(currMinNyquist, sampleRate / 2); + bool success = arModelers[chan]->setParams(arOrder, historyLength, sampleRateMultiple[chan]); + jassert(success); } +} - minNyquist = currMinNyquist; - if (highCut > minNyquist) - { - // push down highCut to make it valid - CoreServices::sendStatusMessage("Lowering Phase Calculator upper passband limit to the Nyquist frequency (" + - String(minNyquist) + " Hz)"); - setHighCut(minNyquist); - ed->refreshHighCut(); - } +void PhaseCalculator::updateScaleFactor() +{ + htScaleFactor = getScaleFactor(band, lowCut, highCut); } void PhaseCalculator::setFilterParameters() { - Array activeInputs = getActiveInputs(); - int nActiveInputs = activeInputs.size(); - + int numInputs = getNumInputs(); + jassert(filters.size() == numInputs); double currLowCut = lowCut, currHighCut = highCut; jassert(currLowCut >= 0 && currLowCut < currHighCut); - for (int activeChan = 0; activeChan < nActiveInputs; ++activeChan) + for (int chan = 0; chan < numInputs; ++chan) { - jassert(activeChan < filters.size()); - int chan = activeInputs[activeChan]; - Dsp::Params params; params[0] = getDataChannel(chan)->getSampleRate(); // sample rate params[1] = 2; // order params[2] = (currHighCut + currLowCut) / 2; // center frequency params[3] = currHighCut - currLowCut; // bandwidth - filters[activeChan]->setParams(params); + filters[chan]->setParams(params); } } @@ -761,24 +877,51 @@ void PhaseCalculator::addActiveChannel() // simple arrays bufferFreeSpace.add(historyLength); chanState.add(NOT_FULL); - lastSample.add(0); + lastPhase.add(0); + lastComputedSample.add(0); // owned arrays - hilbertBuffer.add(new FFTWArray(hilbertLength)); - forwardPlan.add(new FFTWPlan(hilbertLength, hilbertBuffer.getLast(), FFTW_MEASURE)); - backwardPlan.add(new FFTWPlan(hilbertLength, hilbertBuffer.getLast(), FFTW_BACKWARD, FFTW_MEASURE)); historyLock.add(new CriticalSection()); arParamLock.add(new CriticalSection()); arParams.add(new Array()); arParams.getLast()->resize(arOrder); - filters.add(new BandpassFilter()); + htState.add(new Array()); + htState.getLast()->resize(Hilbert::DELAY.at(band) * 2 + 1); +} + +bool PhaseCalculator::validateSampleRate(int chan) +{ + auto e = getEditor(); + bool p, r, a; + e->getChannelSelectionState(chan, &p, &r, &a); + if (!p) { return false; } + + // test whether sample rate is a multiple of HT_FS + float fsMult = getDataChannel(chan)->getSampleRate() / Hilbert::FS; + float fsMultRound = std::round(fsMult); + if (std::abs(fsMult - fsMultRound) < FLT_EPSILON) + { + // leave selected + int fsMultInt = static_cast(fsMultRound); + sampleRateMultiple.set(chan, fsMultInt); + dsOffset.set(chan, fsMultInt); + bool s = arModelers[chan]->setParams(arOrder, historyLength, fsMultInt); + jassert(s); + return true; + } + + // deselect and send warning + deselectChannel(chan); + CoreServices::sendStatusMessage("Channel " + String(chan + 1) + " was deselected because" + + " its sample rate is not a multiple of " + String(Hilbert::FS)); + return false; } void PhaseCalculator::unwrapBuffer(float* wp, int nSamples, int activeChan) { for (int startInd = 0; startInd < nSamples - 1; startInd++) { - float diff = wp[startInd] - (startInd == 0 ? lastSample[activeChan] : wp[startInd - 1]); + float diff = wp[startInd] - (startInd == 0 ? lastPhase[activeChan] : wp[startInd - 1]); if (abs(diff) > 180) { // search forward for a jump in the opposite direction @@ -824,20 +967,20 @@ void PhaseCalculator::unwrapBuffer(float* wp, int nSamples, int activeChan) void PhaseCalculator::smoothBuffer(float* wp, int nSamples, int activeChan) { int actualGL = jmin(GLITCH_LIMIT, nSamples - 1); - float diff = wp[0] - lastSample[activeChan]; + float diff = wp[0] - lastPhase[activeChan]; if (diff < 0 && diff > -180) { // identify whether signal exceeds last sample of the previous buffer within glitchLimit samples. int endIndex = -1; for (int i = 1; i <= actualGL; i++) { - if (wp[i] > lastSample[activeChan]) + if (wp[i] > lastPhase[activeChan]) { endIndex = i; break; } // corner case where signal wraps before it exceeds lastSample - else if (wp[i] - wp[i - 1] < -180 && (wp[i] + 360) > lastSample[activeChan]) + else if (wp[i] - wp[i - 1] < -180 && (wp[i] + 360) > lastPhase[activeChan]) { wp[i] += 360; endIndex = i; @@ -848,10 +991,10 @@ void PhaseCalculator::smoothBuffer(float* wp, int nSamples, int activeChan) if (endIndex != -1) { // interpolate points from buffer start to endIndex - float slope = (wp[endIndex] - lastSample[activeChan]) / (endIndex + 1); + float slope = (wp[endIndex] - lastPhase[activeChan]) / (endIndex + 1); for (int i = 0; i < endIndex; i++) { - wp[i] = lastSample[activeChan] + (i + 1) * slope; + wp[i] = lastPhase[activeChan] + (i + 1) * slope; } } } @@ -878,7 +1021,7 @@ void PhaseCalculator::updateSubProcessorMap() uint16 subProcessorIdx = chanInfo->getSubProcessorIdx(); int procFullId = static_cast(getProcessorFullId(sourceNodeId, subProcessorIdx)); foundFullIds.add(procFullId); - + if (subProcessorMap.contains(procFullId)) { maxUsedIdx = jmax(maxUsedIdx, subProcessorMap[subProcessorIdx]); @@ -1040,16 +1183,18 @@ void PhaseCalculator::calcVisPhases(juce::int64 sdbEndTs) } } -void PhaseCalculator::arPredict(const double* readEnd, double* writeStart, int writeNum, const double* params, int order) +void PhaseCalculator::arPredict(const double* lastSample, double* prediction, + const double* params, int samps, int stride, int order) { - for (int s = 0; s < writeNum; ++s) + for (int s = 0; s < samps; ++s) { // s = index to write output - writeStart[s] = 0; + prediction[s] = 0; for (int ind = s - 1; ind > s - 1 - order; --ind) { // ind = index of previous output to read - writeStart[s] -= params[s - 1 - ind] * (ind < 0 ? readEnd[ind] : writeStart[ind]); + prediction[s] -= params[s - 1 - ind] * + (ind < 0 ? lastSample[(ind + 1) * stride] : prediction[ind]); } } } @@ -1071,7 +1216,7 @@ void PhaseCalculator::hilbertManip(FFTWArray* fftData) // normalize and double positive frequencies FloatVectorOperations::multiply(reinterpret_cast(wp + 1), 2.0 / n, numPosNegFreqDoubles); - + if (hasNyquist) { // normalize but don't double Nyquist frequency @@ -1081,3 +1226,72 @@ void PhaseCalculator::hilbertManip(FFTWArray* fftData) // set negative frequencies to 0 FloatVectorOperations::clear(reinterpret_cast(wp + firstNegFreq), numPosNegFreqDoubles); } + +double PhaseCalculator::getScaleFactor(Band band, double lowCut, double highCut) +{ + double maxResponse = -DBL_MAX; + double minResponse = DBL_MAX; + + Array testFreqs({ lowCut, highCut }); + // also look at any magnitude response extrema that fall within the selected band + for (double freq : Hilbert::EXTREMA.at(band)) + { + if (freq > lowCut && freq < highCut) + { + testFreqs.add(freq); + } + } + + // at each frequency, calculate the filter response + int nCoefs = Hilbert::DELAY.at(band); + for (double freq : testFreqs) + { + double normFreq = freq * Dsp::doublePi / (Hilbert::FS / 2); + std::complex response = 0; + + auto* transf = Hilbert::TRANSFORMER.at(band); + for (int kCoef = 0; kCoef < nCoefs; ++kCoef) + { + double coef = transf[kCoef]; + + // near component + response += coef * std::polar(1.0, -(kCoef * normFreq)); + + // mirrored component + // there is no term for -nCoefs because that coefficient is 0. + response -= coef * std::polar(1.0, -((2 * nCoefs - kCoef) * normFreq)); + } + + double absResponse = std::abs(response); + maxResponse = jmax(maxResponse, absResponse); + minResponse = jmin(minResponse, absResponse); + } + + // scale factor is reciprocal of geometric mean of max and min + return 1 / std::sqrt(minResponse * maxResponse); +} + +double PhaseCalculator::htFilterSamp(double input, Band band, Array& state) +{ + double* state_p = state.getRawDataPointer(); + + // initialize new state entry + int nCoefs = Hilbert::DELAY.at(band); + int order = nCoefs * 2; + jassert(order == state.size() - 1); + state_p[order] = 0; + + // incorporate new input + auto& transf = Hilbert::TRANSFORMER.at(band); + for (int kCoef = 0; kCoef < nCoefs; ++kCoef) + { + double val = input * transf[kCoef]; + state_p[kCoef] += val; // near component + state_p[order - kCoef] -= val; // mirrored component + } + + // output and shift state + double sampOut = state_p[0]; + memmove(state_p, state_p + 1, order * sizeof(double)); + return sampOut; +} diff --git a/Source/PhaseCalculator.h b/Source/PhaseCalculator.h index 3768ab9..c1a747d 100644 --- a/Source/PhaseCalculator.h +++ b/Source/PhaseCalculator.h @@ -45,17 +45,18 @@ accuracy of phase-locked stimulation in real time. #include #include // Filtering #include -#include "FFTWWrapper.h" // Fourier transform -#include "ARModeler.h" // Autoregressive modeling +#include + +#include "FFTWWrapper.h" // Fourier transform +#include "ARModeler.h" // Autoregressive modeling +#include "HTransformers.h" // Hilbert transformers & frequency bands // parameter indices enum Param { - HILBERT_LENGTH, - PAST_LENGTH, - PRED_LENGTH, RECALC_INTERVAL, AR_ORDER, + BAND, LOWCUT, HIGHCUT, OUTPUT_MODE, @@ -82,6 +83,26 @@ enum OutputMode { PH = 1, MAG, PH_AND_MAG, IM }; class PhaseCalculator : public GenericProcessor, public Thread { friend class PhaseCalculatorEditor; + + // ------- constants ------------ + + // default passband width if pushing lowCut down or highCut up to fix invalid range, + // and also the minimum for lowCut. + // (meant to be larger than actual minimum floating-point eps) + static const float PASSBAND_EPS; + + // priority of the AR model calculating thread (0 = lowest, 10 = highest) + static const int AR_PRIORITY = 3; + + // "glitch limit" (how long of a segment is allowed to be unwrapped or smoothed, in samples) + static const int GLITCH_LIMIT = 200; + + // process length for real-time visualization ground truth hilbert transform + // based on evaluation of phase error compared to offline processing + static const int VIS_HILBERT_LENGTH = 65536; + static const int VIS_TS_MIN_DELAY = 40000; + static const int VIS_TS_MAX_DELAY = 48000; + public: PhaseCalculator(); @@ -124,6 +145,13 @@ class PhaseCalculator : public GenericProcessor, public Thread // for visualizer continuous channel void saveCustomChannelParametersToXml(XmlElement* channelElement, int channelNumber, InfoObjectCommon::InfoObjectType channelType) override; void loadCustomChannelParametersFromXml(XmlElement* channelElement, InfoObjectCommon::InfoObjectType channelType) override; + + /* + * Circular distance between angles x and ref (in radians). The "cutoff" is the maximum + * possible positive output; greater differences will be wrapped to negative angles. + * For interpolation, and also RosePlot visualization. + */ + static double circDist(double x, double ref, double cutoff = 2 * Dsp::doublePi); private: @@ -133,15 +161,15 @@ class PhaseCalculator : public GenericProcessor, public Thread void handleEvent(const EventChannel* eventInfo, const MidiMessage& event, int samplePosition = 0) override; - // Sets hilbertLength (which in turn influences predLength and historyLength) - void setHilbertLength(int newHilbertLength); - - // Sets predLength (which in turn influences historyLength) - void setPredLength(int newPredLength); - // Sets arOrder (which in turn influences historyLength and the arModeler) void setAROrder(int newOrder); + // Sets frequency band (which resets lowCut and highCut to defaults for this band) + void setBand(Band newBand, bool force = false); + + // Resets lowCut and highCut to defaults for the current band + void resetCutsToDefaults(); + // Sets lowCut (which in turn influences highCut) void setLowCut(float newLowCut); @@ -155,12 +183,17 @@ class PhaseCalculator : public GenericProcessor, public Thread // VIS_HILBERT_LENGTH, hilbertLength, predictionLength, and arOrder). void updateHistoryLength(); - // Update minimum Nyquist frequency of active channels based on sample rates - void updateMinNyquist(); + // Update the htScaleFactor depending on the lowCut and highCut of the filter. + void updateScaleFactor(); // Update the filters of active channels. From FilterNode code. void setFilterParameters(); + // If given input channel has an incompatible sample rate, deselect it (and send a warning). + // Otherwise, save the downsampling multiplier and initialize downsampling phase. + // Returns the new selection state of the channel. + bool validateSampleRate(int chan); + // Allocate memory for a new active input channel void addActiveChannel(); @@ -195,13 +228,17 @@ class PhaseCalculator : public GenericProcessor, public Thread /* * 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 writeNum future data values starting at location writeStart. - * *** assumes there are at least 'order' existing data points *before* readEnd - * to use to calculate first 'order' data points. readEnd can equal writeStart. + * + * Writes samps future data values to prediction. */ - static void arPredict(const double* readEnd, double* writeStart, int writeNum, - const double* params, int order); + static void arPredict(const double* lastSample, double* prediction, + const double* params, int samps, int stride, int order); /* * hilbertManip: Hilbert transforms data in the frequency domain (including normalization by length of data). @@ -209,15 +246,17 @@ class PhaseCalculator : public GenericProcessor, public Thread */ static void hilbertManip(FFTWArray* fftData); - // ---- customizable parameters ------ + // 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. + static double getScaleFactor(Band band, double lowCut, double highCut); - // number of samples to pass through the Hilbert transform in the main calculation - int hilbertLength; + // Execute the hilbert transformer on one sample and update the state. + static double htFilterSamp(double input, Band band, Array& state); - // number of future values to predict (0 <= predictionLength <= hilbertLength - AR_ORDER) - int predictionLength; + // ---- customizable parameters ------ - // size of historyBuffer ( = max(VIS_HILBERT_LENGTH, hilbertLength - predictionLength)) + // size of historyBuffer ( >= VIS_HILBERT_LENGTH and long enough to train AR model effectively) int historyLength; // time to wait between AR model recalculations in ms @@ -228,6 +267,9 @@ class PhaseCalculator : public GenericProcessor, public Thread OutputMode outputMode; + // frequency band (determines which Hilbert transformer to use) + Band band; + // filter passband float highCut; float lowCut; @@ -242,9 +284,8 @@ class PhaseCalculator : public GenericProcessor, public Thread int numActiveChansAllocated = 0; - // Storage area for filtered data to be read by the main thread to fill hilbertBuffer, - // by the side thread to calculate AR model parameters, and by the visualization event - // handler to calculate acccurate phases at past event times. + // Storage area for filtered data to be used by the side thread to calculate AR model parameters + // and by the visualization event handler to calculate acccurate phases at past event times. AudioBuffer historyBuffer; // Keep track of how much of the historyBuffer is empty (per channel) @@ -253,38 +294,49 @@ class PhaseCalculator : public GenericProcessor, public Thread // Keeps track of each channel's state (see enum definition above) Array chanState; - // Buffers for FFTW processing - OwnedArray hilbertBuffer; - - // Plans for the FFTW Fourier Transform library - OwnedArray forwardPlan; // FFT - OwnedArray backwardPlan; // IFFT - // mutexes for sharedDataBuffer arrays, which are used in the side thread to calculate AR parameters. // since the side thread only READS sharedDataBuffer, only needs to be locked in the main thread when it's WRITTEN to. OwnedArray historyLock; // for autoregressive parameter calculation - ARModeler arModeler; + OwnedArray arModelers; - // keeps track of each channel's last output sample, to be used for smoothing. - Array lastSample; + // for active input channels, equals the sample rate divided by HT_FS + Array sampleRateMultiple; // AR model parameters OwnedArray> arParams; OwnedArray arParamLock; - // so that the warning message only gets sent once per run - bool haveSentWarning; + // storage area for predicted samples (to compensate for group delay) + Array predSamps; + + // state of each hilbert transformer (persistent and temporary) + OwnedArray> htState; + Array htTempState; + + // store indices of current buffer to feed into transformer + Array htInds; + + // store output of hilbert transformer + Array> htOutput; + + // approximate multiplier for the imaginary component output of the HT (depends on filter band) + double htScaleFactor; + + // keep track of each active channel's last non-interpolated ("computed") transformer output + // and the number of samples by which this sample precedes the start of the next buffer + // (ranges from 1 to sampleRateMultiple). + Array> lastComputedSample; + Array dsOffset; + + // keep track of last phase output, for glitch correction + Array lastPhase; // maps full source subprocessor IDs of incoming streams to indices of // corresponding subprocessors created here. HashMap subProcessorMap; - // for filtering - min Nyquist frequency over active channels - // (upper limit for highCut) - float minNyquist; - // delayed analysis for visualization FFTWArray visHilbertBuffer; FFTWPlan visForwardPlan, visBackwardPlan; @@ -316,29 +368,6 @@ class PhaseCalculator : public GenericProcessor, public Thread OwnedArray filters; BandpassFilter visReverseFilter; - // -------static------------ - - // default passband width if pushing lowCut down or highCut up to fix invalid range, - // and also the minimum for lowCut. - // (meant to be larger than actual minimum floating-point eps) - static const float PASSBAND_EPS; - - // priority of the AR model calculating thread (0 = lowest, 10 = highest) - static const int AR_PRIORITY = 3; - - // "glitch limit" (how long of a segment is allowed to be unwrapped or smoothed, in samples) - static const int GLITCH_LIMIT = 200; - - // process length limits (powers of 2) - static const int MIN_HILB_LEN_POW = 9; - static const int MAX_HILB_LEN_POW = 16; - - // process length for real-time visualization ground truth hilbert transform - // based on evaluation of phase error compared to offline processing - static const int VIS_HILBERT_LENGTH = 65536; - static const int VIS_TS_MIN_DELAY = 40000; - static const int VIS_TS_MAX_DELAY = 48000; - JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR(PhaseCalculator); }; diff --git a/Source/PhaseCalculatorCanvas.cpp b/Source/PhaseCalculatorCanvas.cpp index 62bf18c..042e5fd 100644 --- a/Source/PhaseCalculatorCanvas.cpp +++ b/Source/PhaseCalculatorCanvas.cpp @@ -491,7 +491,7 @@ void RosePlot::setReference(double newReference) void RosePlot::addAngle(double newAngle) { - newAngle = circDist(newAngle, 0.0); + newAngle = PhaseCalculator::circDist(newAngle, 0.0); angleData->insert(newAngle); rSum += std::exp(std::complex(0, newAngle)); repaint(); @@ -517,12 +517,8 @@ double RosePlot::getCircMean(bool usingReference) } double reference = usingReference ? referenceAngle : 0.0; - double meanRad = circDist(std::arg(rSum), reference); - // change range to [-90, 270), for ease of use - if (meanRad >= 3 * PI / 2) - { - meanRad -= 2 * PI; - } + // use range of (-90, 270] for ease of use + double meanRad = PhaseCalculator::circDist(std::arg(rSum), reference, 3 * PI / 2); return meanRad * 180 / PI; } @@ -545,12 +541,12 @@ void RosePlot::labelTextChanged(Label* labelThatHasChanged) double doubleInput; double currReferenceDeg = referenceAngle * 180.0 / PI; bool valid = PhaseCalculatorEditor::updateControl(labelThatHasChanged, - -DBL_MAX, DBL_MAX, currReferenceDeg, &doubleInput); + -DBL_MAX, DBL_MAX, currReferenceDeg, doubleInput); if (valid) { // convert to radians - double newReference = circDist(doubleInput * PI / 180.0, 0.0); + double newReference = PhaseCalculator::circDist(doubleInput * PI / 180.0, 0.0); labelThatHasChanged->setText(String(newReference * 180.0 / PI), dontSendNotification); setReference(newReference); canvas->updateStatLabels(); @@ -570,8 +566,10 @@ RosePlot::circularBinComparator::circularBinComparator(int numBinsIn, double ref bool RosePlot::circularBinComparator::operator() (const double& lhs, const double& rhs) const { - int lhsBin = static_cast(std::floor(circDist(lhs, referenceAngle) * numBins / (2 * PI))); - int rhsBin = static_cast(std::floor(circDist(rhs, referenceAngle) * numBins / (2 * PI))); + double lhsDist = PhaseCalculator::circDist(lhs, referenceAngle); + double rhsDist = PhaseCalculator::circDist(rhs, referenceAngle); + int lhsBin = static_cast(std::floor(lhsDist * numBins / (2 * PI))); + int rhsBin = static_cast(std::floor(rhsDist * numBins / (2 * PI))); return lhsBin < rhsBin; } @@ -602,12 +600,6 @@ void RosePlot::reorganizeAngleData() angleData.swapWith(newAngleData); } -double RosePlot::circDist(double x, double ref) -{ - double xMod = std::fmod(x - ref, 2 * PI); - return (xMod < 0 ? xMod + 2 * PI : xMod); -} - void RosePlot::updateAngles() { float step = 2 * PI_F / numBins; @@ -615,7 +607,7 @@ void RosePlot::updateAngles() for (int i = 0; i < numBins; ++i) { binMidpoints.set(i, step * (i + 0.5)); - float firstAngle = static_cast(circDist(PI / 2, step * (i + 1))); + float firstAngle = static_cast(PhaseCalculator::circDist(PI / 2, step * (i + 1))); segmentAngles.set(i, { firstAngle, firstAngle + step }); } } \ No newline at end of file diff --git a/Source/PhaseCalculatorCanvas.h b/Source/PhaseCalculatorCanvas.h index cb99c98..b6a78c6 100644 --- a/Source/PhaseCalculatorCanvas.h +++ b/Source/PhaseCalculatorCanvas.h @@ -89,9 +89,6 @@ class RosePlot AngleDataMultiset(int numBins, double referenceAngle, AngleDataMultiset* dataSource); }; - // circular distance in radians (return value is in [0, 2*pi)) - static double circDist(double x, double ref); - // make binMidpoints and segmentAngles reflect current numBins void updateAngles(); diff --git a/Source/PhaseCalculatorEditor.cpp b/Source/PhaseCalculatorEditor.cpp index 13a0f7e..9bffe20 100644 --- a/Source/PhaseCalculatorEditor.cpp +++ b/Source/PhaseCalculatorEditor.cpp @@ -22,130 +22,91 @@ along with this program. If not, see . #include "PhaseCalculatorEditor.h" #include "PhaseCalculatorCanvas.h" +#include "HTransformers.h" #include // INT_MAX -#include // stoi, stof, stod +#include // FLT_MAX +#include // abs PhaseCalculatorEditor::PhaseCalculatorEditor(PhaseCalculator* parentNode, bool useDefaultParameterEditors) - : VisualizerEditor (parentNode, 325, useDefaultParameterEditors) + : VisualizerEditor (parentNode, 220, useDefaultParameterEditors) , extraChanManager (parentNode) , prevExtraChans (0) { tabText = "Event Phase Plot"; - int filterWidth = 80; + int filterWidth = 120; // make the canvas now, so that restoring its parameters always works. canvas = new PhaseCalculatorCanvas(parentNode); - lowCutLabel = new Label("lowCutL", "Low cut"); - lowCutLabel->setBounds(10, 30, 80, 20); - lowCutLabel->setFont(Font("Small Text", 12, Font::plain)); + bandLabel = new Label("bandL", "Freq range:"); + bandLabel->setBounds(5, 25, 100, 20); + bandLabel->setFont({ "Small Text", 12, Font::plain }); + bandLabel->setColour(Label::textColourId, Colours::darkgrey); + addAndMakeVisible(bandLabel); + + bandBox = new ComboBox("bandB"); + for (int b = 0; b < NUM_BANDS; ++b) + { + bandBox->addItem(Hilbert::BAND_NAME.at(b), b + 1); + } + bandBox->setSelectedId(parentNode->band + 1); + bandBox->setTooltip(FREQ_RANGE_TOOLTIP); + bandBox->setBounds(7, 42, 110, 20); + 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(15, 47, 60, 18); + 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); - highCutLabel = new Label("highCutL", "High cut"); - highCutLabel->setBounds(10, 70, 80, 20); - highCutLabel->setFont(Font("Small Text", 12, Font::plain)); + 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(15, 87, 60, 18); + 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); - hilbertLengthLabel = new Label("hilbertLength", "Buffer length:"); - hilbertLengthLabel->setBounds(filterWidth + 8, 25, 180, 20); - hilbertLengthLabel->setFont(Font("Small Text", 12, Font::plain)); - hilbertLengthLabel->setColour(Label::textColourId, Colours::darkgrey); - addAndMakeVisible(hilbertLengthLabel); - - hilbertLengthBox = new ComboBox("Buffer size"); - hilbertLengthBox->setEditableText(true); - for (int pow = PhaseCalculator::MIN_HILB_LEN_POW; pow <= PhaseCalculator::MAX_HILB_LEN_POW; ++pow) - { - hilbertLengthBox->addItem(String(1 << pow), pow); - } - hilbertLengthBox->setText(String(parentNode->hilbertLength), dontSendNotification); - hilbertLengthBox->setTooltip(HILB_LENGTH_TOOLTIP); - hilbertLengthBox->setBounds(filterWidth + 10, 45, 80, 20); - hilbertLengthBox->addListener(this); - addAndMakeVisible(hilbertLengthBox); - - hilbertLengthUnitLabel = new Label("hilbertLengthUnit", "Samp."); - hilbertLengthUnitLabel->setBounds(filterWidth + 90, 45, 40, 20); - hilbertLengthUnitLabel->setFont(Font("Small Text", 12, Font::plain)); - hilbertLengthUnitLabel->setColour(Label::textColourId, Colours::darkgrey); - addAndMakeVisible(hilbertLengthUnitLabel); - - pastLengthLabel = new Label("pastLengthL", "Past:"); - pastLengthLabel->setBounds(filterWidth + 8, 85, 60, 15); - pastLengthLabel->setFont(Font("Small Text", 12, Font::plain)); - pastLengthLabel->setColour(Label::backgroundColourId, Colour(230, 168, 0)); - pastLengthLabel->setColour(Label::textColourId, Colours::darkgrey); - addAndMakeVisible(pastLengthLabel); - - predLengthLabel = new Label("predLengthL", "Future:"); - predLengthLabel->setBounds(filterWidth + 70, 85, 60, 15); - predLengthLabel->setFont(Font("Small Text", 12, Font::plain)); - predLengthLabel->setColour(Label::backgroundColourId, Colour(102, 140, 255)); - predLengthLabel->setColour(Label::textColourId, Colours::darkgrey); - addAndMakeVisible(predLengthLabel); - - pastLengthEditable = new Label("pastLengthE"); - pastLengthEditable->setEditable(true); - pastLengthEditable->addListener(this); - pastLengthEditable->setText(String(parentNode->hilbertLength - parentNode->predictionLength), dontSendNotification); - pastLengthEditable->setBounds(filterWidth + 8, 102, 60, 18); - pastLengthEditable->setColour(Label::backgroundColourId, Colours::grey); - pastLengthEditable->setColour(Label::textColourId, Colours::white); - addAndMakeVisible(pastLengthEditable); - - predLengthEditable = new Label("predLengthE"); - predLengthEditable->setEditable(true); - predLengthEditable->addListener(this); - predLengthEditable->setText(String(parentNode->predictionLength), dontSendNotification); - predLengthEditable->setBounds(filterWidth + 70, 102, 60, 18); - predLengthEditable->setColour(Label::backgroundColourId, Colours::grey); - predLengthEditable->setColour(Label::textColourId, Colours::white); - addAndMakeVisible(predLengthEditable); - - predLengthSlider = new Slider("predLength"); - predLengthSlider->setLookAndFeel(&v3LookAndFeel); - predLengthSlider->setSliderStyle(Slider::LinearBar); - predLengthSlider->setTextBoxStyle(Slider::NoTextBox, false, 40, 20); - predLengthSlider->setScrollWheelEnabled(false); - predLengthSlider->setBounds(filterWidth + 8, 70, 122, 10); - predLengthSlider->setColour(Slider::thumbColourId, Colour(255, 187, 0)); - predLengthSlider->setColour(Slider::backgroundColourId, Colour(51, 102, 255)); - predLengthSlider->setTooltip(PRED_LENGTH_TOOLTIP); - predLengthSlider->addListener(this); - predLengthSlider->setRange(0, parentNode->hilbertLength, 1); - predLengthSlider->setValue(parentNode->hilbertLength - parentNode->predictionLength, dontSendNotification); - addAndMakeVisible(predLengthSlider); + 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 + 140, 25, 100, 20); - recalcIntervalLabel->setFont(Font("Small Text", 12, Font::plain)); + 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 + 145, 44, 55, 18); + 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); @@ -153,21 +114,21 @@ PhaseCalculatorEditor::PhaseCalculatorEditor(PhaseCalculator* parentNode, bool u addAndMakeVisible(recalcIntervalEditable); recalcIntervalUnit = new Label("recalcU", "ms"); - recalcIntervalUnit->setBounds(filterWidth + 200, 47, 25, 15); - recalcIntervalUnit->setFont(Font("Small Text", 12, Font::plain)); + 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 + 140, 65, 60, 20); - arOrderLabel->setFont(Font("Small Text", 12, Font::plain)); + 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 + 195, 66, 25, 18); + 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); @@ -175,8 +136,8 @@ PhaseCalculatorEditor::PhaseCalculatorEditor(PhaseCalculator* parentNode, bool u addAndMakeVisible(arOrderEditable); outputModeLabel = new Label("outputModeL", "Output:"); - outputModeLabel->setBounds(filterWidth + 140, 87, 70, 20); - outputModeLabel->setFont(Font("Small Text", 12, Font::plain)); + outputModeLabel->setBounds(filterWidth, 87, 70, 20); + outputModeLabel->setFont({ "Small Text", 12, Font::plain }); outputModeLabel->setColour(Label::textColourId, Colours::darkgrey); addAndMakeVisible(outputModeLabel); @@ -187,7 +148,7 @@ PhaseCalculatorEditor::PhaseCalculatorEditor(PhaseCalculator* parentNode, bool u outputModeBox->addItem("IMAG", IM); outputModeBox->setSelectedId(parentNode->outputMode); outputModeBox->setTooltip(OUTPUT_MODE_TOOLTIP); - outputModeBox->setBounds(filterWidth + 145, 105, 76, 19); + outputModeBox->setBounds(filterWidth + 5, 105, 76, 19); outputModeBox->addListener(this); addAndMakeVisible(outputModeBox); @@ -201,21 +162,9 @@ void PhaseCalculatorEditor::comboBoxChanged(ComboBox* comboBoxThatHasChanged) { PhaseCalculator* processor = static_cast(getProcessor()); - if (comboBoxThatHasChanged == hilbertLengthBox) + if (comboBoxThatHasChanged == bandBox) { - int newId = hilbertLengthBox->getSelectedId(); - int newHilbertLength; - if (newId) // one of the items in the list is selected - { - newHilbertLength = (1 << newId); - } - else if (!updateControl(comboBoxThatHasChanged, 1 << PhaseCalculator::MIN_HILB_LEN_POW, - 1 << PhaseCalculator::MAX_HILB_LEN_POW, processor->hilbertLength, &newHilbertLength)) - { - return; - } - - processor->setParameter(HILBERT_LENGTH, static_cast(newHilbertLength)); + processor->setParameter(BAND, static_cast(bandBox->getSelectedId() - 1)); } else if (comboBoxThatHasChanged == outputModeBox) { @@ -227,34 +176,10 @@ void PhaseCalculatorEditor::labelTextChanged(Label* labelThatHasChanged) { PhaseCalculator* processor = static_cast(getProcessor()); - int sliderMax = static_cast(predLengthSlider->getMaximum()); - - if (labelThatHasChanged == pastLengthEditable) - { - int intInput; - bool valid = updateControl(labelThatHasChanged, 0, processor->hilbertLength, - processor->hilbertLength - processor->predictionLength, &intInput); - - if (valid) - { - processor->setParameter(PAST_LENGTH, static_cast(intInput)); - } - } - else if (labelThatHasChanged == predLengthEditable) + if (labelThatHasChanged == recalcIntervalEditable) { int intInput; - bool valid = updateControl(labelThatHasChanged, 0, processor->hilbertLength, - processor->predictionLength, &intInput); - - if (valid) - { - processor->setParameter(PRED_LENGTH, static_cast(intInput)); - } - } - else if (labelThatHasChanged == recalcIntervalEditable) - { - int intInput; - bool valid = updateControl(labelThatHasChanged, 0, INT_MAX, processor->calcInterval, &intInput); + bool valid = updateControl(labelThatHasChanged, 0, INT_MAX, processor->calcInterval, intInput); if (valid) { @@ -264,7 +189,7 @@ void PhaseCalculatorEditor::labelTextChanged(Label* labelThatHasChanged) else if (labelThatHasChanged == arOrderEditable) { int intInput; - bool valid = updateControl(labelThatHasChanged, 1, INT_MAX, processor->arOrder, &intInput); + bool valid = updateControl(labelThatHasChanged, 1, INT_MAX, processor->arOrder, intInput); if (valid) { @@ -274,8 +199,7 @@ void PhaseCalculatorEditor::labelTextChanged(Label* labelThatHasChanged) else if (labelThatHasChanged == lowCutEditable) { float floatInput; - bool valid = updateControl(labelThatHasChanged, PhaseCalculator::PASSBAND_EPS, - processor->minNyquist - PhaseCalculator::PASSBAND_EPS, processor->lowCut, &floatInput); + bool valid = updateControl(labelThatHasChanged, 0.0f, FLT_MAX, processor->lowCut, floatInput); if (valid) { @@ -285,35 +209,31 @@ void PhaseCalculatorEditor::labelTextChanged(Label* labelThatHasChanged) else if (labelThatHasChanged == highCutEditable) { float floatInput; - bool valid = updateControl(labelThatHasChanged, 2 * PhaseCalculator::PASSBAND_EPS, - processor->minNyquist, processor->highCut, &floatInput); - - if (valid) + bool valid = updateControl(labelThatHasChanged, 0.0f, FLT_MAX, processor->highCut, floatInput); + + if (valid) { processor->setParameter(HIGHCUT, floatInput); } } } -void PhaseCalculatorEditor::sliderEvent(Slider* slider) -{ - if (slider == predLengthSlider) - { - int newVal = static_cast(slider->getValue()); - int maxVal = static_cast(slider->getMaximum()); - getProcessor()->setParameter(PRED_LENGTH, static_cast(maxVal - newVal)); - } -} - void PhaseCalculatorEditor::channelChanged(int chan, bool newState) { - auto pc = static_cast(getProcessor()); + auto pc = static_cast(getProcessor()); if (chan < pc->getNumInputs()) { Array activeInputs = pc->getActiveInputs(); - if (newState && activeInputs.size() > pc->numActiveChansAllocated) - { - pc->addActiveChannel(); + if (newState) + { + // check whether sample rate is compatible (and if not, disable channel) + if (!pc->validateSampleRate(chan)) { return; } + + // ensure space allocated for per-active-channel arrays + if (activeInputs.size() > pc->numActiveChansAllocated) + { + pc->addActiveChannel(); + } } if (pc->outputMode == PH_AND_MAG) @@ -332,21 +252,14 @@ void PhaseCalculatorEditor::channelChanged(int chan, bool newState) } else { - // Can just do a partial update - pc->updateMinNyquist(); // minNyquist may have changed depending on active chans - pc->setFilterParameters(); // need to update in case the passband has changed - updateVisualizer(); // update the available continuous channels for visualizer + updateVisualizer(); // update the available continuous channels for visualizer } } } void PhaseCalculatorEditor::startAcquisition() { - GenericEditor::startAcquisition(); - hilbertLengthBox->setEnabled(false); - predLengthSlider->setEnabled(false); - pastLengthEditable->setEnabled(false); - predLengthEditable->setEnabled(false); + bandBox->setEnabled(false); lowCutEditable->setEnabled(false); highCutEditable->setEnabled(false); arOrderEditable->setEnabled(false); @@ -356,11 +269,7 @@ void PhaseCalculatorEditor::startAcquisition() void PhaseCalculatorEditor::stopAcquisition() { - GenericEditor::stopAcquisition(); - hilbertLengthBox->setEnabled(true); - predLengthSlider->setEnabled(true); - pastLengthEditable->setEnabled(true); - predLengthEditable->setEnabled(true); + bandBox->setEnabled(true); lowCutEditable->setEnabled(true); highCutEditable->setEnabled(true); arOrderEditable->setEnabled(true); @@ -433,15 +342,17 @@ void PhaseCalculatorEditor::saveCustomParameters(XmlElement* xml) xml->setAttribute("Type", "PhaseCalculatorEditor"); PhaseCalculator* processor = (PhaseCalculator*)(getProcessor()); - + XmlElement* paramValues = xml->createNewChildElement("VALUES"); - paramValues->setAttribute("hilbertLength", processor->hilbertLength); - paramValues->setAttribute("predLength", processor->predictionLength); 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); + + const Array& validBand = Hilbert::VALID_BAND.at(processor->band); + paramValues->setAttribute("rangeMin", validBand[0]); + paramValues->setAttribute("rangeMax", validBand[1]); } void PhaseCalculatorEditor::loadCustomParameters(XmlElement* xml) @@ -451,12 +362,9 @@ void PhaseCalculatorEditor::loadCustomParameters(XmlElement* xml) forEachXmlChildElementWithTagName(*xml, xmlNode, "VALUES") { // some parameters have two fallbacks for backwards compatability - hilbertLengthBox->setText(xmlNode->getStringAttribute("hilbertLength", - xmlNode->getStringAttribute("processLength", hilbertLengthBox->getText())), sendNotificationSync); - predLengthEditable->setText(xmlNode->getStringAttribute("predLength", - xmlNode->getStringAttribute("numFuture", predLengthEditable->getText())), sendNotificationSync); 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); @@ -475,34 +383,6 @@ void PhaseCalculatorEditor::refreshHighCut() highCutEditable->setText(String(p->highCut), dontSendNotification); } -void PhaseCalculatorEditor::refreshPredLength() -{ - auto p = static_cast(getProcessor()); - int newPredLength = p->predictionLength; - - jassert(predLengthSlider->getMinimum() == 0); - int maximum = static_cast(predLengthSlider->getMaximum()); - jassert(newPredLength >= 0 && newPredLength <= maximum); - - predLengthSlider->setValue(maximum - newPredLength, dontSendNotification); - pastLengthEditable->setText(String(maximum - newPredLength), dontSendNotification); - predLengthEditable->setText(String(newPredLength), dontSendNotification); -} - -void PhaseCalculatorEditor::refreshHilbertLength() -{ - auto p = static_cast(getProcessor()); - int newHilbertLength = p->hilbertLength; - - hilbertLengthBox->setText(String(newHilbertLength), dontSendNotification); - predLengthSlider->setRange(0, newHilbertLength, 1); - - // if possible, maintain prediction length while making past + pred = hilbertLength - int sliderVal = static_cast(predLengthSlider->getValue()); - pastLengthEditable->setText(String(sliderVal), dontSendNotification); - predLengthEditable->setText(String(newHilbertLength - sliderVal), dontSendNotification); -} - void PhaseCalculatorEditor::refreshVisContinuousChan() { auto p = static_cast(getProcessor()); @@ -515,24 +395,61 @@ void PhaseCalculatorEditor::refreshVisContinuousChan() // static utilities -template<> -int PhaseCalculatorEditor::fromString(const char* in) +int PhaseCalculatorEditor::selectBandFromSavedParams(const XmlElement* xmlNode) { - return std::stoi(in); -} + if (xmlNode->hasAttribute("rangeMin") && xmlNode->hasAttribute("rangeMax")) + { + // I don't trust JUCE's string parsing functionality (doesn't indicate failure, etc...) + float rangeMin, rangeMax; + if (readNumber(xmlNode->getStringAttribute("rangeMin"), rangeMin) && + readNumber(xmlNode->getStringAttribute("rangeMax"), rangeMax)) + { + // try to find a matching band. + for (auto keyVal : Hilbert::VALID_BAND) + { + int b = keyVal.first; + const Array& validRange = keyVal.second; -template<> -float PhaseCalculatorEditor::fromString(const char* in) -{ - return std::stof(in); -} + if (validRange[0] == rangeMin && validRange[1] == rangeMax) + { + return b; + } + } + } + } -template<> -double PhaseCalculatorEditor::fromString(const char* in) -{ - return std::stod(in); -} + // no match, or rangeMin and rangeMax were invalid + // try to find a good fit for lowCut and highCut + if (xmlNode->hasAttribute("lowCut") && xmlNode->hasAttribute("highCut")) + { + float lowCut, highCut; + if (readNumber(xmlNode->getStringAttribute("lowCut"), lowCut) && lowCut >= 0 && + readNumber(xmlNode->getStringAttribute("highCut"), highCut) && highCut > lowCut) + { + float midCut = (lowCut + highCut) / 2; + int bestBand = 0; + float bestBandDist = FLT_MAX; + + for (auto keyVal : Hilbert::VALID_BAND) + { + int b = keyVal.first; + const Array& validRange = keyVal.second; + + float midBand = (validRange[0] + validRange[1]) / 2; + float midDist = std::abs(midBand - midCut); + if (validRange[0] <= lowCut && validRange[1] >= highCut && midDist < bestBandDist) + { + bestBand = b; + bestBandDist = midDist; + } + } + return bestBand; + } + } + // just fallback to first band + return 0; +} // -------- ExtraChanManager --------- diff --git a/Source/PhaseCalculatorEditor.h b/Source/PhaseCalculatorEditor.h index 5cc749b..f9d658d 100644 --- a/Source/PhaseCalculatorEditor.h +++ b/Source/PhaseCalculatorEditor.h @@ -24,6 +24,7 @@ along with this program. If not, see . #ifndef PHASE_CALCULATOR_EDITOR_H_INCLUDED #define PHASE_CALCULATOR_EDITOR_H_INCLUDED +#include // string parsing #include #include "PhaseCalculator.h" @@ -44,9 +45,6 @@ class PhaseCalculatorEditor // implements Label::Listener void labelTextChanged(Label* labelThatHasChanged) override; - // overrides GenericEditor - void sliderEvent(Slider* slider) override; - // update display based on current channel void channelChanged(int chan, bool newState) override; @@ -64,16 +62,30 @@ class PhaseCalculatorEditor // display updaters - do not trigger listeners. void refreshLowCut(); void refreshHighCut(); - void refreshPredLength(); - void refreshHilbertLength(); void refreshVisContinuousChan(); private: - // Utilities for parsing entered values + /* + * Select a frequency band when loading based on the saved information. + * Don't want to use the index in case bands are added/removed in the future (and it's less readable). + * Instead, try to match the min and max valid frequencies of the range. If that fails (or if + * the rangeMin and rangeMax attributes don't exist), find the band that encompasses lowCut and highCut + * and puts them closest to the center. Finally, if all else fails, just return the first band. + */ + static int selectBandFromSavedParams(const XmlElement* xmlNode); + /* + * Tries to read a number of type out from input. Returns false if unsuccessful. + * Otherwise, returns true and writes the result to *out. + */ template - static T fromString(const char* in); + static bool readNumber(const String& input, T& out) + { + std::istringstream istream(input.toStdString()); + istream >> out; + return !istream.fail(); + } /* * Return whether the control contained a valid input between min and max, inclusive. @@ -83,26 +95,19 @@ class PhaseCalculatorEditor * In header to make sure specializations not used in PhaseCalculatorEditor.cpp * are still available to other translation units. */ - template static bool updateControl(Ctrl* c, const T min, const T max, - const T defaultValue, T* out) + const T defaultValue, T& out) { - T parsedVal; - try + if (readNumber(c->getText(), out)) { - parsedVal = fromString(c->getText().toRawUTF8()); + out = jmax(min, jmin(max, out)); + c->setText(String(out), dontSendNotification); + return true; } - catch (const std::logic_error&) - { - c->setText(String(defaultValue), dontSendNotification); - return false; - } - - *out = jmax(min, jmin(max, parsedVal)); - c->setText(String(*out), dontSendNotification); - return true; + c->setText(String(defaultValue), dontSendNotification); + return false; } // keep track of the record status of each "extra" channel @@ -131,21 +136,16 @@ class PhaseCalculatorEditor int prevExtraChans; ExtraChanManager extraChanManager; + ScopedPointer