关于ECG信号处理的MATLAB代码的逐行解读

Jean风居住的街道

已于 2024-03-25 16:35:03 修改

阅读量1.3k

点赞数 21

文章标签：信号处理 matlab 人工智能

于 2024-01-17 17:01:57 首次发布

本文链接：https://blog.csdn.net/qq_40781144/article/details/135650283

版权

作为即将要入门ECG信号处理的小小白，即使拿到大佬们的代码并且跑通也还是不懂每一句代码的含义及其背后的原理，这里借助chatgpt的解释把它备份到这里，一来是作为自己的笔记本用时方便查阅，二来帮助像我一样懵懂的萌新省去查阅繁杂资料的麻烦，本文旨在做笔记所以未刻意排版做美观处理，如有错误，欢迎指正
参考文章见：基于MATLAB实现ECG心电信号处理
ps：每句前面序号代表代码所在行数

fp 是通带截止频率（通带结束点），表示在滤波器的通带中保留信号的频率范围
fs 是阻带截止频率（阻带开始点），表示在滤波器的阻带中抑制信号的频率范围
rp 是通带衰减，表示在通带中滤波器允许的最大衰减
rs 是阻带衰减，表示在阻带中滤波器要达到的最小衰减
将角频率转换为弧度/秒。角频率（wp和ws）是频率的一个标准表示，以弧度为单位。这里使用了 2 * pi 的乘法因子
buttord 函数用于计算巴特沃斯滤波器的最佳阶数（n）和截止频率（wn）。它的输入参数是通带角频率、阻带角频率、通带衰减和阻带衰减，以及滤波器类型（这里是 ‘s’ 表示模拟滤波器）。n 是计算得到的最佳滤波器阶数。wn 是对应于通带截止频率的模拟滤波器的截止频率

在MATLAB中，'butter’函数的设计参数中有一个参数用于指定滤波器的类型，即’s’表示模拟滤波器，而’d’表示数字滤波器。

's' 表示使用模拟（analog）滤波器设计方法。模拟滤波器设计通常是在频率域中对连续时间系统进行设计，不考虑离散化和采样。设计完成后可以使用freqs等函数来查看模拟滤波器的频率响应。
'd' 表示使用数字（digital）滤波器设计方法。数字滤波器设计考虑了离散化和采样的影响，通常更适用于数字信号处理应用。设计完成后可以使用freqz等函数来查看数字滤波器的频率响应。

在上面的代码中，'s' 用于模拟滤波器的设计，因此 butter 函数设计的是一个连续时间（模拟）滤波器。在设计完成后，如果你希望在数字系统中应用这个滤波器，你可能需要将其离散化（数字化）。这可以通过使用c2d函数进行。例如，在这个例子中，c2d函数将模拟滤波器的系数转换为数字滤波器的系数，其中fsamp是离散化时的采样频率。

[b, a] = butter(n, wn, 'low', 's');  % 模拟滤波器设计
[b_digital, a_digital] = c2d(b, a, fsamp);  % 离散化为数字滤波器

在MATLAB的Signal Processing Toolbox中，buttap 函数用于设计巴特沃斯模拟滤波器的极点-零点-增益表示法参数。巴特沃斯滤波器是一种特殊类型的滤波器，通常被设计为低通滤波器。

[z,P,k]=buttap(n);

在 buttap(n) 中，n 是滤波器的阶数，它决定了滤波器的复杂性和截止频率的控制。阶数越高，滤波器越复杂，截止频率的控制越尖锐。

巴特沃斯滤波器的特点之一是其极点分布。在 buttap 函数返回的极点-零点-增益表示法中，z 是滤波器的极点，P 是滤波器的零点，k 是增益。对于低通巴特沃斯滤波器，极点主要分布在左半平面，而零点主要分布在原点。这种特殊的极点分布使得低通巴特沃斯滤波器在通带内有一个相对平坦的频率响应。

因此，通过 buttap(n) 函数设计的极点-零点-增益表示法通常对应于低通滤波器。这并不是唯一的设计低通滤波器的方法，但在巴特沃斯滤波器中，它是一种常见的实现方式。

在这行代码中，zp2tf 函数用于将极点-零点-增益（ZPK）表示法转换为传递函数（TF）表示法。在这个特定的情境下，z 是极点，P 是零点，k 是增益。

[bp, ap] = zp2tf(z, P, k);

z 是模拟滤波器的极点。
P 是模拟滤波器的零点。
k 是增益。

bp 是传递函数的分子系数，ap 是传递函数的分母系数。

这里的 bp 和 ap 表示的是传递函数：

在这里插入图片描述

其中，( n ) 是滤波器的阶数。

在这个上下文中，z, P, 和 k 是通过 buttap 函数得到的巴特沃斯滤波器的极点-零点-增益表示法的参数。通过 zp2tf 将这些参数转换为传递函数形式，从而获得模拟滤波器的传递函数。

在这一行代码中，lp2lp 函数被用于将一个低通滤波器（Low Pass，LP）的传递函数转换为另一个低通滤波器。这是通过指定新的截止频率 wp 来实现的。

[bs, as] = lp2lp(bp, ap, wp);

bp 和 ap 是原始低通滤波器的传递函数分子和分母系数。
wp 是新的低通滤波器的截止频率。

bs 和 as 是转换后的新低通滤波器的传递函数分子和分母系数。

在这里，通过 lp2lp 函数，原始低通滤波器的传递函数 (H_a(s)) 被转换为新低通滤波器的传递函数 (H_s(s))。这种转换通常是通过调整原始滤波器的频率响应来实现的，确保新的滤波器具有所需的截止频率 wp。这样做的目的是确保滤波器在新的截止频率下仍然具有适当的性能，以滤除肌电信号。例如，根据具体的应用需求，调整截止频率以适应不同的信号特性。

在这一行代码中，freqs 函数用于计算数字滤波器的频率响应。这个频率响应包括复数值，其中 hs 表示频率响应，ws 表示对应的频率。

[hs, ws] = freqs(bs, as);

bs 和 as 是数字滤波器的传递函数的分子和分母系数。
hs 是数字滤波器在频率域上的响应。
ws 是对应于 hs 中每个频率点的频率。

hs 和 ws 的关系表示了数字滤波器在不同频率下的增益或衰减。

此步骤通常用于了解数字滤波器在频域中的性能，以便进一步分析其对输入信号的影响。可以通过绘制频率响应图来更直观地理解数字滤波器对信号的影响。

在这一行代码中，bilinear 函数用于对数字滤波器进行双线性变换。这个过程是将连续时间（模拟）滤波器的传递函数转换为离散时间（数字）滤波器的传递函数。

[bz, az] = bilinear(bs, as, Fs);

bs 和 as 是连续时间数字滤波器的传递函数的分子和分母系数。
Fs 是数字滤波器的采样频率。
bz 和 az 是通过双线性变换得到的数字滤波器的传递函数的分子和分母系数。

通过这个步骤，模拟滤波器被离散化以适应数字信号的处理。这是因为在实际应用中，信号通常以数字形式进行采样和处理。使用双线性变换可以将连续时间滤波器转换为数字滤波器，以便对离散时间信号进行滤波。

在这段代码中，进行了离散傅立叶变换（FFT）以获取信号的频谱信息。以下是对代码的解释：

N = 512;                      % FFT的点数
n = 0:N-1;                    % 时间序列
mf = fft(M(:, 1), N);         % 进行频谱变换（傅里叶变换）
mag = abs(mf);                % 计算频谱的幅度
f = (0:length(mf)-1) * Fs / length(mf);  % 计算频率轴

N = 512 设置了FFT的点数，这决定了频谱分析的精度。
n = 0:N-1 创建了时间序列，用于表示信号的时间轴。
fft(M(:, 1), N) 对信号 M(:, 1) 进行离散傅立叶变换，得到频谱信息。结果存储在 mf 中。
abs(mf) 计算了频谱的幅度，即信号在频域中的强度。
f = (0:length(mf)-1) * Fs / length(mf) 计算了频率轴，其中 Fs 是采样频率。这个轴表示FFT输出中每个频率点对应的实际频率。

通过这段代码，你可以得到信号的频谱信息，并可以使用 mag 和 f 进行进一步的频谱分析和可视化。

在这段代码中，进行了心电信号的功率谱密度估计，并通过图形表示。以下是对代码的解释：

wn = M(:, 1);                 % 从信号中提取一列作为处理的信号
P = 10 * log10(abs(fft(wn).^2)/N);  % 计算信号的功率谱密度估计
f = (0:length(P)-1) / length(P);   % 计算频率轴

figure
plot(f, P);                   % 绘制功率谱密度曲线
grid
xlabel('归一化频率');
ylabel('功率(dB)');
title('心电信号的功率谱');

wn = M(:, 1) 从原始信号中提取一列，作为需要处理的信号。
fft(wn) 对提取的信号进行傅里叶变换，得到频域表示。
abs(fft(wn).^2)/N 计算频域表示的幅度的平方，并除以FFT的点数，得到功率谱密度估计。
10 * log10(...) 将功率转换为分贝（dB）单位，以便更容易观察。
f = (0:length(P)-1) / length(P) 计算频率轴。
plot(f, P) 绘制功率谱密度曲线。
grid 添加网格线以提高可视化效果。
xlabel('归一化频率') 和 ylabel('功率(dB)') 添加轴标签。
title('心电信号的功率谱') 添加图形标题。

通过这段代码，你可以可视化心电信号的功率谱密度估计，以更好地了解信号在频域上的特性。

这行代码计算了信号 wn 的功率谱密度估计：

P = 10 * log10(abs(fft(wn).^2)/N);

fft(wn)：对信号 wn 进行离散傅里叶变换，得到频域表示。这里使用 fft 函数。
abs(fft(wn).^2)：计算频域表示的幅度的平方。这是为了获取功率谱的估计，因为功率是幅度的平方。
/N：将结果除以FFT的点数，这是为了得到归一化的功率谱密度。
10 * log10(...)：将结果取对数并乘以10，以转换为分贝（dB）单位。这样做是为了更方便地表示信号的强度。

综合起来，这行代码计算了信号在频域中每个频率分量的功率谱密度估计，并将结果以分贝为单位表示。功率谱密度估计用于了解信号在频域中不同频率上的功率分布，提供了一种分析信号频谱特性的方法。

这段代码设计了一个带陷滤波器，其目的是抑制特定频率，即50 Hz的工频。以下是对代码的逐行解释：

Me = 100;               % 滤波器阶数
L = 100;                % 窗口长度
beta = 100;             % 衰减系数
Fs = 1500;              % 采样频率
wc1 = 49/Fs*pi;         % 高通滤波器截止频率，对应51Hz
wc2 = 51/Fs*pi;         % 低通滤波器截止频率，对应49Hz
h = ideal_lp(0.132*pi, Me) - ideal_lp(wc1, Me) + ideal_lp(wc2, Me); % 陷波器冲激响应
w = kaiser(L, beta);
y = h .* rot90(w);      % 50Hz陷波器冲激响应序列
m2 = filter(y, 1, m);

Me 是滤波器的阶数，决定了滤波器的复杂度。
L 是窗口长度，用于设计Kaiser窗口。
beta 是Kaiser窗口的衰减系数，控制窗口的形状。
Fs 是采样频率。
wc1 和 wc2 是高通滤波器和低通滤波器的截止频率，分别对应于51 Hz和49 Hz。这两个截止频率一起定义了带陷滤波器的通带范围，排除了50 Hz的工频。wc1 和 wc2 并不直接对应于高通或低通滤波器的截止频率，而是用于引入零点的频率。这两个零点的存在是为了在设计的带陷滤波器中形成两个抑制带，使其在50 Hz的电力线频率附近产生一个“带阻”效果。
ideal_lp 函数是一个理想低通滤波器的生成函数。在这里，它被用于生成一个理想的带陷滤波器的冲激响应。
kaiser 函数生成一个Kaiser窗口，用于加权滤波器的冲激响应。
rot90(w) 将Kaiser窗口旋转，以使其与滤波器的冲激响应对齐。
h .* rot90(w) 是带陷滤波器的冲激响应，通过Kaiser窗口进行加权。
filter(y, 1, m) 使用带陷滤波器对输入信号 m 进行滤波。这将抑制50 Hz的工频成分，而保留其他频率的信号成分。