第一句子网 > 【心电信号】基于matlab GUI心电信号数字滤波处理【含Matlab源码 484期】

【心电信号】基于matlab GUI心电信号数字滤波处理【含Matlab源码 484期】

时间：2024-03-23 19:44:22

一、心电信号数字滤波处理简介

心电信号作为一种人体的基本生理信号, 是心脏电活动在人体体表的表现, 信号一般比较微弱, 频率在0.05Hz~100Hz范围内, 幅度为10V (胎儿) ~5m V (成人) , 心电信号信噪比和频率都较低, 在心电的采集、放大、检测等过程中, 易受到外界的各种干扰。常见的噪声干扰有:第一是基线漂移, 一般是由人体呼吸和心肌兴奋所引起的, 它的频率低于0.5Hz, 属于低频干扰;其次是肌电干扰, 它是由人体肌肉颤动所致, 它的发生频率具有随机性, 范围在5Hz~2000Hz之间;第三是工频干扰, 它是由室内照明及动力设备影响到人体的分布电容所引起的, 频率为50Hz。消除或减少这些干扰时识别心电信号特征和参数的前提。心电信号噪声来源不同, 频率也存在差异, 正是由于这些差异, 对不同的信号干扰其滤波方法也不同, 滤波可以用硬件实现, 但实现过程相对困难, 也可以用软件编程方法实现, 数字滤波技术成为目前滤除心电干扰的有效手段。

1 程序设计与实现

“心电信号的数字滤波处理”软件所要实现的功能和任务如下所示:

(1) 信号输入:信号源的读取及参数的输入;

(2) 信号滤波:选择信号分析通道, 选择滤波器类型和种类, 根据选择的滤波器类型及参数指标, 用相对应的阶数选择函数返回阶数N及截止频率Wn, 根据N及Wn利用IIR滤波器响应的设计函数对信号进行滤波处理。

(3) 滤波器特性演示:显示所设计的各个滤波器的幅度和相位响应;

(4) 信号显示:包括原始信号的显示和经过各次滤波后信号的显示。

图1 程序设计流程图

2 具体界面设计如下所示

(1) 打开并选择文件:创建打开文件对话框并显示文件存储路径, 数据文件通常为.txt或.dat格式。“选择文件”用按钮 (push button) 实现, 当点击时, 能够打开如图4所示的对话框, 可选择数据文件, 并在文本框 (edit text) 内可显示文件存储路径。

(2) 信号通道选择:由于所采集的心电信号数据是12通道的, 进行心电信号分析时只需选择其中之一, 信号选择通道用下拉菜单 (pop-up menu) 实现。

(3) 选择滤波器类型:滤波器类型共有四种:Butterworth、Chebyshev1、Chebyshev2、Elliptic, 用下拉菜单 (pop-up menu) 实现。

(4) 滤波器功能实现:带阻、高通、低通分别用三个单选按钮 (radio button) 实现, 并用按钮组 (button group) 把三个控件组织在同一区域内。

(5) 参数输入:采样频率 (Fs) 、通带截止频率 (Fp1, Fp2) 、阻带截止频率 (Fs1, Fs2) 、通带波动 (Rp) 、阻带衰减 (Rs) 从界面上输入, 显示这些参数指标的组件为文本标签 (static text) , 显示输入参数的组件为文本框 (edit text) , 当选择带阻时能显示全部参数指标, 当选择低通、高通时, 通过设置属性“visible”, 隐藏通带截止频率 (Fp2) 、阻带截止频率 (Fs2) , 使其编辑框不可见。

(6) 滤波结果显示:滤波前后的波形在相应的坐标轴上显示, 其中第一个坐标轴显示滤波前的原始心电信号signal, 第二个坐标轴显示signal经过带阻滤波后的信号波形s1, 第三个坐标轴显示s1经高通滤波后的信号波形s2, 第四个坐标轴显示s2经过低通滤波后的波形s3, 第五个坐标轴同时显示原始信号signal和经过三次滤波后的波形s3, 以便进行对比分析, 点击每个坐标轴旁边的“Run”按钮可显示坐标轴相对应的信号波形。

(7) 滤波器性能演示及退出界面:用按钮 (push button) 实现, 当点击“滤波器性能演示”时, 能够显示选择的滤波器的幅度及相位响应, 当点击“Quit”时, 退出当前窗口, 返回编辑界面。

二、部分源代码

function varargout = ECG(varargin)% ECG MATLAB code 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 22-May- 21:35:12% Begin initialization code - DO NOT EDITgui_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});endif nargout[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});elsegui_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 ECGhandles.output = hObject;% Update handles structureguidata(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 structurevarargout{1} = handles.output;function edit1_Callback(hObject, eventdata, handles)% hObject handle to edit1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of edit1 as text% str2double(get(hObject,'String')) returns contents of edit1 as a double% --- Executes during object creation, after setting all properties.function edit1_CreateFcn(hObject, eventdata, handles)% hObject handle to edit1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');end% --- Executes on button press in pushbutton1.function pushbutton1_Callback(hObject, eventdata, handles)% hObject handle to pushbutton1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)[FileName PathName]=uigetfile('*.mat','MAT-files(*.mat)','选择文件');str=[PathName FileName];set(handles.edit1,'string',str);Fs=200;N=512;MEAN=0;MIN=0;MAX=0;VAR=0;STD=0;RR=0;global im;global MEAN;global MIN;global MAX;global VAR;global STD;global RR;if strcmp(str,'D:\心电数据\被试2 心电\chenwei1.mat')==1im=1;elseim=2;endqq=get(handles.popupmenu1,'value')bb=get(handles.popupmenu2,'value')switch(im)case 1load('chenwei1.mat');ECG=xin(qq,:);case 2load('fanglipeng1.mat');ECG=xin(qq,:);end% uu1=get(handles.radiobutton7,'value');% uu2=get(handles.radiobutton8,'value');% uu3=get(handles.radiobutton8,'value');% if uu1==1%MAX=max(EGC);%set(handles.edit3,'string',MAX);% end% if uu1==2%MIN=min(EGC);%set(handles.edit4,'string',MIN);% end% if uu1==3%MEAN=mean(EGC); %set(handles.edit5,'string',MEAN);% endifbb==1disp([num2str(bb)]);t=1/50:1/50:length(ECG)/50;MEAN=mean(xin(qq,:));MAX=max(xin(qq,:)); MIN=min(xin(qq,:));VAR=var(xin(qq,:),0,2);STD=std(xin(qq,:),0,2);axes(handles.axes3);grid on;plot(t,ECG);xlabel('时间（s）');ylabel('幅值(V)');title('原始心电信号时域显示');N=512;y=fft(ECG,N);mag=abs(y);f=(0:length(y)-1)*Fs/length(y);axes(handles.axes5);plot(f,mag)xlim([0,100]);title('原始心电信号频谱图')xlabel('频率（Hz）');ylabel('幅值');axes(handles.axes8);histogram(ECG);xlabel('信号幅值');ylabel('个数');title('不同区间心电信号的分布');elseifbb==2disp([num2str(bb)]);%-----------------------------心电信号高通滤波器进行基线漂移纠正--------------------------------%%还可采用椭圆滤波器或中值法（myMedfilt）等进行基线漂移纠正。fp=1.4;fr=0.05;rp=1;rs=30;wp=fp/(Fs/2);wr=fr/(Fs/2);[n,wn]=buttord(wp,wr,rp,rs);[b,a]=butter(n,wn,'high');[hw,w]=freqz(b,a);y2=filter(b,a,ECG);axes(handles.axes1);grid on;t1=0:1/50:(length(y2)-1)/50;plot(t1,y2)title('基线漂移去除时域')xlabel('时间（s）');ylabel('幅值(V)');y2=fft(y2,N);mfb=abs(y2);f=(0:length(y2)-1)*Fs/length(y2);axes(handles.axes2);grid on;plot(f,mfb);xlim([0,100]);title('基线纠正后信号幅频')xlabel('频率（Hz）');ylabel('幅值');elseifbb==3disp([num2str(bb)]);fp=60;fs=100;%通带截止频率，阻带截止频率 rp=1.4;rs=1.6; %通带、阻带衰减 wp=2*pi*fp;ws=2*pi*fs;[n,wn]=buttord(wp,ws,rp,rs,'s'); %’s’是确定巴特沃斯模拟滤波器阶次和3dB 截止模拟频率 [z,P,k]=buttap(n); %设计归一化巴特沃斯模拟低通滤波器，z为极点，p为零点和k为增益 [bp,ap]=zp2tf(z,P,k); %转换为Ha(p),bp为分子系数，ap为分母系数 [bs,as]=lp2lp(bp,ap,wp); %Ha(p)转换为低通Ha(s)并去归一化，bs为分子系数，as为分母系数 [hs,ws]=freqs(bs,as); %模拟滤波器的幅频响应 [bz,az]=bilinear(bs,as,Fs);%对模拟滤波器双线性变换 [h1,w1]=freqz(bz,az); %数字滤波器的幅频响应 m=filter(bz,az,ECG);t1=0:1/50:(length(m)-1)/50;axes(handles.axes1);grid on;axis tight;plot(t1,m)title('低通滤波心电信号时域图');xlabel('时间（s）');ylabel('幅值(V)');y1=fft(m,N);mfa=abs(y1);f=(0:length(y1)-1)*Fs/length(y1);axes(handles.axes2);grid on;axis tight;plot(f,mfa);xlim([0,100]);title('低通滤波心电信号频谱图');xlabel('频率（Hz）');ylabel('幅值');elseifbb==4 disp([num2str(bb)]);%50Hz陷波器：由一个低通滤波器加上一个高通滤波器组成 %而高通滤波器由一个全通滤波器减去一个低通滤波器构成 Me=100;%滤波器阶数 L=100;%窗口长度 beta=100; %衰减系数 Fs=200; wc1=49/Fs*pi;%wc1为高通滤波器截止频率，对应51Hz wc2=50/Fs*pi;%wc2为低通滤波器截止频率，对应49Hz h=ideal_lp(0.132*pi,Me)-ideal_lp(wc1,Me)+ideal_lp(wc2,Me); %h为陷波器冲击响应 w=kaiser(L,beta); y=h.*rot90(w); %y为50Hz陷波器冲击响应序列 y3=filter(y,1,ECG); t1=0:1/50:(length(y3)-1)/50;axes(handles.axes1);grid on;axis tight;plot(t1,y3)title('去除工频干扰时域时域')xlabel('时间（s）');ylabel('幅值(V)');y3=fft(y3,N);mfc=abs(y3);f=(0:length(y3)-1)*Fs/length(y3);axes(handles.axes2);grid on;axis tight;plot(f,mfc);xlim([0,100]);title('去除工频干扰时域频域');xlabel('频率（Hz）');ylabel('幅值');end%%%心率计算%低通滤波fp=60;fs=100;%通带截止频率，阻带截止频率 rp=1.4;rs=1.6; %通带、阻带衰减 wp=2*pi*fp;ws=2*pi*fs;[n,wn]=buttord(wp,ws,rp,rs,'s'); %’s’是确定巴特沃斯模拟滤波器阶次和3dB 截止模拟频率 [z,P,k]=buttap(n); %设计归一化巴特沃斯模拟低通滤波器，z为极点，p为零点和k为增益 [bp,ap]=zp2tf(z,P,k); %转换为Ha(p),bp为分子系数，ap为分母系数 [bs,as]=lp2lp(bp,ap,wp); %Ha(p)转换为低通Ha(s)并去归一化，bs为分子系数，as为分母系数 [hs,ws]=freqs(bs,as); %模拟滤波器的幅频响应 [bz,az]=bilinear(bs,as,Fs);%对模拟滤波器双线性变换 [h1,w1]=freqz(bz,az); %数字滤波器的幅频响应 m=filter(bz,az,ECG);fp=1.4;fr=0.05;rp=1;rs=30;wp=fp/(fs/2);wr=fr/(fs/2);[n,wn]=buttord(wp,wr,rp,rs);[b,a]=butter(n,wn,'high');[hw,w]=freqz(b,a);y2=filter(b,a,m);%高通滤波fp=1.4;fr=0.05;rp=1;rs=30;wp=fp/(Fs/2);wr=fr/(Fs/2);[n,wn]=buttord(wp,wr,rp,rs);[b,a]=butter(n,wn,'high');[hw,w]=freqz(b,a);y2=filter(b,a,m);%带阻滤波Me=100;%滤波器阶数 L=100;%窗口长度 beta=100; %衰减系数 wc1=49/Fs*pi;%wc1为高通滤波器截止频率，对应51Hz wc2=50/Fs*pi;%wc2为低通滤波器截止频率，对应49Hz h=ideal_lp(0.132*pi,Me)-ideal_lp(wc1,Me)+ideal_lp(wc2,Me); %h为陷波器冲击响应 w=kaiser(L,beta); y=h.*rot90(w); %y为50Hz陷波器冲击响应序列 y3=filter(y,1,y2); %------------------------R波检测及心率计算---------------------------------a=1;b=[-0.000563925539225382,-0.000687497480673608,-0.000400105459896818,1.40107308696153e-18,0.000167430218323421,-4.18540736387490e-19,-0.000215512204999462,-2.31378858508372e-19,0.000825087119402599,0.00173867121147225,0.00170041274003132,-2.89185229056478e-18,-0.00281247259984511,-0.00484192899664537,-0.00404159503274541,6.02125250728051e-18,0.00516947529617596,0.00792307431385663,0.00588976752279114,-7.41129393462523e-18,-0.00583817525454588,-0.00764316928687033,-0.00463943332803296,-9.61260744539564e-19,0.00196521978920439,-2.40802770032636e-18,-0.00240551785915901,4.71696771800100e-18,0.00853358817557591,0.0172851206658474,0.0162833499451329,-1.24036605494920e-17,-0.0253107164699268,-0.0426277900469897,-0.0350794664744733,1.73330580157676e-17,0.0449057473944898,0.0701489638765231,0.0540262508439817,-3.98392673063995e-17,-0.0616816640708770,-0.0915638898806427,-0.0672043113722621,-2.75356818775531e-17,0.0700978687033445,0.0996538686301192,0.0700978687033445,-2.75356818775531e-17,-0.0672043113722621,-0.0915638898806427,-0.0616816640708770,-3.98392673063995e-17,0.0540262508439817,0.0701489638765231,0.0449057473944898,1.73330580157676e-17,-0.0350794664744733,-0.0426277900469897,-0.0253107164699268,-1.24036605494920e-17,0.0162833499451329,0.0172851206658474,0.00853358817557591,4.71696771800100e-18,-0.00240551785915901,-2.40802770032636e-18,0.00196521978920439,-9.61260744539564e-19,-0.00463943332803296,-0.00764316928687033,-0.00583817525454588,-7.41129393462523e-18,0.00588976752279114,0.00792307431385663,0.00516947529617596,6.02125250728051e-18,-0.00404159503274541,-0.00484192899664537,-0.00281247259984511,-2.89185229056478e-18,0.00170041274003132,0.00173867121147225,0.000825087119402599,-2.31378858508372e-19,-0.000215512204999462,-4.18540736387490e-19,0.000167430218323421,1.40107308696153e-18,-0.000400105459896818,-0.000687497480673608,-0.000563925539225382];d=filter(b,a,y3);t=(0:(length(d)-1))/Fs;d1=sort(d,'descend');d2=0;for i=201 :2200d2=d1(i)+d2;end;d3=d2/2000;[R_V,R_L]=findpeaks(d(12001:24000),'minpeakdistance',1,'minpeakheight',d3);d3=d2/2000;[R_V,R_L]=findpeaks(d(12001:24000),'minpeakdistance',1,'minpeakheight',d3);%根据位置和采样频率来计算采样区间内的平均心率H_Rate = 60*(length(R_L)-1)/((R_L(length(R_L))-R_L(1))/200);RR=num2str(H_Rate,3);%算出采样的时间区间%R_Time = R_L(length(R_L))/200; % --- Executes on button press in radiobutton7.function radiobutton7_Callback(hObject, eventdata, handles)% hObject handle to radiobutton7 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)global MAX;set(handles.edit3,'string',MAX);% Hint: get(hObject,'Value') returns toggle state of radiobutton7% --- Executes on button press in radiobutton8.function radiobutton8_Callback(hObject, eventdata, handles)% hObject handle to radiobutton8 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)global MIN;set(handles.edit4,'string',MIN);% Hint: get(hObject,'Value') returns toggle state of radiobutton8% --- Executes on button press in radiobutton9.function radiobutton9_Callback(hObject, eventdata, handles)% hObject handle to radiobutton9 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)global MEAN;set(handles.edit5,'string',MEAN);% Hint: get(hObject,'Value') returns toggle state of radiobutton9function edit3_Callback(hObject, eventdata, handles)% hObject handle to edit3 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of edit3 as text% str2double(get(hObject,'String')) returns contents of edit3 as a double% --- Executes during object creation, after setting all properties.function edit3_CreateFcn(hObject, eventdata, handles)% hObject handle to edit3 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');end