ECG学习日记（一）

来源：互联网发布：用c语言表白编辑：程序博客网时间：2024/06/05 20:07

本文是针对学习过程中对ECG相关知识及源码的整合，现有的关于ECG网络文章比较少，希望通过写写该系列拙文简单记录下学习中的点滴吧。。。

参考博客：
1.http://blog.csdn.net/chenyusiyuan/article/details/2040234
2.http://blog.csdn.net/i_weimoli/article/details/53497384

ECG数据库

目前入门的用来做实验的最广泛的是MIT-BIH数据库 —— [ 下载地址 ]

MIT-BIH数据专门的Matlab心电数据读取程序:

clc; clear all;%------ SPECIFY DATA ------------------------------------------------------%------ 指定数据文件 -------------------------------------------------------PATH= 'D:\mycode\ECGData';  % 指定数据的储存路径HEADERFILE= '117.hea';      % .hea 格式，头文件，可用记事本打开ATRFILE= '117.atr';         % .atr 格式，属性文件，数据格式为二进制数DATAFILE='117.dat';         % .dat 格式，ECG 数据SAMPLES2READ=2048;          % 指定需要读入的样本数                            % 若.dat文件中存储有两个通道的信号:                            % 则读入 2*SAMPLES2READ 个数据 %------ LOAD HEADER DATA --------------------------------------------------%------ 读入头文件数据 -----------------------------------------------------%% 示例：用记事本打开的117.hea 文件的数据%%      117 2 360 650000%      117.dat 212 200 11 1024 839 31170 0 MLII%      117.dat 212 200 11 1024 930 28083 0 V2%      # 69 M 950 654 x2%      # None%%-------------------------------------------------------------------------fprintf(1,'\\n$> WORKING ON %s ...\n', HEADERFILE); % 在Matlab命令行窗口提示当前工作状态% % 【注】函数 fprintf 的功能将格式化的数据写入到指定文件中。% 表达式：count = fprintf(fid,format,A,...)% 在字符串'format'的控制下，将矩阵A的实数数据进行格式化，并写入到文件对象fid中。该函数返回所写入数据的字节数 count。% fid 是通过函数 fopen 获得的整型文件标识符。fid=1，表示标准输出（即输出到屏幕显示）；fid=2，表示标准偏差。%signalh= fullfile(PATH, HEADERFILE);    % 通过函数 fullfile 获得头文件的完整路径fid1=fopen(signalh,'r');    % 打开头文件，其标识符为 fid1 ，属性为'r'--“只读”z= fgetl(fid1);             % 读取头文件的第一行数据，字符串格式A= sscanf(z, '%*s %d %d %d',[1,3]); % 按照格式 '%*s %d %d %d' 转换数据并存入矩阵 A 中nosig= A(1);    % 信号通道数目sfreq=A(2);     % 数据采样频率clear A;        % 清空矩阵 A ，准备获取下一行数据for k=1:nosig           % 读取每个通道信号的数据信息    z= fgetl(fid1);    A= sscanf(z, '%*s %d %d %d %d %d',[1,5]);    dformat(k)= A(1);           % 信号格式; 这里只允许为 212 格式    gain(k)= A(2);              % 每 mV 包含的整数个数    bitres(k)= A(3);            % 采样精度（位分辨率）    zerovalue(k)= A(4);         % ECG 信号零点相应的整数值    firstvalue(k)= A(5);        % 信号的第一个整数值 (用于偏差测试)end;fclose(fid1);clear A;%------ LOAD BINARY DATA --------------------------------------------------%------ 读取 ECG 信号二值数据 ----------------------------------------------%if dformat~= [212,212], error('this script does not apply binary formats different to 212.'); end;signald= fullfile(PATH, DATAFILE);            % 读入 212 格式的 ECG 信号数据fid2=fopen(signald,'r');A= fread(fid2, [3, SAMPLES2READ], 'uint8')';  % matrix with 3 rows, each 8 bits long, = 2*12bitfclose(fid2);% 通过一系列的移位（bitshift）、位与（bitand）运算，将信号由二值数据转换为十进制数M2H= bitshift(A(:,2), -4);        %字节向右移四位，即取字节的高四位M1H= bitand(A(:,2), 15);          %取字节的低四位PRL=bitshift(bitand(A(:,2),8),9);     % sign-bit   取出字节低四位中最高位，向右移九位PRR=bitshift(bitand(A(:,2),128),5);   % sign-bit   取出字节高四位中最高位，向右移五位M( : , 1)= bitshift(M1H,8)+ A(:,1)-PRL;M( : , 2)= bitshift(M2H,8)+ A(:,3)-PRR;if M(1,:) ~= firstvalue, error('inconsistency in the first bit values'); end;switch nosigcase 2    M( : , 1)= (M( : , 1)- zerovalue(1))/gain(1);    M( : , 2)= (M( : , 2)- zerovalue(2))/gain(2);    TIME=(0:(SAMPLES2READ-1))/sfreq;case 1    M( : , 1)= (M( : , 1)- zerovalue(1));    M( : , 2)= (M( : , 2)- zerovalue(1));    M=M';    M(1)=[];    sM=size(M);    sM=sM(2)+1;    M(sM)=0;    M=M';    M=M/gain(1);    TIME=(0:2*(SAMPLES2READ)-1)/sfreq;otherwise  % this case did not appear up to now!    % here M has to be sorted!!!    disp('Sorting algorithm for more than 2 signals not programmed yet!');end;clear A M1H M2H PRR PRL;fprintf(1,'\\n$> LOADING DATA FINISHED \n');%------ LOAD ATTRIBUTES DATA ----------------------------------------------atrd= fullfile(PATH, ATRFILE);      % attribute file with annotation datafid3=fopen(atrd,'r');A= fread(fid3, [2, inf], 'uint8')';fclose(fid3);ATRTIME=[];ANNOT=[];sa=size(A);saa=sa(1);i=1;while i<=saa    annoth=bitshift(A(i,2),-2);    if annoth==59        ANNOT=[ANNOT;bitshift(A(i+3,2),-2)];        ATRTIME=[ATRTIME;A(i+2,1)+bitshift(A(i+2,2),8)+...                bitshift(A(i+1,1),16)+bitshift(A(i+1,2),24)];        i=i+3;    elseif annoth==60        % nothing to do!    elseif annoth==61        % nothing to do!    elseif annoth==62        % nothing to do!    elseif annoth==63        hilfe=bitshift(bitand(A(i,2),3),8)+A(i,1);        hilfe=hilfe+mod(hilfe,2);        i=i+hilfe/2;    else        ATRTIME=[ATRTIME;bitshift(bitand(A(i,2),3),8)+A(i,1)];        ANNOT=[ANNOT;bitshift(A(i,2),-2)];   end;   i=i+1;end;ANNOT(length(ANNOT))=[];       % last line = EOF (=0)ATRTIME(length(ATRTIME))=[];   % last line = EOFclear A;ATRTIME= (cumsum(ATRTIME))/sfreq;ind= find(ATRTIME <= TIME(end));ATRTIMED= ATRTIME(ind);ANNOT=round(ANNOT);ANNOTD= ANNOT(ind);%------ DISPLAY DATA ------------------------------------------------------figure(1); clf, box on, hold onplot(TIME, M(:,1),'r');if nosig==2    plot(TIME, M(:,2),'b');end;for k=1:length(ATRTIMED)    text(ATRTIMED(k),0,num2str(ANNOTD(k)));end;xlim([TIME(1), TIME(end)]);xlabel('Time / s'); ylabel('Voltage / mV');string=['ECG signal ',DATAFILE];title(string);fprintf(1,'\\n$> DISPLAYING DATA FINISHED \n');% -------------------------------------------------------------------------fprintf(1,'\\n$> ALL FINISHED \n');

matlab2016b运行上述程序，可以得到该数据下ECG：
这里写图片描述

心电图干扰

心电图的干扰因素主要有以下三种：

工频干扰
基线漂移
肌电干扰

在读取心电图数据后，肌电信号的滤除的matlab源码：

clc;%------------------------------低通滤波器滤除肌电信号------------------------------Fs=1500;                        %采样频率fp=80;fs=100;                    %通带截止频率，阻带截止频率rp=1.4;rs=1.6;                    %通带、阻带衰减wp=2*pi*fp;ws=2*pi*fs;   [n,~]=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,M(:,1));figurefreqz(bz,az);title('巴特沃斯低通滤波器幅频曲线');figuresubplot(2,1,1);plot(TIME,M(:,1));xlabel('t(s)');ylabel('mv');title('原始心电信号波形');grid;subplot(2,1,2);plot(TIME,m);xlabel('t(s)');ylabel('mv');title('低通滤波后的时域图形');grid;N=512;n=0:N-1;mf=fft(M(:,1),N);               %进行频谱变换（傅里叶变换）mag=abs(mf);f=(0:length(mf)-1)*Fs/length(mf);  %进行频率变换figuresubplot(2,1,1)plot(f,mag);axis([0,1500,1,50]);grid;      %画出频谱图xlabel('频率(HZ)');ylabel('幅值');title('心电信号频谱图');mfa=fft(m,N);                    %进行频谱变换（傅里叶变换）maga=abs(mfa);fa=(0:length(mfa)-1)*Fs/length(mfa);  %进行频率变换subplot(2,1,2)plot(fa,maga);axis([0,1500,1,50]);grid;  %画出频谱图xlabel('频率(HZ)');ylabel('幅值');title('低通滤波后心电信号频谱图');wn=M(:,1);P=10*log10(abs(fft(wn).^2)/N);f=(0:length(P)-1)/length(P);figureplot(f,P);gridxlabel('归一化频率');ylabel('功率(dB)');title('心电信号的功率谱');

运行得到处理后的结果：

接着是对工频干扰抑制的matlab源码：

%-----------------带陷滤波器抑制工频干扰-------------------%50Hz陷波器：由一个低通滤波器加上一个高通滤波器组成%而高通滤波器由一个全通滤波器减去一个低通滤波器构成Me=100;               %滤波器阶数L=100;                %窗口长度beta=100;             %衰减系数Fs=1500;wc1=49/Fs*pi;     %wc1为高通滤波器截止频率，对应51Hzwc2=51/Fs*pi     ;%wc2为低通滤波器截止频率，对应49Hzh=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陷波器冲击响应序列m2=filter(y,1,m);figuresubplot(2,1,1);plot(abs(h));axis([0 100 0 0.2]);xlabel('频率(Hz)');ylabel('幅度(mv)');title('陷波器幅度谱');grid;N=512;P=10*log10(abs(fft(y).^2)/N);f=(0:length(P)-1);subplot(2,1,2);plot(f,P);xlabel('频率(Hz)');ylabel('功率(dB)');title('陷波器功率谱');grid;figuresubplot (2,1,1); plot(TIME,m);xlabel('t(s)');ylabel('幅值');title('原始信号');grid;subplot(2,1,2);plot(TIME,m2);xlabel('t(s)');ylabel('幅值');title('带阻滤波后信号');grid;figureN=512;subplot(2,1,1);plot(abs(fft(m))*2/N);axis([0 100 0 1]);xlabel('t(s)');ylabel('幅值');title('原始信号频谱');grid;subplot(2,1,2);plot(abs(fft(m2))*2/N);axis([0 100 0 1]);xlabel('t(s)');ylabel('幅值');title('带阻滤波后信号频谱');grid;  %理想低通滤波器%截止角频率wc，阶数Mefunction hd=ideal_lp(wc,Me)alpha=(Me-1)/2;n=[0:Me-1];p=n-alpha+eps;              %eps为很小的数，避免被0除hd=sin(wc*p)./(pi*p);       %用Sin函数产生冲击响应end

运行结果为：

进一步接着基线漂移的纠正的matlab源码：

%------------------IIR零相移数字滤波器纠正基线漂移-------------------Wp=1.4*2/Fs;     %通带截止频率 Ws=0.6*2/Fs;     %阻带截止频率 devel=0.005;    %通带纹波 Rp=20*log10((1+devel)/(1-devel));   %通带纹波系数  Rs=20;                          %阻带衰减 [N Wn]=ellipord(Wp,Ws,Rp,Rs,'s');   %求椭圆滤波器的阶次 [b a]=ellip(N,Rp,Rs,Wn,'high');       %求椭圆滤波器的系数 [hw,w]=freqz(b,a,512);   result =filter(b,a,m2); figurefreqz(b,a);figuresubplot(211); plot(TIME,m2); xlabel('t(s)');ylabel('幅值');title('原始信号');gridsubplot(212); plot(TIME,result); xlabel('t(s)');ylabel('幅值');title('线性滤波后信号');gridfigureN=512subplot(2,1,1);plot(abs(fft(m2))*2/N);xlabel('频率(Hz)');ylabel('幅值');title('原始信号频谱');grid;subplot(2,1,2);plot(abs(fft(result))*2/N);xlabel('频率(Hz)');ylabel('幅值');title('线性滤波后');grid;ubplot(2,1,2);plot(abs(fft(result))*2/N);xlabel('线性滤波后信号频谱');ylabel('幅值');grid;

运行结果为：

此外，运用小波来对ECG进行处理是一种行之有效的方法，代码如下：

clc;E=M(:,2);E=E';n=size(E);s=E(1:2000);%小波分解[C L]=wavedec(E,3,'db5');cD1=detcoef(C,L,1);cD2=detcoef(C,L,2);cD3=detcoef(C,L,3);thr1=thselect(cD1,'rigrsure');thr2=thselect(cD2,'rigrsure');thr3=thselect(cD3,'rigrsure');TR=[thr1,thr2,thr3];SORH='s';[XC,CXC,LXC,PERF0,PERF2]=wdencmp('lvd',E,'db5',3,TR,SORH);N=n(2);x=E;y=XC;F=0;M=0;for ii=1:N    m(ii)=(x(ii)-y(ii))^2;    t(ii)=y(ii)^2;    f(ii)=t(ii)/m(ii);    F=F+f(ii);    M=M+m(ii);end;SNR=10*log10(F);MSE=M/N;SM=SNR/MSE;subplot(2,1,1);plot(s(1:1000));title('原始信号')subplot(2,1,2);plot(XC(1:1000));title('处理后的信号')SNR,MSE

运行效果：

阅读全文

0 0