% Minimal MATLAB-compatible scaffold for reproducing PPG respiratory-rate analysis.
% Paper family: Pimentel et al. / Charlton et al., robust respiratory-rate
% estimation from pulse oximeters / systematic assessment of ECG and PPG RR algorithms.
%
% Dataset: BIDMC PPG and Respiration Dataset from PhysioNet.
% Expected files under ./data:
%   bidmc_##_Signals.csv, bidmc_##_Breaths.csv
%
% Note: the current OpenClaw runtime has no MATLAB/Octave binary installed, so the
% numerically verified companion implementation is run_ppg_rr_reproduction.py.
% This script mirrors the same pipeline in MATLAB style for use in MATLAB.

clear; clc;
root = fileparts(mfilename('fullpath'));
data_dir = fullfile(root, 'data');
fs = 125;
win_sec = 64;
step_sec = 32;
rows = [];

for subj = 1:10
    sig_file = fullfile(data_dir, sprintf('bidmc_%02d_Signals.csv', subj));
    br_file  = fullfile(data_dir, sprintf('bidmc_%02d_Breaths.csv', subj));
    if ~exist(sig_file, 'file') || ~exist(br_file, 'file'), continue; end
    T = readtable(sig_file, 'VariableNamingRule', 'preserve');
    t = T{:,1};
    ppg = T.PLETH;
    resp = T.RESP;
    B = readtable(br_file, 'VariableNamingRule', 'preserve');
    breaths = mean(table2array(B(:,1:min(2,width(B)))), 2, 'omitnan') ./ fs;
    start = 0;
    while start + win_sec <= t(end)
        stop = start + win_sec;
        r = estimate_window(t, ppg, resp, breaths, start, stop, fs);
        if ~isempty(r)
            r.subject = subj;
            rows = [rows; struct2table(r)]; %#ok<AGROW>
        end
        start = start + step_sec;
    end
end

writetable(rows, fullfile(root, 'window_results_matlab.csv'));
summary_methods = {'ppg_raw_bw_rr','ppg_rr_fused','ip_rr','ppg_BW_trough_baseline','ppg_AM_pulse_amplitude','ppg_FM_inter_pulse_interval'};
fprintf('Method\tN\tBias\tMAE\tLoA low\tLoA high\n');
for i = 1:numel(summary_methods)
    m = summary_methods{i};
    if ~ismember(m, rows.Properties.VariableNames), continue; end
    err = rows.(m) - rows.ref_rr;
    err = err(isfinite(err));
    if isempty(err), continue; end
    bias = mean(err);
    sd = std(err);
    fprintf('%s\t%d\t%.3f\t%.3f\t%.3f\t%.3f\n', m, numel(err), bias, mean(abs(err)), bias-1.96*sd, bias+1.96*sd);
end

function r = estimate_window(t, ppg, resp, breaths, start, stop, fs)
    i0 = max(1, floor(start*fs)+1); i1 = min(numel(ppg), floor(stop*fs));
    if i1-i0 < 20*fs, r = []; return; end
    bs = breaths(breaths >= start & breaths <= stop);
    if numel(bs) < 2, r = []; return; end
    ref_rr = 60*(numel(bs)-1)/(bs(end)-bs(1));
    resp_w = zscore_local(detrend_ma(resp(i0:i1), round(8*fs)));
    [ip_rr, ~] = dft_peak_bpm(resp_w, fs);
    ppg_raw_bw = zscore_local(detrend_ma(ppg(i0:i1), round(8*fs)));
    [ppg_raw_bw_rr, ~] = dft_peak_bpm(ppg_raw_bw, fs, 4, 30);
    streams = ppg_features(ppg(i0:i1), fs, start);
    names = fieldnames(streams);
    vals = []; qvals = []; estimates = struct();
    for k = 1:numel(names)
        s = streams.(names{k});
        [grid, y] = resample_linear(s.t, s.v, 5, start, stop);
        y = zscore_local(detrend_ma(y(:), 40));
        [rr, q] = dft_peak_bpm(y, 5);
        estimates.(['ppg_' names{k}]) = rr;
        vals(end+1) = rr; %#ok<AGROW>
        qvals(end+1) = q; %#ok<AGROW>
    end
    if isempty(vals)
        fused = NaN;
    elseif numel(vals) >= 2 && std(vals) <= 4
        fused = mean(vals);
    else
        [~, idx] = max(qvals); fused = vals(idx);
    end
    r = struct('start_s', start, 'stop_s', stop, 'ref_rr', ref_rr, 'ip_rr', ip_rr, 'ppg_raw_bw_rr', ppg_raw_bw_rr, 'ppg_rr_fused', fused);
    for k = 1:numel(names)
        r.(['ppg_' names{k}]) = estimates.(['ppg_' names{k}]);
    end
end

function streams = ppg_features(ppg, fs, start)
    x = zscore_local(detrend_ma(ppg(:), round(1.5*fs)));
    peaks = find_peaks_simple(x, fs);
    troughs = [];
    for i = 1:numel(peaks)-1
        [~, j] = min(x(peaks(i):peaks(i+1)));
        troughs(end+1) = peaks(i)+j-1; %#ok<AGROW>
    end
    streams = struct();
    if numel(troughs) > 5
        streams.BW_trough_baseline = struct('t', start + troughs(:)/fs, 'v', x(troughs));
    end
    amp_t = []; amp_v = [];
    for i = 1:numel(peaks)
        prev = troughs(troughs < peaks(i));
        if ~isempty(prev)
            tr = prev(end);
            amp_t(end+1) = start + peaks(i)/fs; %#ok<AGROW>
            amp_v(end+1) = x(peaks(i))-x(tr); %#ok<AGROW>
        end
    end
    if numel(amp_t) > 5
        streams.AM_pulse_amplitude = struct('t', amp_t(:), 'v', amp_v(:));
    end
    if numel(peaks) > 6
        streams.FM_inter_pulse_interval = struct('t', start + peaks(2:end)'/fs, 'v', diff(peaks(:))/fs);
    end
end

function peaks = find_peaks_simple(x, fs)
    y = movmean(x, max(1, round(0.08*fs)));
    cand = find(y(2:end-1) > y(1:end-2) & y(2:end-1) >= y(3:end)) + 1;
    cand = cand(y(cand) >= mean(y) + 0.10*std(y));
    min_dist = round(0.30*fs);
    peaks = [];
    for i = 1:numel(cand)
        if isempty(peaks) || cand(i)-peaks(end) >= min_dist
            peaks(end+1) = cand(i); %#ok<AGROW>
        elseif y(cand(i)) > y(peaks(end))
            peaks(end) = cand(i);
        end
    end
end

function y = detrend_ma(x, n)
    y = x(:) - movmean(x(:), max(3,n));
end

function y = zscore_local(x)
    y = (x(:)-mean(x(:))) ./ max(std(x(:)), eps);
end

function [grid, out] = resample_linear(t, v, fs_new, start, stop)
    grid = (start:1/fs_new:stop)';
    out = interp1(t(:), v(:), grid, 'linear', 'extrap');
end

function [bpm, q] = dft_peak_bpm(x, fs, lo_bpm, hi_bpm)
    if nargin < 3, lo_bpm = 4; end
    if nargin < 4, hi_bpm = 60; end
    x = x(:) - mean(x(:)); n = numel(x);
    x = x .* hann(n);
    lo = max(1, ceil((lo_bpm/60)*n/fs)); hi = min(floor(n/2), floor((hi_bpm/60)*n/fs));
    pows = zeros(hi-lo+1,1);
    for k = lo:hi
        basis = exp(-1j*2*pi*k*(0:n-1)'/n);
        pows(k-lo+1) = abs(sum(x.*basis)).^2;
    end
    [mx, idx] = max(pows);
    k = lo+idx-1;
    bpm = k*fs/n*60;
    q = mx/(median(pows)+eps);
end
