【心電訊號】基於matlab心電訊號ECG濾波處理【含Matlab原始碼 1067期】
阿新 • • 發佈:2021-06-28
一、簡介
基於matlab心電訊號ECG濾波處理
二、原始碼
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-2020 21:35:12 % 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; 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')==1 im=1; else im=2; end qq=get(handles.popupmenu1,'value') bb=get(handles.popupmenu2,'value') switch(im) case 1 load('chenwei1.mat'); ECG=xin(qq,:); case 2 load('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); % end if bb==1 disp([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('不同區間心電訊號的分佈'); elseif bb==2 disp([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('幅值'); elseif bb==3 disp([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('幅值'); elseif bb==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 :2200 d2=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 radiobutton9 function 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
三、執行結果
四、備註
版本:2014a
完整程式碼或代寫加1564658423