MATLAB读取MIT-BIH数据库的ECG数据利用差分和小波算法进行PQRST波定位

从MIT-BIH数据库读取ECG数据

首先需要正确读取MIT-BIH的数据文件（通常包括 .dat, .hea, .atr 文件）。

1. 数据文件简介：

   • .dat 文件：存储ECG信号数据。
   • .hea 文件：头文件，包含通道数、采样频率、信号增益等信息。
   • .atr 文件：注释文件，包含由专家标记的心搏类型等信息。

2. 读取212格式数据：
   MIT-BIH数据库的部分数据采用一种称为"212"的格式存储，每三个字节存储两个12位的采样值。以下是一个参考读取函数，此函数基于网络上流传较广且被多次引用的 rddata.m 脚本思路：

   function [ecg1, ecg2, time] = read_mit_bih_212(filename, samples2read)
   % READ_MIT_BIH_212 读取MIT-BIH 212格式的ECG数据
   % 输入：
   %   filename - 数据文件名（不包含扩展名）
   %   samples2read - 需要读取的样本数（每个通道）
   % 输出：
   %   ecg1, ecg2 - 两个通道的ECG信号
   %   time - 时间向量（秒）
   
   % 读取头文件信息
   headerfile = [filename, '.hea'];
   fid = fopen(headerfile, 'r');
   if fid == -1
       error('无法打开头文件：%s', headerfile);
   end
   
   % 解析头文件获取基本信息，例如采样频率和通道数
   z = fgetl(fid);
   A = sscanf(z, '%*s %d %d %d', [1, 3]);
   num_channels = A(1);       % 信号通道数
   sfreq = A(2);              % 数据采样频率
   total_samples = A(3);      % 每个通道的信号样本数
   fclose(fid);
   
   % 根据头文件信息，确定信号格式（例如212格式）和增益等
   % 此处可能需要根据具体头文件内容进行更详细的解析
   
   % 读取.dat文件
   datafile = [filename, '.dat'];
   fid = fopen(datafile, 'rb');
   if fid == -1
       error('无法打开数据文件：%s', datafile);
   end
   
   % 读取二进制数据，212格式每三个字节存储两个12位的采样值
   raw_data = fread(fid, [3, samples2read], 'uint8')';
   fclose(fid);
   
   % 将212格式的数据转换为物理信号值
   % 提取每个信号的高4位和低8位
   M2H = bitshift(raw_data(:, 2), -4);        % 第二个信号的高4位
   M1H = bitand(raw_data(:, 2), 15);          % 第一个信号的高4位
   
   % 处理符号位（12位数据，最高位为符号位）
   PRL = bitshift(bitand(raw_data(:, 2), 8), 9);   % 第一个信号的符号位
   PRR = bitshift(bitand(raw_data(:, 2), 128), 5); % 第二个信号的符号位
   
   % 组合形成两个通道的原始整数信号
   ecg1_raw = bitshift(M1H, 8) + raw_data(:, 1) - PRL;
   ecg2_raw = bitshift(M2H, 8) + raw_data(:, 3) - PRR;
   
   % 此处应根据头文件中的增益和零点信息，将原始整数信号转换为物理值（例如mV）
   % gain1 = ... ; zerovalue1 = ... ; 需要从.hea文件解析
   % ecg1 = (ecg1_raw - zerovalue1) / gain1;
   % ecg2 = (ecg2_raw - zerovalue2) / gain2;
   % 为演示，我们暂时直接使用原始整数信号，并假设增益和零点
   ecg1 = ecg1_raw;
   ecg2 = ecg2_raw;
   
   % 创建时间向量
   dt = 1 / sfreq;
   time = (0 : dt : (samples2read-1)*dt)';
   end
   

   注意：上述代码中的增益（gain）和零点（zerovalue）需要您根据具体使用的MIT-BIH记录的头文件（.hea）进行解析和赋值。直接使用原始整数信号 ecg1_raw 和 ecg2_raw 进行后续处理在算法演示上是可行的，但若需物理量纲（如mV），则必须进行转换。

PQRST波定位算法

ECG信号通常包含噪声，在特征波定位之前，预处理（去噪）非常重要。

1. 预处理（滤波）
   一个常见的做法是使用带通滤波器，例如通带频率1.2-49Hz，可以去除基线漂移和高频噪声。

   % 假设信号已读取到变量 ecg1 中，采样频率为 Fs
   orden = 200; % 滤波器阶数
   Fp1 = 1.2; Fp2 = 49; % 通带频率（Hz）
   Fp1_n = Fp1/(Fs/2); Fp2_n = Fp2/(Fs/2); % 归一化频率
   b = fir1(orden, [Fp1_n, Fp2_n], 'bandpass'); % FIR滤波器设计
   filtered_ecg = filtfilt(b, 1, ecg1); % 零相位滤波
   

2. 基于差分算法的QRS波检测
   差分算法计算量小，适合实时检测。

   function [r_peaks, qrs_onsets, qrs_offsets] = detect_qrs_diff(ecg, Fs)
   % 基于差分算法检测QRS复合波
   % 1. 微分（近似一阶差分）
   diff_ecg = diff(ecg);
   % 2. 平方（增强高频成分，突出QRS波）
   squared_ecg = diff_ecg .^ 2;
   % 3. 移动平均积分（平滑信号）
   window_width = round(0.08 * Fs); % 80ms窗宽
   ma_ecg = movmean(squared_ecg, window_width);
   
   % 4. 自适应阈值检测R峰
   % 寻找ma_ecg中的峰值，使用findpeaks
   [peaks, locs] = findpeaks(ma_ecg, 'MinPeakHeight', mean(ma_ecg)*2, ... 
                              'MinPeakDistance', round(0.6 * Fs)); % 最小峰间隔0.6秒
   
   r_peaks = locs + 1; % 补偿diff操作导致的位置偏移
   
   % 5. 搜索Q点和S点（在R峰前后的小范围内寻找原始ECG信号的最低点）
   q_locs = zeros(size(r_peaks));
   s_locs = zeros(size(r_peaks));
   search_win = round(0.05 * Fs); % R峰前后50ms搜索窗
   
   for i = 1:length(r_peaks)
       % 搜索Q点（R峰之前）
       start_idx = max(1, r_peaks(i) - search_win);
       [~, min_idx] = min(ecg(start_idx:r_peaks(i)));
       q_locs(i) = start_idx + min_idx - 1;
       
       % 搜索S点（R峰之后）
       end_idx = min(length(ecg), r_peaks(i) + search_win);
       [~, min_idx] = min(ecg(r_peaks(i):end_idx));
       s_locs(i) = r_peaks(i) + min_idx - 1;
   end
   
   % QRS起始和结束通常可近似为Q点之前和S点之后的一个点，或通过更复杂的逻辑判断
   qrs_onsets = q_locs; % 简化处理
   qrs_offsets = s_locs; % 简化处理
   end
   

   提示：差分算法对QRS波群的提取准确率较高，但参数（如阈值、窗宽）可能需要根据具体信号调整。

3. 基于小波算法的P波、T波检测
   小波变换擅长分析非平稳信号，能有效检测P波和T波。

   function [p_peaks, t_peaks] = detect_p_t_wave(ecg, Fs, r_peaks, qrs_onsets, qrs_offsets)
   % 检测P波和T波
   % 1. 使用平稳小波变换（SWT）或提升小波变换
   % 此处以SWT为例，使用'sym4'小波，分解到第5层
   num_levels = 5;
   wname = 'sym4';
   [swa, swd] = swt(ecg, num_levels, wname); % 平稳小波变换
   
   % 2. P波检测：通常在QRS起始点之前，关注低频近似分量
   % 使用第5层近似分量(swa(5,:))，它包含了P波的主要信息
   p_peaks = [];
   for i = 1:length(qrs_onsets)
       % 在QRS起始点前150ms到50ms范围内搜索P波
       search_start = max(1, qrs_onsets(i) - round(0.15 * Fs));
       search_end = max(1, qrs_onsets(i) - round(0.05 * Fs));
       if search_end > search_start
           segment = swa(5, search_start:search_end);
           [p_max, p_idx] = max(segment);
           % 设置阈值判断
           if p_max > 0.5 * std(swa(5, :)) % 阈值需调整
               p_peaks(end+1) = search_start + p_idx - 1;
           end
       end
   end
   
   % 3. T波检测：通常在QRS结束点之后，也关注低频近似分量
   t_peaks = [];
   for i = 1:length(qrs_offsets)
       % 在QRS结束点后50ms到400ms范围内搜索T波
       search_start = min(length(ecg), qrs_offsets(i) + round(0.05 * Fs));
       search_end = min(length(ecg), qrs_offsets(i) + round(0.4 * Fs));
       if search_end > search_start
           segment = swa(5, search_start:search_end);
           [t_max, t_idx] = max(segment);
           % 设置阈值判断
           if t_max > 0.5 * std(swa(5, :)) % 阈值需调整
               t_peaks(end+1) = search_start + t_idx - 1;
           end
       end
   end
   end
   

   注意：小波变换检测P波和T波有一定难度，因为它们的形态和幅度多变，且容易受噪声影响。上述阈值 (0.5 * std(swa(5, :))) 是示例，需要根据实际信号仔细调整。使用平稳小波变换（SWT） 可以降低平移敏感性，提高检测精度。

参考代码 MATLAB读取MIT库的ECG数据程序 www.3dddown.com/cna/64125.html

主程序流程示例

将上述模块整合起来：

clear all; close all; clc;

% 1. 读取数据
filename = '117'; % 例如MIT-BIH记录117
samples2read = 360 * 10; % 例如读取10秒的数据，假设采样频率360Hz
[ecg1, ecg2, time] = read_mit_bih_212(filename, samples2read);
Fs = 360; % 根据头文件信息设置采样频率

% 2. 预处理（滤波）
orden = 200;
Fp1 = 1.2; Fp2 = 49;
Fp1_n = Fp1/(Fs/2); Fp2_n = Fp2/(Fs/2);
b = fir1(orden, [Fp1_n, Fp2_n], 'bandpass');
filtered_ecg = filtfilt(b, 1, ecg1);

% 3. QRS波检测 (差分算法)
[r_peaks, qrs_onsets, qrs_offsets] = detect_qrs_diff(filtered_ecg, Fs);

% 4. P波和T波检测 (小波算法)
[p_peaks, t_peaks] = detect_p_t_wave(filtered_ecg, Fs, r_peaks, qrs_onsets, qrs_offsets);

% 5. 可视化结果
figure;
plot(time, filtered_ecg); hold on;
plot(time(r_peaks), filtered_ecg(r_peaks), 'rv', 'MarkerFaceColor', 'r', 'DisplayName', 'R峰');
plot(time(p_peaks), filtered_ecg(p_peaks), 'g^', 'MarkerFaceColor', 'g', 'DisplayName', 'P峰');
plot(time(t_peaks), filtered_ecg(t_peaks), 'bs', 'MarkerFaceColor', 'b', 'DisplayName', 'T峰');
legend show;
xlabel('时间 (秒)');
ylabel('幅度 (mV)');
title('ECG信号PQRST波定位');
grid on;

叮一下

• 数据路径：确保MATLAB的当前工作目录或搜索路径包含您的MIT-BIH数据文件。
• 参数调整：示例中的算法参数（如阈值、窗宽）可能需要根据您读取的具体MIT-BIH记录和预处理效果进行调整。
• 算法优化：提供的算法是基础实现。在实际研究中，可能需要结合更多策略（如使用ECG的斜率信息、模板匹配、机器学习等）来提高检测的鲁棒性和准确性，尤其是在处理有噪声或心律不齐的信号时。
• 性能评估：可以使用MIT-BIH数据库附带的专家标记注释文件（.atr）来定量评估您算法的检测准确率。