@@ -15,6 +15,11 @@
 %   saturation - Raw intensity at floor (0) or ceiling (device max) (raw)
 %   sci        - Scalp coupling index via cardiac cross-correlation (raw)
 %   cardiac    - Cardiac peak presence in power spectrum (raw)
+%
+%   sci and cardiac are skipped (channels not penalized, with a
+%   skipReason) when the sampling rate cannot reach the cardiac band or
+%   the recording spans fewer than 5 cycles of CardiacBand(1) (10 s for
+%   the default band). takizawa likewise skips recordings under 10 s.
 %   cov        - Coefficient of variation of raw signal (raw)
 %   takizawa   - 4-rule hemoglobin QC (filtered Hb)
 %
@@ -79,18 +84,21 @@
 nRawCols = size(data.raw, 2);
 
 %% Resolve channel layout from Device or fallback
-hasDev = isfield(data, 'device') && ~isempty(data.device);
+hasDev = isfield(data, 'device') && isa(data.device, 'pf2.Device');
 if hasDev
     dev = data.device;
     nChannels = dev.nChannels;
     wavelengths = dev.wavelengthSet;
     nWl = numel(wavelengths);
-    % Build per-channel column map from Device
+    % Build per-channel column map from Device. Channel k is the k-th
+    % entry of the channel list (as in processFNIRS2), which differs from
+    % the raw channel number on merged probes with non-contiguous numbering
     chNums = dev.channelNumbers();
     wlVec = dev.wavelengths();
+    chList = dev.channelList();
     chColMapWl = cell(1, nChannels); % columns per wavelength (excl dark)
     for ch = 1:nChannels
-        wlCols = find(chNums == ch & wlVec > 0);
+        wlCols = find(chNums == chList(ch) & wlVec > 0);
         chColMapWl{ch} = wlCols;
     end
 elseif isfield(data, 'info') && isfield(data.info, 'synthetic') ...
@@ -331,17 +339,30 @@
             % Need Nyquist > ~1.0 Hz (fs > ~2 Hz) for reliable detection.
             nyq = fs / 2;
             minUsableCardiacHz = 1.0;
+            % A spectral peak needs several beats: at least 5 cycles of
+            % the band's lower edge (10 s for the default band), the same
+            % minimum pf2.qc.sci and pf2.qc.powerSpectrum apply
+            minCardiacDuration = 5 / opts.CardiacBand(1);
+            durationSec = size(data.raw, 1) / fs;
+            cardiacSkipReason = '';
             if nyq <= minUsableCardiacHz
-                % Cannot meaningfully detect cardiac peaks at this fs
+                cardiacSkipReason = sprintf(...
+                    'Sampling rate (%.1f Hz) too low for cardiac peak detection (Nyquist=%.1f Hz, need >%.1f Hz)', ...
+                    fs, nyq, minUsableCardiacHz);
+            elseif durationSec < minCardiacDuration
+                cardiacSkipReason = sprintf(...
+                    'Recording too short (%.1f s < %.0f s) for cardiac peak detection', ...
+                    durationSec, minCardiacDuration);
+            end
+            if ~isempty(cardiacSkipReason)
+                % Cannot meaningfully detect cardiac peaks
                 report.cardiac.detected = false(1, nChannels);
                 report.cardiac.snr = nan(1, nChannels);
                 report.cardiac.freq = nan(1, nChannels);
                 report.cardiac.pass = true(1, nChannels);
                 report.cardiac.threshold = opts.CardiacSNR;
                 report.cardiac.skipped = true;
-                report.cardiac.skipReason = sprintf(...
-                    'Sampling rate (%.1f Hz) too low for cardiac peak detection (Nyquist=%.1f Hz, need >%.1f Hz)', ...
-                    fs, nyq, minUsableCardiacHz);
+                report.cardiac.skipReason = cardiacSkipReason;
                 passMatrix(ci, :) = true;
             else
                 psdResult = pf2.qc.powerSpectrum(data, ...
@@ -408,6 +429,9 @@
             report.takizawa.ruleNames = tkReport.ruleNames;
             report.takizawa.pass = tkReport.pass;
             report.takizawa.skipped = tkReport.skipped;
+            if tkReport.skipped
+                report.takizawa.skipReason = tkReport.skipReason;
+            end
             passMatrix(ci, :) = tkReport.pass;
     end
 end
 
@@ -43,12 +43,13 @@ function report(qcReport)
     checkName = checks{ci};
     if isfield(qcReport, checkName) && isfield(qcReport.(checkName), 'skipped') ...
             && qcReport.(checkName).skipped
-        reason = '';
-        if isfield(qcReport.(checkName), 'skipReason')
-            reason = [' — ' qcReport.(checkName).skipReason];
+        % The reason says why (low sampling rate, short recording, ...)
+        reason = ' (not applicable to this recording)';
+        if isfield(qcReport.(checkName), 'skipReason') ...
+                && ~isempty(qcReport.(checkName).skipReason)
+            reason = [': ' qcReport.(checkName).skipReason];
         end
-        fprintf('  Note: %s check SKIPPED (not applicable at this fs)%s\n', ...
-            capitalize(checkName), reason);
+        fprintf('  Note: %s check SKIPPED%s\n', capitalize(checkName), reason);
     end
 end
 fprintf('\n');
 
@@ -152,6 +152,11 @@
 
 nCh = numel(channels);
 layout = lower(char(opts.Layout));
+if nCh == 0
+    % No channels (e.g. all masked bad): tiled has no subplots to title,
+    % so draw a single empty axes instead
+    layout = 'overlay';
+end
 
 % Physiological bands for shading
 bands = struct();
@@ -244,10 +249,13 @@ function drawBandShading(ax, freqs, psdMatrix, bands, bandColors, bandLabels)
 
 fMin = min(freqs);
 fMax = max(freqs);
-yLimits = [min(psdMatrix(psdMatrix > 0), [], 'all') * 0.1, ...
-           max(psdMatrix(:), [], 'all') * 10];
-if isempty(yLimits) || any(isnan(yLimits)) || yLimits(1) >= yLimits(2)
+% An all-zero (or empty) spectrum, e.g. from a recording too short for a
+% PSD, has no positive power to scale a log axis to
+posPower = psdMatrix(psdMatrix > 0 & isfinite(psdMatrix));
+if isempty(posPower)
     yLimits = [1e-10, 1];
+else
+    yLimits = [min(posPower) * 0.1, max(posPower) * 10];
 end
 
 bandNames = {'mayer', 'respiratory', 'cardiac'};
 
@@ -19,7 +19,8 @@
 %                    (default: 'auto')
 %   'Channels'     - [1 x C] channel indices to analyze. Default: all good
 %                    channels from fchMask. (default: [])
-%   'WindowLength' - PSD window length in seconds (default: 10)
+%   'WindowLength' - PSD window length in seconds; clamped to between
+%                    3 samples and the recording length (default: 10)
 %   'Overlap'      - Fractional overlap between windows, 0 to <1
 %                    (default: 0.5)
 %   'FreqRange'    - [1x2] frequency range to return in Hz. Upper bound is
@@ -28,7 +29,11 @@
 %
 % Inputs:
 %   data - fNIRS data struct. Must contain .fs and the selected signal
-%          field (.raw, .HbO, .HbR, or .OD).
+%          field (.raw, .HbO, .HbR, or .OD). For raw data, an optional
+%          .device (pf2.Device) selects each channel's first non-dark
+%          wavelength column, with channel k taken as the k-th entry of
+%          the device channel list; without it, columns are assumed
+%          interleaved.
 %
 % Outputs:
 %   result - Struct with fields:
@@ -39,9 +44,16 @@
 %            .fs        - Sampling rate
 %
 %            When DetectPeaks = true, also includes:
-%            .cardiac      - struct with .freq, .power, .snr, .detected
-%            .respiratory  - struct with .freq, .power, .detected
-%            .mayer        - struct with .freq, .power, .detected
+%            .cardiac      - struct with .freq, .power, .snr, .detected,
+%                            .skipped, .skipReason
+%            .respiratory  - struct with .freq, .power, .detected,
+%                            .skipped, .skipReason
+%            .mayer        - struct with .freq, .power, .detected,
+%                            .skipped, .skipReason
+%            A band is skipped (detected = false, freq/power NaN) when
+%            the recording spans fewer than 5 cycles of the band's lower
+%            edge: under 10 s for cardiac, 50 s for respiratory, 100 s
+%            for Mayer.
 %
 % Algorithm:
 %   1. Select signal from data struct based on Signal parameter
@@ -108,22 +120,39 @@
 
 %% Select columns from signal data
 if strcmp(signalType, 'raw')
-    % For raw data, resolve which columns to use per channel
-    % Use first wavelength of each channel pair (alternating layout)
+    % For raw data, resolve the first wavelength column of each channel
     nRawCols = size(sigData, 2);
-    if isfield(data, 'info') && isfield(data.info, 'synthetic') ...
-            && isfield(data.info.synthetic, 'wavelengths')
-        nWl = numel(data.info.synthetic.wavelengths);
-    elseif isfield(data, 'probeinfo')
-        nWl = 2;  % Default assumption
+    if isfield(data, 'device') && isa(data.device, 'pf2.Device')
+        % Device layout: skip dark columns (wavelength 0) and handle
+        % non-interleaved column orders. Channel k is the k-th entry of
+        % the device channel list (as in processFNIRS2), so non-contiguous
+        % numbering such as merged probes still maps by position.
+        chNums = data.device.channelNumbers();
+        wlVec = data.device.wavelengths();
+        chList = data.device.channelList();
+        colIndices = nan(1, numel(channels));
+        for k = 1:numel(channels)
+            if channels(k) < 1 || channels(k) > numel(chList)
+                continue;
+            end
+            c = find(chNums == chList(channels(k)) & wlVec > 0, 1);
+            if ~isempty(c) && c <= nRawCols
+                colIndices(k) = c;
+            end
+        end
     else
-        nWl = 2;  % Default: 2 wavelengths
+        if isfield(data, 'info') && isfield(data.info, 'synthetic') ...
+                && isfield(data.info.synthetic, 'wavelengths')
+            nWl = numel(data.info.synthetic.wavelengths);
+        else
+            nWl = 2;  % Default: 2 interleaved wavelengths
+        end
+        colIndices = (channels - 1) * nWl + 1;
+        colIndices(colIndices > nRawCols) = NaN;
     end
-    % Map channel index to first wavelength column
-    colIndices = (channels - 1) * nWl + 1;
-    colIndices(colIndices > nRawCols) = [];
-    channels = channels(1:numel(colIndices));
-    sigMatrix = sigData(:, colIndices);
+    keep = ~isnan(colIndices);
+    channels = channels(keep);
+    sigMatrix = sigData(:, colIndices(keep));
 else
     % HbO, HbR, OD — columns map directly to channels
     validCh = channels(channels <= size(sigData, 2));
@@ -139,10 +168,29 @@
 freqRange(2) = min(freqRange(2), nyquist);
 
 %% Compute PSD using Welch's method
-windowSamples = round(opts.WindowLength * fs);
-overlapSamples = round(windowSamples * opts.Overlap);
+% Clamp the window to [3, recording length]: pwelch rejects windows longer
+% than the signal, and the welchPSD Hann window has zero power under 3
+nSamples = size(sigMatrix, 1);
+windowSamples = min(max(round(opts.WindowLength * fs), 3), nSamples);
+overlapSamples = min(round(windowSamples * opts.Overlap), windowSamples - 1);
 nfft = max(256, 2^nextpow2(windowSamples));
 
+if nChannels == 0 || nSamples < 3
+    % Nothing to analyze (e.g. every channel masked bad): return empty
+    % per-channel results instead of indexing into an empty matrix
+    f = (0:nfft/2)' * fs / nfft;
+    f = f(f >= freqRange(1) & f <= freqRange(2));
+    result.psd = zeros(numel(f), nChannels);
+    result.freqs = f;
+    result.channels = channels;
+    result.signal = signalType;
+    result.fs = fs;
+    if opts.DetectPeaks
+        result = detectAllBands(result, nSamples / fs);
+    end
+    return;
+end
+
 % Check for pwelch availability
 hasPwelch = ~isempty(which('pwelch'));
 
@@ -182,9 +230,7 @@
 
 %% Detect physiological peaks
 if opts.DetectPeaks
-    result.cardiac = detectBandPeak(f, psdMatrix, [0.5, 2.5], 'cardiac');
-    result.respiratory = detectBandPeak(f, psdMatrix, [0.1, 0.5], 'respiratory');
-    result.mayer = detectBandPeak(f, psdMatrix, [0.05, 0.15], 'mayer');
+    result = detectAllBands(result, nSamples / fs);
 end
 
 end
@@ -292,6 +338,37 @@
 end
 
 
+function result = detectAllBands(result, durationSec)
+% DETECTALLBANDS Run peak detection for each physiological band
+%
+% A band needs the recording to span at least minCycles cycles of its
+% lower edge; shorter recordings cannot resolve a peak there, so the band
+% is reported as skipped rather than as a spurious detection or miss.
+
+minCycles = 5;
+bandDefs = {'cardiac', [0.5, 2.5]; 'respiratory', [0.1, 0.5]; 'mayer', [0.05, 0.15]};
+
+for b = 1:size(bandDefs, 1)
+    name = bandDefs{b, 1};
+    limits = bandDefs{b, 2};
+    minDuration = minCycles / limits(1);
+    if durationSec < minDuration
+        bandResult = detectBandPeak(result.freqs, zeros(size(result.psd)), limits, name);
+        bandResult.skipped = true;
+        bandResult.skipReason = sprintf( ...
+            'Recording too short (%.1f s < %.0f s) for %s band [%.2f, %.2f] Hz', ...
+            durationSec, minDuration, name, limits(1), limits(2));
+    else
+        bandResult = detectBandPeak(result.freqs, result.psd, limits, name);
+        bandResult.skipped = false;
+        bandResult.skipReason = '';
+    end
+    result.(name) = bandResult;
+end
+
+end
+
+
 function bandResult = detectBandPeak(freqs, psdMatrix, bandLimits, bandName)
 % DETECTBANDPEAK Find peak in a physiological frequency band
 %
 
@@ -36,9 +36,11 @@
 %          .raw  - [T x C_raw] raw light intensity
 %          .fs   - Sampling frequency (Hz)
 %          Wavelength layout resolved from (in order):
-%            1. Explicit Wavelengths/ChannelNumbers parameters
-%            2. data.probeinfo.Probe{1}.TableCh
-%            3. data.info.synthetic.wavelengths (alternating ch1_wl1, ch1_wl2, ...)
+%            1. data.device (pf2.Device); channel i is the i-th entry of
+%               the device channel list
+%            2. Explicit Wavelengths/ChannelNumbers parameters
+%            3. data.probeinfo.Probe{1}.TableCh
+%            4. data.info.synthetic.wavelengths (alternating ch1_wl1, ch1_wl2, ...)
 %
 % Outputs:
 %   result - Struct with fields:
@@ -48,6 +50,12 @@
 %            .threshold - Threshold used
 %            .cardiacBand - [low, high] Hz
 %            .fs        - Sampling rate
+%            .skipped   - true when SCI could not be computed (sampling
+%                         rate below the cardiac band, or a recording
+%                         shorter than 5 cycles of CardiacBand(1), i.e.
+%                         10 s for the default band); .sci is then NaN and
+%                         .isGood true so channels are not penalized
+%            .skipReason - Reason for the skip (only when skipped)
 %
 % Algorithm:
 %   1. Resolve wavelength layout from data or explicit parameters
@@ -124,6 +132,28 @@
     cardiacBand(2) = clampedHigh;
 end
 
+%% Check recording length: the cardiac correlation needs several beats
+% (at least 5 cycles of the band's lower edge, 10 s for the default band),
+% and the bandpass filter needs more samples than its padding length
+minCycles = 5;
+minDuration = minCycles / cardiacBand(1);
+nSamples = size(data.raw, 1);
+durationSec = nSamples / fs;
+minSamples = 3 * (2 * opts.FilterOrder) + 1;
+if durationSec < minDuration || nSamples < minSamples
+    result.sci = nan(1, nChannels);
+    result.isGood = true(1, nChannels);
+    result.channels = 1:nChannels;
+    result.threshold = opts.Threshold;
+    result.cardiacBand = cardiacBand;
+    result.fs = fs;
+    result.skipped = true;
+    result.skipReason = sprintf( ...
+        'Recording too short (%.1f s < %.0f s) for cardiac band [%.1f, %.1f] Hz', ...
+        durationSec, minDuration, cardiacBand(1), cardiacBand(2));
+    return;
+end
+
 %% Compute SCI for each channel
 sciValues = zeros(1, nChannels);
 
@@ -187,12 +217,15 @@
     wl = data.device.wavelengths();
     chNums = data.device.channelNumbers();
 
-    uniqueCh = unique(chNums(chNums > 0 & ~isnan(chNums)));
-    nChannels = numel(uniqueCh);
+    % Channel i is the i-th entry of the device channel list (as in
+    % processFNIRS2), which can differ from the unique raw channel numbers
+    % on merged probes
+    chList = data.device.channelList();
+    nChannels = numel(chList);
     channelMap = zeros(nChannels, 2);
 
     for i = 1:nChannels
-        cols = find(chNums == uniqueCh(i) & wl > 0 & ~isnan(wl));
+        cols = find(chNums == chList(i) & wl > 0 & ~isnan(wl));
         if numel(cols) >= 2
             channelMap(i, :) = cols(1:2);
         end