Skip to content

Commit b9c2aa8

Browse files
AdrianCurtinAdrianCurtin
committed
Return NaN from wcoherence for missing or constant input
With cached CWTs, the coefficients of a series containing NaN are NaN, and the clamp min(wcoh, 1) turned them into a coherence of exactly 1 because min ignores NaN. Hyperscanning analyses that cache per-channel CWTs and include a fully rejected channel or ROI were inflated toward 1 for every pair involving it. Computed internally, an all-NaN series was zero-filled and gave 0. wcoherence now returns NaN when either input has fewer than two finite samples or is constant, and the clamp leaves NaN in place. Add WaveletCoherenceTest, which covers both paths, constant input, a short gap and the [0, 1] bounds; the missing- and constant-input tests fail without the fix.
1 parent 599df3a commit b9c2aa8

3 files changed

Lines changed: 109 additions & 2 deletions

File tree

‎+pf2_base/+tests/+unit/WaveletCoherenceTest.m‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
classdef WaveletCoherenceTest < matlab.unittest.TestCase
2+
% WAVELETCOHERENCETEST Tests for pf2_base.wavelet.wcoherence
3+
%
4+
% Covers:
5+
% - bounds: identical signals give coherence near 1, values lie in [0, 1]
6+
% - missing input: an all-NaN series gives NaN, both when the CWTs are
7+
% computed internally and when cached CWTs are supplied (the cached
8+
% coefficients of a NaN series are NaN, and the clamp to 1 previously
9+
% turned them into perfect coherence)
10+
% - constant input gives NaN
11+
% - a series with a short gap stays finite and in range
12+
%
13+
% Example:
14+
% results = runtests('pf2_base.tests.unit.WaveletCoherenceTest');
15+
16+
properties
17+
fs = 2
18+
n = 1200
19+
band = [0.009 0.06]
20+
end
21+
22+
methods (Test)
23+
24+
function testIdenticalSignalsNearOne(tc)
25+
rng(1);
26+
x = cumsum(randn(tc.n, 1));
27+
r = pf2_base.wavelet.wcoherence(x, x, tc.fs, 'FreqRange', tc.band);
28+
tc.verifyGreaterThan(r.value, 0.99);
29+
tc.verifyLessThanOrEqual(r.value, 1);
30+
end
31+
32+
function testValuesWithinUnitInterval(tc)
33+
rng(2);
34+
r = pf2_base.wavelet.wcoherence(randn(tc.n, 1), randn(tc.n, 1), tc.fs, ...
35+
'FreqRange', tc.band);
36+
w = r.wcoh(isfinite(r.wcoh));
37+
tc.verifyGreaterThanOrEqual(min(w), 0);
38+
tc.verifyLessThanOrEqual(max(w), 1);
39+
end
40+
41+
function testAllNaNInputGivesNaN(tc)
42+
rng(3);
43+
r = pf2_base.wavelet.wcoherence(randn(tc.n, 1), nan(tc.n, 1), tc.fs, ...
44+
'FreqRange', tc.band);
45+
tc.verifyTrue(isnan(r.value));
46+
end
47+
48+
function testCachedCwtOfNaNSeriesGivesNaN(tc)
49+
rng(4);
50+
X = randn(tc.n, 2);
51+
Y = randn(tc.n, 2);
52+
Y(:, 2) = NaN;
53+
args = {'FreqRange', [0.0045 0.12], 'VoicesPerOctave', 8};
54+
Wx = pf2_base.wavelet.cwt(X, tc.fs, args{:});
55+
Wy = pf2_base.wavelet.cwt(Y, tc.fs, args{:});
56+
slice = @(W, c) struct('coeffs', W.coeffs(:, :, c), 'freqs', W.freqs, ...
57+
'coi', W.coi, 'scales', W.scales);
58+
59+
rGood = pf2_base.wavelet.wcoherence(X(:, 1), Y(:, 1), tc.fs, ...
60+
'FreqRange', tc.band, 'CwtX', slice(Wx, 1), 'CwtY', slice(Wy, 1));
61+
rNaN = pf2_base.wavelet.wcoherence(X(:, 1), Y(:, 2), tc.fs, ...
62+
'FreqRange', tc.band, 'CwtX', slice(Wx, 1), 'CwtY', slice(Wy, 2));
63+
64+
tc.verifyTrue(isfinite(rGood.value));
65+
tc.verifyLessThan(rGood.value, 1);
66+
tc.verifyTrue(isnan(rNaN.value));
67+
end
68+
69+
function testConstantInputGivesNaN(tc)
70+
rng(5);
71+
r = pf2_base.wavelet.wcoherence(randn(tc.n, 1), zeros(tc.n, 1), tc.fs, ...
72+
'FreqRange', tc.band);
73+
tc.verifyTrue(isnan(r.value));
74+
end
75+
76+
function testShortGapStaysFinite(tc)
77+
rng(6);
78+
y = randn(tc.n, 1);
79+
y(100:160) = NaN;
80+
r = pf2_base.wavelet.wcoherence(randn(tc.n, 1), y, tc.fs, 'FreqRange', tc.band);
81+
tc.verifyTrue(isfinite(r.value));
82+
tc.verifyGreaterThanOrEqual(r.value, 0);
83+
tc.verifyLessThanOrEqual(r.value, 1);
84+
end
85+
86+
end
87+
end

‎+pf2_base/+wavelet/wcoherence.m‎

Lines changed: 21 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,10 @@
3737
%
3838
% Outputs:
3939
% result - Struct with fields:
40-
% .value - Mean WCT magnitude in FreqRange (scalar, COI-masked)
40+
% .value - Mean WCT magnitude in FreqRange (scalar, COI-masked).
41+
% NaN when either input has no finite samples or is
42+
% constant, or when the supplied CWT coefficients are NaN;
43+
% coherence is undefined for those inputs.
4144
% .pvalue - NaN (use permutation test for significance)
4245
% .method - 'wcoherence'
4346
% .windowed - false
@@ -97,6 +100,9 @@
97100
y = y(:);
98101
T = length(x);
99102

103+
% Coherence is undefined for a missing or constant input
104+
degenerate = isDegenerate(x) || isDegenerate(y);
105+
100106
% NaN handling: linear interpolation
101107
x = fillNaN(x);
102108
y = fillNaN(y);
@@ -155,7 +161,13 @@
155161

156162
% --- Wavelet coherence ---
157163
wcoh = abs(sWxy).^2 ./ (sWxx .* sWyy + eps(class(sWxx)));
158-
wcoh = min(wcoh, 1);
164+
% Clamp rounding overshoot; min(NaN, 1) would return 1 and turn
165+
% NaN coefficients (e.g. a cached CWT of a missing series) into
166+
% perfect coherence
167+
wcoh(wcoh > 1) = 1;
168+
if degenerate
169+
wcoh(:) = NaN;
170+
end
159171

160172
% --- Build time vector ---
161173
times = (0:T-1)' / fs;
@@ -251,6 +263,13 @@
251263
end
252264

253265

266+
function tf = isDegenerate(v)
267+
% True when v has fewer than two finite samples or no variation
268+
f = v(isfinite(v));
269+
tf = numel(f) < 2 || max(f) == min(f);
270+
end
271+
272+
254273
function v = fillNaN(v)
255274
% Linear interpolation of NaN values
256275
nanIdx = isnan(v);

‎CHANGELOG.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ Groundwork for running headless single-subject processing under GNU Octave.
2929
- `pf2_base.normalizeMarkerDict` parses `{code, label}` cell dictionaries element by element, so numeric codes keep one row each under Octave.
3030
- `pf2_base.external.zp2sos` orders complex roots with `sortrows` instead of `sort(..., 'ComparisonMethod', 'real')`, which Octave lacks; the resulting second-order sections are identical.
3131
- `pf2.data.split` and `pf2.data.resample` failed on aux signals stored as tables with their own time column (for example a derived HRV series): MATLAB R2025b's `table` rejects a single char name for `'VariableNames'`. The time-column name is now passed as a cell array. v1.1.0 is affected the same way.
32+
- `pf2_base.wavelet.wcoherence` returns NaN when either input has no finite samples or is constant. With cached CWTs (`'CwtX'`/`'CwtY'`) the coefficients of a series containing NaN are NaN, and the clamp `min(wcoh, 1)` turned them into a coherence of exactly 1, because `min` ignores NaN; computed internally, an all-NaN series was zero-filled and gave 0. The clamp now leaves NaN in place. Hyperscanning analyses that cache per-channel CWTs and include a fully rejected channel or ROI were inflated toward 1 for every pair involving it. `WaveletCoherenceTest` covers both paths.
3233
- `pf2.import.importHitachiMES` no longer requests the unused second output of `strsplit`, and `pf2.data.split` assigns its time-size variables directly; both failed under Octave.
3334
- `loadDeviceCfg` and `pf2_base.gui.updateCurrentDevice` assign whole table columns instead of growing them by indexed assignment, and no longer call `sprintf` on a categorical or `nansum` on a logical.
3435
- `bvoxy` declares its positional `coefs` and `isOD` options before its name-value parameters. Argument handling under MATLAB is unchanged.

0 commit comments

Comments
 (0)