開始
│
├─ 參數設置
│   ├─ 設置采樣率、持續時間、基礎心率
│   ├─ 配置心率變異性參數
│   └─ 設置噪聲參數
│
├─ 心電信號合成
│   ├─ 生成RR間期序列
│   ├─ 構建心搏模板
│   ├─ 合成干凈ECG信號
│   └─ 添加噪聲干擾
│
├─ 信號預處理
│   ├─ 高通濾波去除基線漂移
│   ├─ 陷波濾波去除電源干擾
│   └─ 帶通濾波增強QRS波
│
├─ QRS波檢測
│   ├─ 差分運算增強高頻成分
│   ├─ 平方運算強化QRS波
│   ├─ 移動積分平滑信號
│   ├─ 自適應閾值檢測候選點
│   └─ 局部最大值精煉R峰位置
│
├─ 心率變異性分析
│   ├─ 計算時域指標
│   ├─ 計算頻域指標
│   ├─ 計算非線性指標
│   └─ 房顫樣不規則性篩查
│
├─ 結果可視化
│   ├─ 繪制原始和濾波后信號
│   ├─ 標記R峰位置
│   ├─ 顯示QRS增強信號
│   ├─ 繪制RR間期序列
│   ├─ 展示功率譜密度
│   └─ 呈現龐加萊圖
│
└─ 結果報告
    ├─ 輸出數值指標
    └─ 提供房顫篩查結果

function ecg_hrv_pipeline()
% ECG + HRV 高級處理流程（生物電工程項目）
% ---------------------------------------------------------------
% 本腳本功能：
% 1) 合成包含心率變異性和常見噪聲的逼真心電圖信號
% 2) 去噪（基線和電源干擾），增強QRS波，檢測R峰
% 3) 計算心率變異的時域、頻域和非線性指標
% 4) 篩查房顫樣不規則性（簡單啟發式方法）
% 5) 可視化每個處理步驟
%
% 運行方法：
%   - 保存為 ecg_hrv_pipeline.m
%   - 在MATLAB中按運行（F5）
%
% 如果信號處理工具箱可用，我們使用butter/filtfilt進行零相位IIR濾波
% 否則，我們回退到簡單的卷積方法


%% -------------------- 參數設置 --------------------
fs      = 360;              % 采樣率（赫茲）
dur_s   = 180;              % 持續時間（秒）
hr_bpm  = 70;               % 標稱心率（每分鐘心跳次數）
af_mode = false;            % 設置為true可生成不規則的房顫樣心律


% 心率變異性（RR間期變異性）設置
sdnn_ms     = 60;           % 目標SDNN（毫秒），用于變異性真實性
lf_hz       = 0.1;          % 低頻調制（約梅耶波）
hf_hz       = 0.25;         % 高頻調制（呼吸相關）


% 噪聲設置
add_baseline_wander = true; % 是否添加基線漂移
baseline_amp_mV     = 0.15; % 基線漂移幅度（毫伏）
baseline_hz         = 0.35; % 基線漂移頻率（赫茲）


add_powerline = true;       % 是否添加電源干擾
powerline_freq = 50;        % 根據地區選擇50或60赫茲
powerline_amp_mV = 0.05;    % 電源干擾幅度（毫伏）


rng(42); % 設置隨機種子以確保結果可重復


%% -------------------- 合成ECG信號 --------------------
N = fs*dur_s; % 計算總樣本數
t = (0:N-1)/fs; % 生成時間向量


% 生成具有LF/HF變異性或房顫樣抖動的RR間期序列
[rr_s, r_times] = generate_rr_series(dur_s, hr_bpm, sdnn_ms, lf_hz, hf_hz, af_mode);


% 通過累加具有形態學特征的心搏（高斯小波）構建合成ECG
ecg_clean = synth_ecg_from_rr(t, r_times, fs);


% 添加逼真噪聲
ecg_noisy = ecg_clean; % 初始化含噪ECG
if add_baseline_wander % 添加基線漂移
    ecg_noisy = ecg_noisy + baseline_amp_mV * sin(2*pi*baseline_hz*t);
end
if add_powerline % 添加電源干擾
    ecg_noisy = ecg_noisy + powerline_amp_mV * sin(2*pi*powerline_freq*t);
end
% 添加白色傳感器噪聲
ecg_noisy = ecg_noisy + 0.02*randn(size(t));


%% -------------------- 去噪處理 --------------------
% 1) 去除基線漂移（高通濾波～0.5 Hz）
ecg_hp = hp_filter(ecg_noisy, fs, 0.5);


% 2) 電源干擾陷波濾波（50/60 Hz）
ecg_notch = notch_filter(ecg_hp, fs, powerline_freq);


% 3) 帶通濾波聚焦QRS波（經典5-15 Hz；我們稍寬一些）
ecg_bp = bp_filter(ecg_notch, fs, 5, 25);


%% -------------------- QRS波增強 + R峰檢測 --------------------
% Pan-Tompkins風格：差分 -> 平方 -> 移動積分
diff_sig = [0 diff(ecg_bp)]; % 一階差分
sq_sig   = diff_sig.^2; % 平方運算增強高頻成分
win_ms   = 120; % 積分窗口約120毫秒
win_n    = max(1, round(fs*win_ms/1000)); % 計算窗口樣本數
int_sig  = movmean(sq_sig, win_n); % 移動平均積分


% 帶不應期的自適應閾值檢測
[locs_R, thr] = detect_r_peaks(int_sig, fs);


% 在原始濾波ECG上精煉R峰位置：在檢測到的索引周圍搜索局部最大值
search_radius = round(0.05*fs); % 搜索半徑（樣本數）
locs_R = refine_r_peaks(ecg_bp, locs_R, search_radius);


%% -------------------- 心率變異性指標計算 --------------------
RR  = diff(locs_R)/fs; % RR間期（秒）
NN  = RR(~isoutlier(RR, 'median')); % 使用中位數法剔除異常值，獲取正常竇性心律間期
if numel(NN) < 3 % 檢查是否有足夠的數據點
    warning('NN間期數量太少——請嘗試更長的持續時間或調整參數');
end


% 時域指標
metrics.meanNN_ms = mean(NN)*1000; % 平均NN間期（毫秒）
metrics.sdnn_ms   = std(NN)*1000; % NN間期標準差（毫秒）
metrics.rmssd_ms  = sqrt(mean(diff(NN).^2))*1000; % 相鄰NN間期差值的均方根（毫秒）
metrics.pnn50     = 100*mean(abs(diff(NN)) > 0.050); % 相鄰NN間期差值大于50毫秒的百分比


% 頻域指標：通過對插值后的心率信號進行Welch譜分析
[lfnu, hfnu, lf_hz_est, hf_hz_est, lf_power, hf_power, total_power] = hrv_freq(NN, locs_R(1:end-1)/fs);


metrics.lfnu = lfnu; metrics.hfnu = hfnu; % 歸一化低頻和高頻功率
metrics.lf_power = lf_power; metrics.hf_power = hf_power; metrics.total_power = total_power; % 絕對功率
metrics.lf_center_hz = lf_hz_est; metrics.hf_center_hz = hf_hz_est; % 頻帶中心頻率


% 非線性指標：龐加萊圖（SD1/SD2）
[SD1, SD2] = poincare_SD1_SD2(NN);
metrics.SD1_ms = SD1*1000; metrics.SD2_ms = SD2*1000; % 龐加萊圖指標


% 簡單房顫樣不規則性篩查（啟發式方法）：
%   - RR間期的高變異系數，低自相關性，升高的pNN50，降低的SD1/SD2比率穩定性
cv_rr   = std(RR)/mean(RR); % RR間期變異系數
ac_rr   = autocorr_stat(RR, 1); % 一階自相關系數
metrics.cv_rr = cv_rr; metrics.ac_lag1 = ac_rr; % 存儲指標


% 房顫篩查條件判斷
af_flag = (cv_rr > 0.15) && (metrics.pnn50 > 15) && (ac_rr < 0.35);
metrics.af_screen_positive = af_flag; % 房顫篩查結果


%% -------------------- 結果報告與可視化 --------------------
figure('Color','w','Name','ECG & HRV處理流程','Position',[80 80 1200 800]); % 創建圖形窗口


% 子圖1：原始和含噪ECG對比
subplot(4,1,1);
plot(t, ecg_noisy, 'Color', [0.2 0.2 0.2]); hold on;
plot(t, ecg_clean, 'r');
legend('含噪ECG','干凈ECG（真實值）'); xlabel('時間（秒）'); ylabel('毫伏');
title('含噪聲的合成ECG信號');


% 子圖2：濾波后ECG和檢測到的R峰
subplot(4,1,2);
plot(t, ecg_notch); hold on;
stem(locs_R/fs, ecg_notch(locs_R), 'g','filled','MarkerSize',3);
xlabel('時間（秒）'); ylabel('毫伏');
title('濾波后ECG + 檢測到的R峰');


% 子圖3：QRS增強信號和閾值
subplot(4,1,3);
plot(t, int_sig); hold on; yline(thr, '--r');
xlabel('時間（秒）'); ylabel('任意單位');
title('QRS增強（積分信號）+ 檢測閾值');


% 子圖4：RR間期序列（心率變異性）
subplot(4,1,4);
plot(locs_R(1:end-1)/fs, RR*1000, '.-'); grid on;
xlabel('時間（秒）'); ylabel('RR間期（毫秒）');
title('RR間期序列（心率變異性）');


% 在命令窗口打印指標
disp('--- 心率變異性指標 ---');
disp(metrics);


% 第二個圖形窗口：頻域分析和非線性分析
figure('Color','w','Name','心率變異性頻域和非線性分析','Position',[120 120 1100 400]);
% 子圖1：功率譜密度
subplot(1,3,1);
[pxx,f] = hrv_psd(NN, locs_R(1:end-1)/fs);
plot(f, pxx); xlim([0 0.5]); grid on;
xlabel('頻率（赫茲）'); ylabel('功率（毫秒2/赫茲）'); title('心率變異性功率譜密度（Welch方法）');


% 子圖2：龐加萊圖
subplot(1,3,2);
NN1 = NN(1:end-1); NN2 = NN(2:end);
plot(NN1*1000, NN2*1000, '.'); grid on; axis equal;
xlabel('NN_n（毫秒）'); ylabel('NN_{n+1}（毫秒）');
title(sprintf('龐加萊圖（SD1=%.1f毫秒, SD2=%.1f毫秒）', metrics.SD1_ms, metrics.SD2_ms));


% 子圖3：歸一化LF/HF功率
subplot(1,3,3);
bar([metrics.lfnu, metrics.hfnu]); set(gca,'XTickLabel',{'LFnu','HFnu'}); ylim([0 100]); grid on;
title('歸一化LF/HF功率（%）');


% 簡單文本面板輸出
fprintf('\n房顫樣不規則性篩查：%s\n', ternary(metrics.af_screen_positive, '陽性（啟發式）', '陰性'));
fprintf('平均NN間期 = %.0f毫秒 | SDNN = %.0f毫秒 | RMSSD = %.0f毫秒 | pNN50 = %.1f%%\n', ...
    metrics.meanNN_ms, metrics.sdnn_ms, metrics.rmssd_ms, metrics.pnn50);
fprintf('LFnu = %.1f%% | HFnu = %.1f%% | RR間期變異系數 = %.2f | 滯后1自相關 = %.2f\n\n', ...
    metrics.lfnu, metrics.hfnu, metrics.cv_rr, metrics.ac_lag1);


end


%% ================== 輔助函數 ==================


function [rr_s, r_times] = generate_rr_series(dur_s, hr_bpm, sdnn_ms, lf_hz, hf_hz, af_mode)
% 生成具有LF/HF調制或房顫樣抖動的RR間期序列
mean_rr = 60/hr_bpm; % 平均RR間期（秒）
tstep = 0.25; % 調制時間步長
tmod = 0:tstep:dur_s; % 調制時間向量


if af_mode % 房顫模式
    % 房顫樣：更大、更隨機的變異性，帶有1/f^α噪聲
    base = mean_rr + 0.30*randn(size(tmod));
    base = max(0.3, base); % 確保RR間期不小于0.3秒
else % 正常模式
    % 準正弦LF/HF調制
    lf = 0.05*sin(2*pi*lf_hz*tmod + 2*pi*rand);
    hf = 0.03*sin(2*pi*hf_hz*tmod + 2*pi*rand);
    base = mean_rr + lf + hf + 0.01*randn(size(tmod));
end


% 縮放到目標SDNN
rr_interp = interp1(tmod, base, 0:mean_rr:dur_s+10, 'linear', 'extrap');
rr_s = rr_interp(:);
scale = (sdnn_ms/1000) / std(rr_s);
rr_s = (rr_s - mean(rr_s))*scale + mean_rr;
rr_s(rr_s < 0.3) = 0.3; % 確保RR間期不小于0.3秒


r_times = cumsum(rr_s); % 累積RR間期得到R峰時間
r_times = r_times(r_times < dur_s); % 截斷到持續時間范圍內
rr_s = diff([0; r_times]); % 從R峰時間重新計算RR間期
end


function ecg = synth_ecg_from_rr(t, r_times, fs)
% 簡單心搏模板：使用高斯波形疊加模擬P、Q、R、S、T波
% 幅度（毫伏）和寬度調整以符合真實情況
ecg = zeros(size(t)); % 初始化ECG信號
% 模板參數相對于R峰時間（秒）
comp = [ ... %   延遲(秒)   幅度(毫伏)  寬度(秒)
          -0.20       0.10     0.05;   % P波
          -0.04      -0.15     0.015;  % Q波
           0.00       1.00     0.012;  % R波
           0.02      -0.25     0.015;  % S波
           0.25       0.35     0.08];  % T波
for k=1:numel(r_times) % 對每個R峰時間
    for c=1:size(comp,1) % 對每個波形成分
        mu = r_times(k) + comp(c,1); % 成分中心時間
        A  = comp(c,2); % 成分幅度
        w  = comp(c,3); % 成分寬度
        ecg = ecg + A * exp(-0.5*((t-mu)/w).^2); % 添加高斯波形
    end
end
% 添加輕微形態變異性
ecg = ecg + 0.005*randn(size(ecg));
end


function y = hp_filter(x, fs, fc)
% 高通濾波～去除基線漂移
if has_spt() % 檢查是否有信號處理工具箱
    [b,a] = butter(2, fc/(fs/2), 'high'); % 設計2階高通巴特沃斯濾波器
    y = filtfilt(b,a,x); % 零相位濾波
else % 無信號處理工具箱的回退方法
    % 通過移動中位數和移動平均進行簡單去趨勢
    w = max(1, round(fs*0.6)); % 窗口大小
    x_med = movmedian(x, w); % 移動中位數
    y = x - movmean(x_med, w); % 減去移動平均
end
end


function y = notch_filter(x, fs, f0)
% 電源干擾陷波濾波
Q = 30; % 品質因數
if has_spt() % 檢查是否有信號處理工具箱
    w0 = f0/(fs/2); % 歸一化中心頻率
    bw = w0/Q; % 帶寬
    [b,a] = iirnotch(w0, bw); % 設計陷波濾波器
    y = filtfilt(b,a,x); % 零相位濾波
else % 無信號處理工具箱的回退方法
    % 粗略方法：減去窄帶正弦擬合
    t = (0:numel(x)-1)/fs; % 時間向量
    ref = sin(2*pi*f0*t); % 參考正弦波
    alpha = (ref*x.')/(ref*ref.'); % 計算最佳幅度
    y = x - alpha*ref; % 減去估計的干擾
end
end


function y = bp_filter(x, fs, f1, f2)
% 帶通濾波聚焦QRS波能量
if has_spt() % 檢查是否有信號處理工具箱
    [b,a] = butter(3, [f1 f2]/(fs/2), 'bandpass'); % 設計3階帶通巴特沃斯濾波器
    y = filtfilt(b,a,x); % 零相位濾波
else % 無信號處理工具箱的回退方法
    % 移動平均高通強調
    y = x - movmean(x, round(fs*0.15)); % 簡單高通濾波
end
end


function tf = has_spt()
% 檢測信號處理工具箱是否存在（butter/filtfilt）
tf = exist('butter','file')==2 && exist('filtfilt','file')==2;
end


function [locs, thr] = detect_r_peaks(int_sig, fs)
% 在積分后的QRS增強信號上使用自適應閾值檢測
int_sig = int_sig(:); % 確保列向量
% 使用百分位數估計噪聲和信號水平
noise_est = prctile(int_sig, 60); % 噪聲估計（第60百分位數）
sig_est   = prctile(int_sig, 98); % 信號估計（第98百分位數）
thr = 0.3*noise_est + 0.7*sig_est; % 自適應閾值


% 不應期約200毫秒
ref = round(0.2*fs); % 不應期樣本數
cand = find(int_sig > thr); % 找到超過閾值的候選點


if isempty(cand) % 如果沒有候選點
    locs = [];
    return;
end


locs = []; % 初始化R峰位置
k = 1; % 初始化索引
while k <= numel(cand) % 遍歷所有候選點
    idx = cand(k); % 當前候選點索引
    % 形成一個不應期窗口
    w_end = min(numel(int_sig), idx+ref); % 窗口結束索引
    [~, localmax] = max(int_sig(idx:w_end)); % 在窗口內找到局部最大值
    locs(end+1) = idx + localmax - 1; % 存儲局部最大值位置
    % 跳過不應期
    k = find(cand > idx+ref, 1, 'first'); % 找到下一個不應期外的候選點
    if isempty(k), break; end % 如果沒有更多候選點，退出循環
end
locs = unique(locs); % 確保位置唯一
end


function locs_ref = refine_r_peaks(ecg, locs, rad)
% 在帶通濾波后的ECG上精煉R峰位置為局部最大值
locs_ref = locs; % 初始化精煉后的位置
N = numel(ecg); % ECG信號長度
for i=1:numel(locs) % 對每個檢測到的位置
    a = max(1, locs(i)-rad); % 搜索窗口開始
    b = min(N, locs(i)+rad); % 搜索窗口結束
    [~, mx] = max(ecg(a:b)); % 在窗口內找到最大值
    locs_ref(i) = a + mx - 1; % 更新位置
end
% 移除過于接近的峰（<200毫秒）
minsep = round(0.2 * (length(ecg)/max(1,ceil(max(locs)/ (length(ecg)))))); % 最小間隔（未使用）
% 確保嚴格遞增且間隔合適
locs_ref = sort(locs_ref); % 排序
if numel(locs_ref) > 1 % 如果有多個峰
    RR = diff(locs_ref); % 計算RR間期（樣本數）
    keep = [true, RR > round(0.2* (length(ecg)/ (ecg_time_length(ecg))))]; % 保留合適間隔的峰
    locs_ref = locs_ref(keep); % 應用篩選
end
end


function T = ecg_time_length(ecg)
% 輔助函數：用于精煉間隔
T = numel(ecg); % 僅用于避免額外參數
end


function [lfnu, hfnu, lf_center, hf_center, lf_power, hf_power, total_power] = hrv_freq(NN, tNN)
% 將NN(t)插值到均勻網格并計算功率譜密度；積分LF/HF頻帶
% 頻帶（成人）：VLF 0.003-0.04 Hz（此處忽略），LF 0.04-0.15，HF 0.15-0.40
if numel(NN) < 4 % 檢查是否有足夠的數據點
    lfnu=NaN; hfnu=NaN; lf_center=NaN; hf_center=NaN;
    lf_power=NaN; hf_power=NaN; total_power=NaN; return;
end
fs_hrv = 4; % 4 Hz插值頻率
tu = tNN(1):1/fs_hrv:tNN(end); % 均勻時間網格
x  = interp1(tNN, NN, tu, 'pchip'); % 三次Hermite插值


[pxx, f] = pwelch(x - mean(x), 256, 128, 1024, fs_hrv); % Welch功率譜估計


lf_idx = f>=0.04 & f<0.15; % LF頻帶索引
hf_idx = f>=0.15 & f<=0.40; % HF頻帶索引


lf_power = trapz(f(lf_idx), pxx(lf_idx)); % LF絕對功率
hf_power = trapz(f(hf_idx), pxx(hf_idx)); % HF絕對功率
total_power = trapz(f(f>=0.04 & f<=0.40), pxx(f>=0.04 & f<=0.40)); % 總功率


lfnu = 100 * lf_power / (lf_power + hf_power); % 歸一化LF功率
hfnu = 100 * hf_power / (lf_power + hf_power); % 歸一化HF功率


lf_center = centroid(f(lf_idx), pxx(lf_idx)); % LF中心頻率
hf_center = centroid(f(hf_idx), pxx(hf_idx)); % HF中心頻率
end


function [pxx,f] = hrv_psd(NN, tNN)
% 輔助函數：顯示功率譜密度
if numel(NN) < 4 % 檢查是否有足夠的數據點
    pxx = zeros(1,512); f = linspace(0,0.5,512); return;
end
fs_hrv = 4; % 4 Hz插值頻率
tu = tNN(1):1/fs_hrv:tNN(end); % 均勻時間網格
x  = interp1(tNN, NN, tu, 'pchip'); % 三次Hermite插值
[pxx,f] = pwelch(x - mean(x), 256, 128, 1024, fs_hrv); % Welch功率譜估計
end


function c = centroid(f, p)
% 計算頻譜中心頻率
if isempty(f) || isempty(p) || sum(p)<=0 % 檢查輸入有效性
    c = NaN; return;
end
c = sum(f(:).*p(:))/sum(p(:)); % 加權平均頻率
end


function [SD1, SD2] = poincare_SD1_SD2(NN)
% SD1：短期變異性，SD2：長期變異性
NN1 = NN(1:end-1); % 前一個NN間期
NN2 = NN(2:end); % 后一個NN間期
diffs = (NN2 - NN1)/sqrt(2); % 差值分量
sums  = (NN2 + NN1)/sqrt(2); % 和分量
SD1 = std(diffs); % SD1：差值分量的標準差
SD2 = std(sums); % SD2：和分量的標準差
end


function ac1 = autocorr_stat(x, lag)
% 計算滯后自相關系數
x = x(:) - mean(x); % 去均值
if numel(x) < lag+1 % 檢查是否有足夠的數據點
    ac1 = NaN; return;
end
ac1 = sum( x(1:end-lag).*x(1+lag:end) ) / sum( x.^2 ); % 自相關計算
end


function out = ternary(cond, a, b)
% 三元操作符函數
if cond, out=a; else, out=b; end
end