Matlab心电信号预处理
Matlab心电信号预处理
一、内容
在网上下载心电信号的公开数据库，实现对心电信号的预处理，包括噪声去除、肌电干扰的去除、工频干扰的抑制、基线漂移的纠正等。
二、实验原理
1、肌电干扰的滤除
肌电干扰由人体肌肉颤动引起，发生率具有随机性，频率范围为20-500Hz，其主要成分的频率与肌肉的类型有关，一般在30-300Hz，而心电信号的频率主要集中在5~20HZ，所以选择低通滤波器来滤除肌电干扰。
巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波得图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。
2、工频干扰的抑制
工频干由于供电网络无所不在，因此50Hz的工频干扰是最普遍的，也是心电信号的主要干扰来源。50HZ陷波器的软件设计方法多种多样，常见方法有小波变换滤波、自适应滤波、模板匹配滤波等，但都需要手工计算获得滤波器的参数，运算比较复杂。
巴特沃斯带阻滤波器又称简单整系数带阻滤波器，其原理为一个全通网络，减去一个具有相同延迟和增益的窄带线性相位FIR滤波器，得到一个具有尖锐陷波特性的陷波滤波器。阻带下限截止频率fc1 = 49 Hz，阻带上限截止频率fc2 = 51 Hz，就可以消除50 Hz 的工频干扰。
3、基线漂移的纠正
某些数字信号中会含有基线干扰信号（低频噪音）。因此可以使用低通滤波器首先提取出低频噪音，然后再用原始信号减去低频噪音就可以得到去除了基线漂移的心电信号。
三、实现过程
1、读取ECG信号
MIT-BIH数据库中每个记录是从24小时的心电记录中选取了具有代表性的30分钟左右的数据片段组成数据库。每条一记录包括2导联数据,整个数据库大约有109494拍,都是经过至少两个以上的心电图专家手工独立标注的;记录以每通道每秒360采样点的规格进行数字化，具有11位分辨率，表示范围超过10毫伏。每个数据库记录包括三个文件,分别是头文件扩展名为.hea、数据文件扩展名为.dat和注释文件扩展名为.atr。
（1）头文件的识读
头文件的内容是由一行或多行ASCII码字符组成，并且至少包含-一个记录行,通常还有信号技术规范行、片段技术规范行(对于多片段数据记录)和信息注释行。
（2）数据文件的识读
MIT-BIH数据库中的数据存储格式有Format8、Format16、Format80、 Format212、 Format310 等8种，每一种格式中都是将来自两个或多个信号采样得到的数据交替存储。本文采用的是用于两个信号的采样数据交替存储的212格式,即用12位的二进制数据来存储，但标准计算机数据格式没有12 位的，就用3个8位的二进制数来表示2个12位的数。使用的时候只需用二进制读取文件的方式，再把数据整合就得到可用数据流，在Matlab平台下进行描图等一系列操作。
（3）注释文件的识读
记录了心电专家对相应心电信号的诊断信息，格式为MIT.它是一种紧簇性格式,每个注释的长度占用偶数个字节的空间，多数情况下是占用两个字节。每一个注释单元的前两个字节的第一个字节为最低有效位，16位中的最高6位表示了注释类型代码，剩余10位说明了该注释点的发生时间(其值为该注释点到前一注释点的间隔)或为辅助信息(附加信息的长度)。
数据文件读取代码：
%------ 具体数据------------------------------------------------------
PATH= 'D:\Matlab2016\bin\ECG'; % path, 这里就是写刚才你保存的数据地址
HEADERFILE= '100.hea'; % 文件格式为文本格式
ATRFILE= '100.atr'; % attributes-file 文件以二进制格式
DATAFILE='100.dat'; % data-file
SAMPLES2READ=660000; % 读取的数据样本点数为660000
% in case of more than one signal:
% 2*SAMPLES2READ samples are read
%------ LOAD HEADER DATA --------------------------------------------------
signald= fullfile(PATH, DATAFILE); % data in format 212
fid2=fopen(signald,'r');
A= fread(fid2, [3, SAMPLES2READ], 'uint8')'; % matrix with 3 rows, each 8 bits long, = 2*12bit
fclose(fid2);
%=----------------------载入二进制数据-----------------------------------------
M2H= bitshift(A(:,2), -4); %字节向右移四位，即取字节的高四位
M1H= bitand(A(:,2), 15); %取字节的低四位
M( : , 1)= bitshift(M1H,8)+ A(:,1); %低四位向左移八位
M( : , 2)= bitshift(M2H,8)+ A(:,3); %高四位向左移八位
M = M-1024; %这个M就是咱们解码出来的数据
m=0.005*(M);
t = (1:1600)/360;
ecg=m(1:1600,1);
plot (t, m(1:1600,1)) %绘制一段心电图形
2、去除肌电干扰
肌电干扰属于高频干扰，设计选用巴特沃斯数字低通滤波器。具体步骤是首先设计巴特沃斯型的模拟低通滤波器,再将该模拟滤波器转换为数字滤波器。设计用到的技术指标:心电信号通带截止频率wp;阻带截止频率ws,;通带最大衰减系数rp, ;阻带最小衰减系数rs。对于肌电干扰，rp=1.4，rs=1.6。
%------------------------------低通滤波器滤除肌电信号--------------------------
Fs=1500; %采样频率
fp=80;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,bp为分子系数，ap为分母系数
[bs,as]=lp2lp(bp,ap,wp) %Ha转换为低通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);
figure
freqz(bz,az);title('巴特沃斯低通滤波器幅频曲线');
figure
subplot(2,1,1);
%TIME = range(length(data));
TIME=1:length(ecg);
plot(t,ecg);
xlabel('时间（s）');ylabel('电压（mv）');title('原始心电信号波形');grid;
subplot(2,1,2);
plot(t,m);
xlabel('时间（s）');ylabel('电压（mv）');title('低通滤波后的时域图形');grid;
N=1500;
n=0:N-1;
mf=fft(ecg,N); %进行频谱变换（傅里叶变换）
mag=abs(mf);
f=(0:length(mf)-1)*Fs/length(mf); %进行频率变换
figure
subplot(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=ecg;
p=10*log10(abs(fft(wn).^2)/N);
f=(0:length(p)-1)/length(p); figure
plot(f,p);grid
xlabel('归一化频率');ylabel('功率(dB)');title('心电信号的功率谱');
3、去除工频干扰
在设计数字滤波器之前，根据实际的应用情况及滤波器的复杂程度，需要确定ECG信号工频滤波的技术指标。应用dlsim函数求解离散系统的响应。y=dlsim(b，a, x); 求输入信号为x时系统的响应。b和a分别表示系统函数H(z)中，由对应的分子项和分母项系数所构成的数组。陷波器系统函数为：
%数字陷波器
f0=50;
n=0:N-1;
%陷波器的设计
apha=-2*cos(2*pi*f0/Fs);
beta=0.96;
b=[1 apha 1];
a=[1 apha*beta beta^2];
figure(1);
freqz(b,a,N,Fs);%陷波器特性显示
y=dlsim(b,a,m);%陷波器滤波处理
%对信号进行频域变换。
y1=fft(y,N);
fa1=(0:length(y1)-1)*Fs/length(y1);
y2=y1.*conj(y1)/N;
figure(2);%滤除前后的信号对比。
subplot(2,2,1);plot(t,m);grid;
xlabel('时间（s）');ylabel('电压（mv）');title('原始心电信号波形');
subplot(2,2,3);plot(t,y);grid;
xlabel('时间（s）');ylabel('电压（mv）');title('50Hz陷波器滤波后的时域图形');
subplot(2,2,2);plot(fa1,maga);axis([0,100,1,100]);grid;
xlabel('频率(HZ)');ylabel('幅值');title('原始心电信号频谱图');
subplot(2,2,4);plot(fa1,abs(y1));axis([0,100,1,100]);grid;
xlabel('频率(HZ)');ylabel('幅值');title('50Hz陷波器滤波后频谱图');
3、基线漂移的矫正
fmaxd_1=5;%截止频率为5Hz
fmaxn_1=fmaxd_1/(Fs/2);
[B,A]=butter(1,fmaxn_1,'low');
figure(4)
freqz(B,A);
ecg_low=filtfilt(B,A,y);%通过5Hz低通滤波器的信号
ecg1=y-ecg_low; %去除这一段信号，得到去基线漂移的信号`在这里插入代码片
ecg2=fft(ecg1,N);
figure(5)
subplot(211)
plot(t,y);
xlabel('时间（s）');ylabel('电压（mv）');title('原始心电信号波形');
subplot(212)
plot(t,ecg1);
xlabel('时间（s）');ylabel('电压（mv）');title('去除基线漂移后的ECG数据');
subplot(2,1,1);plot(fa1,abs(y1));axis([0,100,1,100]);grid;
xlabel('频率(HZ)');ylabel('幅值');title('原始心电信号频谱图');
subplot(2,1,2);plot(fa1,abs(ecg2));axis([0,100,1,100]);grid;
xlabel('频率(HZ)');ylabel('幅值');title('去除基线漂移频谱图');
