Fix BLACKMAN_HARRIS window using Nuttall's coefficients

.github/workflows/run_unix.yml

@@ -374,6 +374,8 @@ jobs:
       run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/csp.py
     - name: Windowing Python
       run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/windowing.py
+    - name: Blackman-Harris Window Regression Python
+      run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/blackman_harris_window.py
     - name: MNE Python
       run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/brainflow_to_mne.py
     - name: BandPower Python

.github/workflows/run_windows.yml

@@ -267,6 +267,9 @@ jobs:
     - name: Windowing Python Test
       run: python %GITHUB_WORKSPACE%\python_package\examples\tests\windowing.py
       shell: cmd
+    - name: Blackman-Harris Window Regression Python Test
+      run: python %GITHUB_WORKSPACE%\python_package\examples\tests\blackman_harris_window.py
+      shell: cmd
     - name: BrainFlow to MNE Python Test
       run: python %GITHUB_WORKSPACE%\python_package\examples\tests\brainflow_to_mne.py
       shell: cmd

docs/UserAPI.rst

@@ -74,7 +74,7 @@ PSD and Welch return one-sided power density in input-units squared per Hz, norm
 by ``sampling_rate * sum(window**2)``. DC and Nyquist are not doubled. This corrects
 the historical underestimation with nonrectangular windows and integer overflow for
 large FFTs. Absolute powers therefore change. Windows are periodic; the four-term
-Blackman-Harris coefficients are 0.355768, 0.487396, 0.144232, and 0.012604.
+Blackman-Harris coefficients are 0.35875, 0.48829, 0.14128, and 0.01168.
 
 Public filters require positive finite cutoffs strictly below Nyquist, ascending
 band edges, order 1-8, and positive finite Chebyshev ripple. Zero-phase methods

python_package/examples/tests/blackman_harris_window.py

@@ -0,0 +1,34 @@
+import numpy as np
+
+from brainflow.data_filter import DataFilter, WindowOperations
+
+
+def four_term_window(coefficients, window_len):
+    a0, a1, a2, a3 = coefficients
+    phase = 2.0 * np.pi * np.arange(window_len) / window_len
+    return a0 - a1 * np.cos(phase) + a2 * np.cos(2 * phase) - a3 * np.cos(3 * phase)
+
+
+def main():
+    # Four-term Blackman-Harris from F. J. Harris (1978), periodic like the other windows.
+    # scipy.signal.windows.blackmanharris(n, sym=False) gives the same values.
+    blackman_harris = (0.35875, 0.48829, 0.14128, 0.01168)
+    # Nuttall's continuous first derivative window, which BLACKMAN_HARRIS used to return.
+    nuttall = (0.355768, 0.487396, 0.144232, 0.012604)
+
+    for window_len in (8, 9, 20, 64, 256):
+        window = np.array(DataFilter.get_window(WindowOperations.BLACKMAN_HARRIS.value, window_len))
+        expected = four_term_window(blackman_harris, window_len)
+        assert np.allclose(window, expected, rtol=0, atol=1e-12), (window_len, window, expected)
+        # The two windows differ by up to about 0.006, far above the tolerance above.
+        assert not np.allclose(window, four_term_window(nuttall, window_len), rtol=0, atol=1e-4), window_len
+
+    # Blackman-Harris does not reach zero at its first sample, unlike Nuttall's window.
+    window = DataFilter.get_window(WindowOperations.BLACKMAN_HARRIS.value, 8)
+    assert np.isclose(window[0], 6e-05, rtol=0, atol=1e-12), window[0]
+
+    print('blackman-harris window regression passed')
+
+
+if __name__ == '__main__':
+    main()

src/data_handler/inc/window_functions.h

@@ -31,12 +31,13 @@ inline void hamming_function (int window_len, double *wind)
     }
 }
 
+// four-term Blackman-Harris, F. J. Harris (1978)
 inline void blackman_harris_function (int window_len, double *wind)
 {
     for (int i = 0; i < window_len; i++)
     {
-        wind[i] = 0.355768 - 0.487396 * cos (2.0 * M_PI * i / window_len) +
-            0.144232 * cos (4.0 * M_PI * i / window_len) -
-            0.012604 * cos (6.0 * M_PI * i / window_len);
+        wind[i] = 0.35875 - 0.48829 * cos (2.0 * M_PI * i / window_len) +
+            0.14128 * cos (4.0 * M_PI * i / window_len) -
+            0.01168 * cos (6.0 * M_PI * i / window_len);
     }
 }