心电信号处理简介
1 引 言
ECG是一种基本的人体生理信号, 具有重要的临床诊断价值。其特点是信号微弱, 信噪比小, 一般正常人的心电信号频率在0.05～100 Hz范围内, 幅度为10 μV (胎儿) ～5 mV (成人) 。在检测心电信号时, 易受到仪器、人体活动等方面的影响, 所检测的心电信号常伴有干扰。心电信号的干扰主要有以下3种:工频干扰、基线漂移和肌电干扰。心电信号的干扰对心电数据分析和处理有很大的影响。因而, 利用数字滤波技术滤除心电信号中存在的各种干扰, 准确检测心电信号, 是心电信号处理的重要工作。
自适应信号处理是近十年来发展迅速的一种信号处理方法。小波分析也是近年来发展起来的一种新的数学理论和方法，并在许多领域得到了成功的应用。作为一种新的时频分析方法，小波分析被称为“数学显微镜”，因为它具有多分辨率分析的特点，可以关注信号的任何细节，进行多分辨率时频域分析。
国内刊物对ECG信号的消噪问题研究虽多, 但一般只是针对某种噪声干扰, 或仅借助于经典滤波方法, 且因肌电干扰较难滤除而对其论述甚少。我们则综合利用LMS自适应算法及小波理论等现代数字信号处理方法, 对心电信号检测过程中引入的上述3种噪声干扰, 系统而有针对性的设计了相应滤波器——自适应陷波器、小波变换滤波器及自适应信号分离器, 来滤除50 Hz工频干扰、基线漂移、肌电干扰等噪声, 特别是对于肌电干扰的滤除进行了尝试。经仿真与实验结果表明, 上述滤波器能有效地消除相应噪声, 得到滤波效果较好的心电信号。
2 滤波器的算法实现
2.1 工频干扰与自适应陷波器
自适应滤波器的典型应用是自适应噪声对消器。噪声对消使用在信号很弱, 或者信号不可检测的噪声场中, 将从一个或多个传感器取得的辅助输入或者参考输入加以过滤, 并从包含信号和噪声的原始输入中减去。结果, 原始噪声就受到衰减, 或者由于对消而被除掉。图1为自适应噪声对消器的原理图。它有两个输入:原始输入与参考输入。原始输入为受干扰信号dj (dj=s+v0) ;参考输入是信号v1, 它与干扰信号v0相关但与信号s不相关。图中自适应滤波器AF (Adaptive filter) 接受误差信号ej的控制调整滤波器权向量Wj, 使输出yj趋于等于dj中与它相关的v0, 于是ej作为dj与yj之差就非常接近或等于信号被检测信号s[3]。
图1 自适应噪声对消器
现在来证明这一结论。设s、v0、v1都是平稳随机过程并具有零均值 (因而yj也如此) 。由于
ej=dj-yj=s+v0-yj (1)
取ej的自相关得:
如前所述, 自适应过程就是自动调节滤波器权向量Wj使得E[ej2]=min的过程。上式第一项为信号功率E[s2], 它与Wj无关。第三项等于零, 因为s与v0、v1及yj均不相关, 所以要使E[ej2]最小等价于式 (2) 第二项最小, 即:
E[e2j]min⇔E[ (v0-yi) 2]min (3)
由式 (1) 得:
v0-yj=ej-s (4)
所以当E[ (v0-yj) 2]被最小化时, E[ (ej-s) 2]也被最小化了, 即ej以最小均方误差趋于s, 可能的最好情况为yj=v0, 则ej=s。
因此, 自适应滤波器可以从噪声中提取信号。虽然上述结果用维纳滤波器也能实现, 但是设计维纳滤波器需要预先知道s与v0或v1的统计特性, 而自适应滤波器不需要, 并且当信号或噪声统计特性变化时, 自适应滤波器也能自适应地调节它的冲击响应特性来适应这一变化。
自适应噪声抵消器的一个典型应用是自适应陷波器。如果信号中的噪声是单色干扰 (频率为ω0的正弦波干扰) , 则消除这种干扰的方法是应用陷波器。 用自适应滤波器组成的陷波器与一般固定网络的陷波器比较有下列优点:
(1) 能够自适应地准确跟踪干扰信号的频率;
(2) 容易控制带宽。
图2给出了一个具有两个权的单频干扰对消器组成的陷波器。原始输入为:任意信号s (t) 与单频干扰Acos (ω0t+a) 的叠加, 经采样后送入dj端, 故有
dj=sj+Acos (ω0jT+a) (5)
图2 自适应陷波器
参考输入为一标准正弦波cos (ω0t) , 经采样后送入x1j及x2j端, 其中后者经过90°相移, 因而有x1j=cos (ω0jT) , x2j=sin (ω0jT) , 两个权值w1j及w2j使得组合后的正弦波幅值和相位都可以调整, 因为两个权表示有两个自由度调整。经组合相加后得到yj其幅度和相位都可以与原始输入中的干扰分量相同, 使输出ej中的单频干扰得以抵消, 达到陷波的目的。由于心电信号中的工频干扰为频率固定的50 Hz正弦波信号, 故可用自适应陷波器滤除。
2.2 基线漂移与小波变换滤波器
基线漂移是低频干扰信号, 频率范围为0.05～2.00 Hz。常用的FIR和IIR滤波器截止频率固定, 在噪声频率超过其截止频率时, 无法发挥滤波作用。小波分析由于具有多分辨分析的特点而被誉为“数学显微镜”。因此, 可以利用小波的这一特性进行基线漂移的滤除。
小波变换是由加窗傅立叶变换发展起来的, 具有时频局域化的特性。小波变换可以将信号分解成不同的频段, 这为处理ECG信号提供了一个有力的手段。连续小波定义如下:
其中:a≠0为尺度因子;b为时移因子;函数Ψa, b (t) 称做小波。
在实际计算过程中, 连续小波必须加以离散化, 最常用的方法是对尺度进行二进制离散, 即令a=2j, b取整数。设x (n) 为离散信号, 按正交小波基在第j层上展开如下:
x (n) =D2j[x (n) ]+A2j[x (n) ] (7)
其中:D2j为细节信号, 代表第j层上高频分量;A2j为逼近信号, 表示该层的低频分量。
在小波变换中应选择合适的母小波易于信号的分解和重构。选用Coif小波作为分析小波。
2.3 肌电干扰与自适应信号分离器
肌电干扰是由于人体运动、肌肉收缩而引起的干扰, 影响心电信号的肌电信号主要是骨骼肌产生的动作电位, 它的频带是2～2 000 Hz。同工频干扰一样, 肌电干扰也会掩盖心电波形中的细节信息, 使得心电波形难以辨认。
参考输入是原始信号的k步延时的自适应对消器可以组成自适应预测系统、谱线增强系统以及信号分离系统, 下面讨论如何用自适应信号分离器来滤除肌电干扰[3]。
图3表示一个用作信号分离目的的系统, 当输入中包括两种成分:宽带信号与周期信号时, 为了分离这两种信号, 可以一方面将该输入送入原始输入dj端, 另一方面把它延时足够长时间后送入AF的参考输入v1端。经过延时后宽带成分已与原来的输入不相关, 而周期性成分延时前后保持强相关。图中s0为周期信号, v0为宽带信号。
于是在ej输出中将周期成分抵消只存在宽带成分, 在yj输出中只存在周期成分, 此时, AF自动调节滤波器权向量Wj, 以达到对周期成分起选通作用。算法推导同2.1。
图3 自适应信号分离器
2.部分源码
function varargout = ECG(varargin)
% ECG M-file for ECG.fig
% ECG, by itself, creates a new ECG or raises the existing
% singleton*.
%
% H = ECG returns the handle to a new ECG or the handle to
% the existing singleton*.
%
% ECG('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in ECG.M with the given input arguments.
%
% ECG('Property','Value',...) creates a new ECG or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before ECG_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to ECG_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES
% Edit the above text to modify the response to help ECG
% Last Modified by GUIDE v2.5 06-Mar-2022 22:10:59
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @ECG_OpeningFcn, ...
'gui_OutputFcn', @ECG_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before ECG is made visible.
function ECG_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to ECG (see VARARGIN)
% Choose default command line output for ECG
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes ECG wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = ECG_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
varargout{1} = handles.output;
% --- Executes on button press in open_file.
function open_file_Callback(hObject, eventdata, handles)
% hObject handle to open_file (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
[filename,pathname,q]=uigetfile('*.txt','select the Data file');%打开文件
if q>0
file=fullfile(pathname,filename);
original_signal=load(file);
end
original_signal=original_signal(:,1);
[m,n]=size(original_signal);
t=zeros(m,1);%横坐标转换为时间s
for i=1:m
t(i,1)=i/500;
end
original_signal=original_signal./1000;%纵坐标转换为mV
y_top=max(original_signal)+0.1;%纵坐标范围
y_bottom=min(original_signal)-0.1;
axes(handles.axes1);
plot(t(:,1),original_signal(:,1));
axis([0,10,y_bottom,y_top]);grid on;
handles.t=t;%横坐标时间存入句柄
handles.original_signal=original_signal;%原始信号存入句柄
guidata(hObject,handles);
% --- Executes on button press in spectrum.
function spectrum_Callback(hObject, eventdata, handles)
% hObject handle to spectrum (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
denoised_signal=handles.denoised_signal;
noised_signal=handles.noised_signal;
original_signal=handles.original_signal;
mm=512;
spect_A=plot_spectrum(original_signal,mm);
spect_B=plot_spectrum(noised_signal,mm);
spect_C=plot_spectrum(denoised_signal,mm);
axes(handles.axes4);
plot(spect_A,'r','LineWidth',2);grid on;
axes(handles.axes5);
plot(spect_B,'r','LineWidth',2);grid on;
axes(handles.axes6);
plot(spect_C,'r','LineWidth',2);grid on;
guidata(hObject,handles);
% --- Executes on button press in R_detect.
function R_detect_Callback(hObject, eventdata, handles)
% hObject handle to R_detect (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
denoised_signal=handles.denoised_signal;
[RR,r1,t1]=R_detect(denoised_signal);
axes(handles.axes3);
hold on;
plot(RR(:,1)./500,RR(:,2),'c+','MarkerEdgeColor','r','MarkerSize',16);
set(handles.R_R,'String',t1);
set(handles.heart_rate,'String',r1);
guidata(hObject,handles);
% --- Executes on button press in clear_all.
function clear_all_Callback(hObject, eventdata, handles)
% hObject handle to clear_all (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
set(handles.frequence,'String',' ');
set(handles.amplitude,'String',' ');
set(handles.edit8,'String',' ');
set(handles.edit11,'String',' ');
set(handles.SNR1,'String',' ');
set(handles.MSE1,'String',' ');
set(handles.SNR2,'String',' ');
set(handles.MSE2,'String',' ');
set(handles.heart_rate,'String',' ');
set(handles.R_R,'String',' ');
set(handles.popupmenu1,'Value',1);
set(handles.radiobutton1,'Value',0);
set(handles.radiobutton2,'Value',1);
set(handles.radiobutton3,'Value',0);
cla(handles.axes1);
cla(handles.axes2);
cla(handles.axes3);
cla(handles.axes4);
cla(handles.axes5);
cla(handles.axes6);
clear all;
三、运行结果
四、matlab版本及参考文献
1 matlab版本
2014a
2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.
[4]张泾周,张光磊,戴冠中.自适应算法与小波变换在心电信号滤波中的应用[J].生物医学工程学杂志。2006,(05)
文章出处登录后可见！
