Merge branch 'hilbert-transformer'

README.md

@@ -1,26 +1,27 @@
-# 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.
 
-<img src="PC_front.png" width="328" /> <img src="PC_vis.png" width="350" />
+<img src="ht_pc.png" width="200" /><img src="ht_pc_menu1.png" width="200"/><img src="ht_pc_menu2.png" width="200"/>
+<img src="PC_vis.png" width="350" />
 
 ## Installation
 
 (Instructions here are for Windows and Linux - Mac support may be coming soon.)
 
 ### 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.
 

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<double> j_per;
     Array<double> j_pef;
     Array<double> j_h;