diff --git a/.github/workflows/run_unix.yml b/.github/workflows/run_unix.yml index 557001dae..38d69c721 100644 --- a/.github/workflows/run_unix.yml +++ b/.github/workflows/run_unix.yml @@ -354,6 +354,8 @@ jobs: run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/transforms.py - name: Downsampling Python run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/downsampling.py + - name: PSD Welch Nyquist Overlap Python + run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/psd_welch_nyquist_overlap.py - name: ICA Python run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/ica.py - name: CSP Python diff --git a/.github/workflows/run_windows.yml b/.github/workflows/run_windows.yml index 244c3f0e5..8ca866a15 100644 --- a/.github/workflows/run_windows.yml +++ b/.github/workflows/run_windows.yml @@ -235,6 +235,9 @@ jobs: - name: Downsampling Python Test run: python %GITHUB_WORKSPACE%\python_package\examples\tests\downsampling.py shell: cmd + - name: PSD Welch Nyquist Overlap Python Test + run: python %GITHUB_WORKSPACE%\python_package\examples\tests\psd_welch_nyquist_overlap.py + shell: cmd - name: CSP Python Test run: python %GITHUB_WORKSPACE%\python_package\examples\tests\csp.py shell: cmd diff --git a/csharp_package/brainflow/brainflow/data_filter.cs b/csharp_package/brainflow/brainflow/data_filter.cs index b8f14f429..ea14748ec 100644 --- a/csharp_package/brainflow/brainflow/data_filter.cs +++ b/csharp_package/brainflow/brainflow/data_filter.cs @@ -751,7 +751,7 @@ public static Tuple get_psd (double[] data, int start_pos, i /// /// data for log PSD /// FFT Size - /// FFT Window overlap, must be between 0 and nfft + /// FFT Window overlap, must be >= 0 and < nfft /// sampling rate /// window function /// Tuple of ampls and freqs arrays @@ -1265,7 +1265,7 @@ public static unsafe Tuple get_psd (double[,] data, int row_ /// data for log PSD /// /// FFT Size - /// FFT Window overlap, must be between 0 and nfft + /// FFT Window overlap, must be >= 0 and < nfft /// sampling rate /// window function /// Tuple of ampls and freqs arrays diff --git a/java_package/brainflow/src/main/java/brainflow/DataFilter.java b/java_package/brainflow/src/main/java/brainflow/DataFilter.java index cff1fc7ba..13bb57db3 100644 --- a/java_package/brainflow/src/main/java/brainflow/DataFilter.java +++ b/java_package/brainflow/src/main/java/brainflow/DataFilter.java @@ -960,7 +960,7 @@ public static Pair get_psd (double[] data, int start_pos, in * * @param data data to process * @param nfft size of FFT, must be even - * @param overlap overlap between FFT Windows, must be between 0 and nfft + * @param overlap overlap between FFT Windows, must be at least 0 and less than nfft * @param sampling_rate sampling rate * @param window window function * @return pair of ampl and freq arrays @@ -988,7 +988,7 @@ public static Pair get_psd_welch (double[] data, int nfft, i * * @param data data to process * @param nfft size of FFT, must be even - * @param overlap overlap between FFT Windows, must be between 0 and nfft + * @param overlap overlap between FFT Windows, must be at least 0 and less than nfft * @param sampling_rate sampling rate * @param window window function * @return pair of ampl and freq arrays diff --git a/python_package/brainflow/data_filter.py b/python_package/brainflow/data_filter.py index d98a3dda2..4e3ff35bb 100644 --- a/python_package/brainflow/data_filter.py +++ b/python_package/brainflow/data_filter.py @@ -1105,7 +1105,7 @@ def get_psd_welch(cls, data, nfft: int, overlap: int, sampling_rate: int, window :type data: NDArray[Shape["*"], Float64] :param nfft: FFT Window size, must be even :type nfft: int - :param overlap: overlap of FFT Windows, must be between 0 and nfft + :param overlap: overlap of FFT Windows, must be >= 0 and < nfft :type overlap: int :param sampling_rate: sampling rate :type sampling_rate: int diff --git a/python_package/examples/tests/psd_welch_nyquist_overlap.py b/python_package/examples/tests/psd_welch_nyquist_overlap.py new file mode 100644 index 000000000..96844c203 --- /dev/null +++ b/python_package/examples/tests/psd_welch_nyquist_overlap.py @@ -0,0 +1,71 @@ +import threading + +import numpy as np + +from brainflow.data_filter import DataFilter, WindowOperations +from brainflow.exit_codes import BrainFlowError, BrainFlowExitCodes + +HANG_TIMEOUT_SECONDS = 60 + + +def call_with_timeout(fn): + """Run fn on a worker thread so a non-terminating call fails instead of hanging. + + ctypes releases the GIL around the native call, so the join below still returns while + the worker spins. Without this, a regression of the overlap guard would stall the whole + CI job rather than reporting a failure. + """ + outcome = {} + + def run(): + try: + fn() + outcome['returned'] = True + except BrainFlowError as err: + outcome['exit_code'] = err.exit_code + + worker = threading.Thread(target=run, daemon=True) + worker.start() + worker.join(HANG_TIMEOUT_SECONDS) + return worker.is_alive(), outcome + + +def main(): + nfft = 32 + sampling_rate = 128 + segments = 4 + window = WindowOperations.HANNING.value + + rng = np.random.default_rng(7) + data = rng.standard_normal(nfft * segments) + + # overlap == nfft used to leave the segment cursor stationary, so the averaging loop + # never advanced and the call spun forever. The valid range is 0 <= overlap < nfft. + still_running, outcome = call_with_timeout( + lambda: DataFilter.get_psd_welch(np.copy(data), nfft, nfft, sampling_rate, window) + ) + assert not still_running, 'get_psd_welch did not terminate with overlap equal to nfft' + assert outcome.get('exit_code') == BrainFlowExitCodes.INVALID_ARGUMENTS_ERROR.value, outcome + + # With no overlap the segments are disjoint, so each Welch bin is exactly the mean of + # the per-segment PSD bins and can be checked against get_psd directly. + ampl, _ = DataFilter.get_psd_welch(np.copy(data), nfft, 0, sampling_rate, window) + per_segment = np.array( + [ + DataFilter.get_psd(np.copy(data[i * nfft:(i + 1) * nfft]), sampling_rate, window)[0] + for i in range(segments) + ] + ) + expected = np.mean(per_segment, axis=0) + + assert ampl.shape[0] == nfft // 2 + 1, ampl.shape + # The averaging loop stopped at nfft / 2, so the final Nyquist bin kept the running sum + # over all segments instead of the average and read `segments` times too large. + assert np.allclose(ampl[-1], expected[-1]), (ampl[-1], expected[-1]) + assert np.allclose(ampl, expected), (ampl, expected) + + print('psd welch nyquist and overlap regression passed') + + +if __name__ == '__main__': + main() diff --git a/src/data_handler/data_handler.cpp b/src/data_handler/data_handler.cpp index 106c614dd..eae13e19a 100644 --- a/src/data_handler/data_handler.cpp +++ b/src/data_handler/data_handler.cpp @@ -1239,7 +1239,7 @@ int get_psd_welch (double *data, int data_len, int nfft, int overlap, int sampli int window_function, double *output_ampl, double *output_freq) { if ((data == NULL) || (data_len < 1) || (nfft & (nfft - 1)) || (output_ampl == NULL) || - (output_freq == NULL) || (sampling_rate < 1) || (overlap < 0) || (overlap > nfft)) + (output_freq == NULL) || (sampling_rate < 1) || (overlap < 0) || (overlap >= nfft)) { data_logger->error ("Please review your arguments."); return (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR; @@ -1270,7 +1270,7 @@ int get_psd_welch (double *data, int data_len, int nfft, int overlap, int sampli return (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR; } // average data - for (int i = 0; i < nfft / 2; i++) + for (int i = 0; i < nfft / 2 + 1; i++) { output_ampl[i] /= counter; }