心电信号基于matlab瞬时抑制心电信号IIR滤波含Matlab源码 1533期

Posted 2021-12-24 紫极神光

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【心电信号】基于matlab心电信号PTT+HRV+PRV【含Matlab源码 1551期】

获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、心电信号简介

0 引言
心电信号是人类最早研究的生物信号之一, 相比其他生物信号更易于检测, 且具有直观的规律。心电图的准确分析对心脏病的及早治疗有重大的意义。人体是一个复杂精密的系统, 有许多不可抗的外界因素, 得到纯净的心电信号非常困难。可以采用神经网络算法去除心电信号的噪声, 但这种方法存在训练难度大、耗时长的缺点。小波变换在处理非线性、非平稳且奇异点较多的信号时具有一定的优越性, 近年来许多学者使用其对心电信号进行研究。

1 心电信号简介
心电信号由以下几个波段组成, 一个典型的心电图如图1所示。

图1 典型心电图
(1) P波:反映心房肌在除极过程中的电位变化过程;
(2) P-R间期:反映的是激动从窦房结通过房室交界区到心室肌开始除极的时限;
(3) QRS波群:反映心室肌除极过程的电位变化;
(4) T波:代表心室肌复极过程中所引起的电位变化;
(5) S-T段:从QRS波群终点到达T波起点间的一段水平线[2];
(6) Q-T间期:心室从除极到复极的时间[3];
(7) U波:代表动作电位的后电位。
由于心电信号十分微弱, 且低频, 极易受到干扰, 不同的干扰源的噪声虽是随机的, 但来自同一个干扰源的噪声往往具有同一类特征。分析干扰的来源, 针对不同的来源使用合适的处理方法, 是数据采集重点考虑的一个问题。常见干扰有3种: (1) 工频干扰; (2) 基线漂移; (3) 肌电干扰。其中已经证明小波变换在抑制心电信号的工频干扰方面具有较大优势。具体噪声频带如表1所示。
表1 心电信号以及主要噪声频带

三、部分源代码

%% Preliminaries.
clc;
clear all;
close all;


%% Define execution parameters.

% Noise (AC interference) parameters.
A0 = 250;   % (Peak) amplitude [uV].
fd = 50;    % Frequency [Hz].
phi = 15;   % Phase [deg].


% Filtering parameters.
BW = 0.8;   % Notch bandwidth [Hz].
M = 10;     % Number of initial samples to consider for the transient
            % suppression technique.


%% Load data. 

% Data is contained in "ecg.mat".
% The original data was obtained from the MIT-BIH Polysomnographic
% Database (http://www.physionet.org/cgi-bin/atm/ATM).
% ecgOriginal   Original (clean) ECG signal (given by Eq. 3).
% fs            Sampling frequency [Hz].
load ecg;

s_n = ecgOriginal;              % To follow paper's notation.

% Original data consists of 10 s recording.
% For simplicity, we will consider only the first 5 s.
s_n = s_n(1:end/2);


%% Calculate important parameters.

N = numel(s_n);                 % Number of samples.
Ts = 1/fs;                      % Sampling period [s].
t = linspace(0, (N-1)*Ts, N);   % Time vector [s].

w0 = 2 * pi * (fd/fs);          % Notch frequency [rad/s].
omega = 2 * pi * (BW/fs);       % Bandwidth [rad/s].


%% Contaminate original signal.

% Create noise, i.e. AC interference (Eq. 4).
d_n = A0 * sin((2 * pi * fd) .* t + deg2rad(phi));

% Add the original signal with the noise (Eq, 3).
x_n = s_n + d_n;


%% Filtering process 1 (conventional).

% 1. Calculate a1 and a2 coefficients using Eq. 2.
a1 = (2 * cos(w0)) / (1 + tan(omega/2));
a2 = (1 - tan(omega/2)) / (1 + tan(omega/2));

% 2. Choose arbitrary initial conditions (x[-1], x[-2], y[-1], y[-2]).
x_1 = 0;    % x[-1].
x_2 = 0;	% x[-2].
y_1 = 0;	% y[-1].
y_2 = 0;	% y[-2].

% 3. From n = 0 to N,  calculate the output, given by Eq. 5.

% Manually calculate the first two samples (for sake of clarity).
y_n(1) = 0.5 * ((1 + a2)* x_n(1) - 2*a1*x_1    + (1 + a2)*x_2) + (a1*y_1)    - (a2*(y_2));



%% Filtering process 2 (notch filtering with transient state suppression).
y_n = ecgTransientSuppression(x_n, fs, fd, BW, M);
ecgFilt2 = y_n;


%% Plots.

figure('Name', 'Input Signal');
subplot(3,1,1);
plot(t,s_n,'b');
title('Original (Clean) ECG');
ylabel('Amplitude [\\muV]')
subplot(3,1,2);
plot(t,d_n,'b');
title('Noise (AC Interference)');
ylabel('Amplitude [\\muV]')
subplot(3,1,3);
plot(t,x_n,'b');
title('ECG + Noise');
xlabel('Time [s]');
ylabel('Amplitude [\\muV]')

figure('Name','Filtering Comparison');
subplot(3,1,1);
plot(t,s_n,'b');
title('Original (Clean) ECG');
ylabel('Amplitude [\\muV]')
subplot(3,1,2);
plot(t,ecgFilt1,'b');
title('Filtered ECG (Conventional Method)');
ylabel('Amplitude [\\muV]')
subplot(3,1,3);
plot(t,ecgFilt2,'b');
title('Filtered ECG (Transient Suppression Method)');
xlabel('Time [s]');
ylabel('Amplitude [\\muV]')

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.
[4]焦运良，邢计元，靳尧凯.基于小波变换的心电信号阈值去噪算法研究[J].信息技术与网络安全. 2019,38(05)