ECG 訊號偵測 QRS 流程

流程

1. 資料預處理
2. Pan & Tompkin 演算法找 R 波
3. 重新定位 R 波位置

訊號來源說明

這裡使用的訊號是台大的腦磁圖儀 ( MEG )

500多個通道中的 ECG 通道

而其中 ECG003、ECG004 通道是有資料的

資料預處理

這裡使用 FieldTrip 來處理 FIF 檔

取出 ECG003、ECG004 通道

並透過訊號事件來去除前後的雜訊

%% 設定資料
fname = 'exp2shortrun1_sss.fif';
%% 清除設定
cfg = [];
%% 設定選擇的資料
cfg.dataset = fname;
%% 設定訊號通道，這裡用'ECG003','ECG004'
cfg.channel    = {'ECG003','ECG004'}; %%設定訊號通道，這裡用ECG
%% 資料預處理
data_orig = ft_preprocessing(cfg);
%% 偵測訊號事件位置
event = ft_read_event(cfg.dataset);
%% 找到起始事件位置
sample_onset = [event(find(strcmp('STI001', {event.type}))).sample]';
%% 找到最後事件位置
sample_offset = [event(find(strcmp('STI002', {event.type}))).sample]';
%% 清除設定
cfg = [];
%% 選擇訊號長度(這裡單位為秒)
%% 原始資料長度為 1 ~ 349 秒
%% 底下擷取長度為 24 ~ 340 秒
cfg.toilim = [(sample_onset+20000)/data_orig.fsample sample_offset/data_orig.fsample]; 
%% 根據上方設定的訊號長度，重新調整訊號長度
data_new = ft_redefinetrial(cfg, data_orig);

選擇欲偵測的訊號

這裡使用 ECG003 通道進行示範

%% 選擇通道
Choose_ECG_Channel = 1;
%% 重組陣列，並劃出原始訊號
if Choose_ECG_Channel == 1
    ECG003 = [];
    ECG003(:,1)=data_new.trial{1,1}(1,:);
    figure(100);
    plot(ECG003); %% 畫出 ECG003
elseif Choose_ECG_Channel == 2
    ECG004 = []; 
    ECG004(:,1)=data_new.trial{1,1}(2,:);
    figure(101);
    plot(ECG004); %% 畫出 ECG004
end

檢視原始訊號

放大檢視訊號

Pan & Tompkin 演算法找 R 波

Pan & Tompkin Matlab Code 來源

https://www.mathworks.com/matlabcentral/fileexchange/45840-complete-pan-tompkins-implementation-ecg-qrs-detector

這裡有三個重要的參數

1. Fs： 取樣率
2. qrs_amp_raw：QRS 波的振幅

3. qrs_i_raw：QRS 波的索引位置

   取樣率務必設定成跟預處理過後的取樣率一樣

%% pan_tompkin 尋找R波
Fs = 1000; %% 設定取樣率

if Choose_ECG_Channel == 1
    figure(102);
    [qrs_amp_raw,qrs_i_raw,delay] = pan_tompkin(ECG003,Fs,1);
elseif Choose_ECG_Channel == 2
    figure(103);
    [qrs_amp_raw,qrs_i_raw,delay] = pan_tompkin(ECG004,Fs,1);
end

Pan & Tompkin 的執行流程

流程圖

這裡可以看到 QRS 波的位置都被找出來了

放大檢視 QRS 波

放大後可以發現濾波後的訊號，找到的 R 波位置

對應原始的訊號，有往左偏的現象 ( 延遲 )

所以後面會在寫個小程式來修復這個問題

重新定位 R 波位置

這裡我是寫成副程式

以方便重複使用

核心理念：為了解決濾波後的 R 波位置有延遲現象，因此透過一個小範圍來尋找波峰值

透過 max 函數尋找最大值

function [r_new_index,r_new_value]=ReLocate_R_Point(data,Search_Range,R_Index)

r_new_value = [];
r_new_index = [];

for i = 1:length(R_Index)

    index_range = R_Index(i)-Search_Range:R_Index(i)+Search_Range;
    value = data(index_range);
    [row,col] = find(value == max(value),1);
    r_new_value(i) = value(row);
    r_new_index(i) = index_range(row);

end
end

使用語法

如果重新找到的位置，還是有所偏移的話，可以把搜尋範圍調大一點

[r_new_index,r_new_value] = ReLocate_R_Point(ECG003,15, qrs_i_raw);
  %索引值       %振幅                         %訊號 %搜尋範圍  %Pan_Tompkin找到的索引值 

主程式範例

%% 重新定位R波
if Choose_ECG_Channel == 1
    [r_new_index,r_new_value] = ReLocate_R_Point(ECG003,15,qrs_i_raw);
    plot(ECG003,'-o','MarkerIndices',r_new_index,...
    'MarkerFaceColor','red',...
    'MarkerSize',4);
elseif Choose_ECG_Channel == 2
    [r_new_index,r_new_value] = ReLocate_R_Point(ECG004,15,qrs_i_raw);
    plot(ECG004,'-o','MarkerIndices',r_new_index,...
    'MarkerFaceColor','red',...
    'MarkerSize',4);
end

重新找 R 波的結果

放大檢視 重新找 R 波的結果

到這裡就差不多了，在接下來就看你要做甚麼分析了

像 HRV ( 心跳變異率 ) 、心率之類的