Skip to content

Fix get_psd_welch Nyquist bin averaging and infinite loop when overlap equals nfft - #851

Open
adityasingh2400 wants to merge 4 commits into
brainflow-dev:masterfrom
adityasingh2400:fix-psd-welch-nyquist-overlap
Open

Fix get_psd_welch Nyquist bin averaging and infinite loop when overlap equals nfft#851
adityasingh2400 wants to merge 4 commits into
brainflow-dev:masterfrom
adityasingh2400:fix-psd-welch-nyquist-overlap

Conversation

@adityasingh2400

Copy link
Copy Markdown

get_psd_welch has two independent defects. Both are in the same function, so they are here together.

The Nyquist bin is never averaged. The accumulation loop sums nfft / 2 + 1 bins across every segment, but the averaging loop divides only nfft / 2 of them. The Nyquist bin at index nfft / 2 is summed and never divided by the segment count, so it comes back inflated by exactly the number of segments. For a stationary signal that means the value grows without bound as you feed in more data, which is the opposite of what a Welch average is for. get_band_power integrates up to data_len - 1, so any band whose upper edge reaches the Nyquist frequency inherits the error.

overlap == nfft hangs the process. Argument validation accepts it, and the segment loop then advances by pos += (nfft - overlap), which is zero. pos never moves, the loop condition never goes false, and the call spins forever inside the library with no way for the caller to recover.

The fix averages nfft / 2 + 1 bins instead of nfft / 2, and rejects overlap >= nfft instead of overlap > nfft. Two characters each. The Python, C# and Java docstrings said the valid range was 0 to nfft inclusive, so they are corrected too.

Verified by building libDataHandler.dylib on macOS arm64 and calling the C API directly, comparing against a reference Welch average computed by calling get_psd on each segment and averaging every bin including Nyquist.

Nyquist bin, nfft=64, overlap=0, sampling_rate=256, signal containing a tone exactly at the 128 Hz Nyquist frequency:

data_len segments before after reference
64 1 0.25000118 0.25000118 0.25000118
128 2 0.50000000 0.25000000 0.25000000
256 4 1.00000000 0.25000000 0.25000000
384 6 1.50000000 0.25000000 0.25000000
512 8 2.00000000 0.25000000 0.25000000

The before to reference ratio is 1.0000, 2.0000, 3.0000 up to 8.0000, exactly the segment count. Every other bin already matched the reference to 0.000e+00 absolute difference before the fix, which isolates the defect to the single Nyquist bin.

Propagation into get_band_power, 6 segments:

band before after reference
4 to 8 Hz (theta) 0.13650396 0.13650396 0.13650396
8 to 13 Hz (alpha) 0.13690234 0.13690234 0.13690234
30 to 50 Hz (gamma) 0.00000038 0.00000038 0.00000038
100 to 128 Hz (reaches Nyquist) 3.52343750 1.02343750 1.02343750

Bands that stop short of Nyquist are untouched. The band that reaches it was 3.44x too large.

For the hang, get_psd_welch(data_len=256, nfft=64, overlap=64, ...) under timeout 15 was killed at 15016 ms with exit 124 before, and returns 13 (INVALID_ARGUMENTS_ERROR) in 9 ms after. Overlaps 0, 8, 16, 32 and 63 all still return STATUS_OK with unchanged values.

All 9 signal processing examples that CI runs pass on this branch. clang-format reports no changes on the touched files, and dotnet build is clean at 0 warnings and 0 errors. The Java change is javadoc only and was not compiled since there is no mvn here, but I checked it with javac -Xdoclint:html after rewording a bare < that would otherwise trip doclint.

Disclosure: this change was prepared with AI assistance. I have reviewed and tested it.

…== nfft

Two independent defects in the same function.

The accumulation loop sums nfft / 2 + 1 bins but the averaging loop only
divides nfft / 2 of them, so the Nyquist bin at index nfft / 2 is summed
across every segment and never divided by the segment count. Its value is
inflated by exactly the number of segments, which means it grows with the
length of the input for a stationary signal. get_band_power inherits the
error for any band that reaches the Nyquist frequency.

Argument validation accepts overlap == nfft. The segment loop then advances
by nfft - overlap, which is zero, so pos never moves and the function spins
forever. Rejecting overlap >= nfft turns a hang into INVALID_ARGUMENTS_ERROR.

Also updates the Python, C# and Java docstrings, which described the valid
range as 0 to nfft inclusive.

Copy link
Copy Markdown
Member

Both native changes look correct, and the reference comparison convincingly isolates the unaveraged Nyquist bin.

Please add automated coverage before merge:

  • assert that overlap == nfft returns INVALID_ARGUMENTS_ERROR without hanging;
  • for multiple segments, assert that the Nyquist bin equals the average of the corresponding per-segment PSD bins.

These are compact regressions for a process hang and silent numerical error, so the example-script coverage alone is not sufficient. It would also be worth checking that all public-language documentation consistently describes the valid range as 0 <= overlap < nfft.

Asserts overlap equal to nfft returns INVALID_ARGUMENTS_ERROR, and runs
that call on a worker thread so a regression of the guard fails rather
than hanging CI. Also asserts every Welch bin including Nyquist equals
the mean of the per-segment PSD bins.
@adityasingh2400

Copy link
Copy Markdown
Author

Added python_package/examples/tests/psd_welch_nyquist_overlap.py and wired it into run_unix.yml. Both checks you asked for are in.

One thing worth calling out. The overlap == nfft call runs on a worker thread with a join timeout, not inline. Inline, a regression of that guard would hang the CI job forever rather than fail, which I confirmed: on master the test process had to be killed at 200 seconds. With the watchdog it fails in 60 with get_psd_welch did not terminate with overlap equal to nfft. ctypes releases the GIL, so the join returns while the native call spins.

The Nyquist check uses no overlap so the segments are disjoint, which lets every Welch bin be compared against the mean of the per-segment get_psd bins rather than only the last one.

On the docs: python, C#, java and swift all state 0 <= overlap < nfft after this change, and swift additionally guards it at the call site. nodejs, julia, matlab, rust, R and cpp do not document the range at all, so there is nothing inconsistent to correct there, only absent. Say the word if you want it added to those too.

@adityasingh2400

Copy link
Copy Markdown
Author

Follow-up: I had only wired the new test into run_unix.yml. The existing python DSP tests also run in run_windows.yml, so I added it there too, matching the shell: cmd style used by the downsampling.py step. run_alpine.yml does not run any python example tests, so I left it alone.

Both workflow files still parse.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants