5,177
社区成员
发帖
与我相关
我的任务
分享基于matlab编程的心电信号特征提取以及分析
基于matlab编程的心电信号特征提取以及分析，应用小波阈值去燥法，小波熵算法，哪位大神帮帮忙编程出来，万分感谢！！
...全文
请发表友善的回复…
发表回复
SUNNYGU0127 2020-02-27
- 打赏
- 举报
<br />楼主我想问一下你，基于小波变换的非平稳信号特征提取的代码如何写？
qq_41554528 2017-12-28
- 打赏
- 举报
PATH= 'E:\xdsj'; %指定数据的储存路径
HEADERFILE= '109.hea'; %.hea 格式，头文件，可用记事本打开
DATAFILE='109.dat'; %.dat 格式，ECG 数据
SAMPLES2READ=2048; %指定需要读入的样本数,若.dat文件中存储有两个通道的信号:则读入 2*SAMPLES2READ 个数据
%头文件中的信息
fprintf(1,'\\n$> WORKING ON %s ...\n', HEADERFILE); % 在Matlab命令行窗口提示当前工作状态
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;
%读取data数据
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*12bit 矩阵A共有SAMPLES2READ行、3列，每列数据都是以uint8格式读入，注意这时数据通过uint8的读入方式已经成为十进制数了
fclose(fid2);
M2H= bitshift(A(2,:), -4); % 字节向右移四位，即取字节的高四位，属于信号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;% 将M1H、M2H分别左移8位，即乘以2^8，再分别加上A(:,1)，A(:,2)，
M( 2 , :)= bitshift(M2H,8)+ A( 2 , : )-PRR;% 由于左移时把符号位也移动了，要减去符号位的值
%if M(1,:) ~= firstvalue, error('inconsistency in the first bit values'); end;
%switch nosig
%case 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 转置为2行SAMPLES2READ列
%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');
figure(1); clf, box on, hold on
plot(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');
%心电去噪
qz=sgolayfilt(M(1,:),2,15);
figure(2)
plot(TIME, qz);
xlim([TIME(1), TIME(end)]);
%X=M(1,:);
X=qz;
% AF1
% ataf1=0.3*max(X);
% for n=2:2047
% B(n)=X(n+1)-X(n-1);
% end
% for m=1:36
%
% num0=0;%将满足条件的数的下标存入数组a1中，num0为数组的下标
% num00=0;%将不满足条件的数的下标存入数组b1中，num00为数组的下标
% for i=2:2011
% if(B(i)>0.05&&B(i+1)>0.05&&B(i+2)>0.05)
% for j=1:36
% if (B(i+j)<-0.03&&B(i+j+1)<-0.03)
% num0=num0+1;
% a1(num0)=i;
% break,
% end
%
% end
% end
% end
%AF1
ataf1=-0.1;
for n=2:2047
Y(n)=X(n+1)-X(n-1);
end
% AF1=[];
% NAF1=[];
% for n=2:2037
% if X(n)>ataf1&Y(n-1)>0.15&Y(n+1)<-0.15
% AF1=[AF1,n];
% else NAF1=[NAF1,n];
% end
% end
% figure(2)
% plot(TIME, qz,'r');
% hold on
% plot(TIME(AF1),X(AF1),'+');
% plot(TIME(NAF1),X(NAF1),'.');
%用Y的阈值来找到QRS波_yyp
yup=find(Y>0.05);
ydn=find(Y<-0.05);
figure(11)
plot(TIME, qz,'b');
hold on
plot(TIME(yup),X(yup),'k+');
plot(TIME(ydn),X(ydn),'r+');
%用X的阈值来找到R波_yyp
xhi=find(X>0.4*mean(X));
figure(12)
plot(TIME, qz,'b');
hold on
plot(TIME(xhi),X(xhi),'r+');
%AF2
ataf2=0.4*max(X);%阈值
for n=1:2048
if X(n)<0
Y0(n)=-X(n);
else
Y0(n)=X(n);
end
end
for i=1:2048
if Y0(i)<ataf2
Y1(i)=ataf2;
else
Y1(i)=Y0(i);
end
end
AF2=[]; %将满足条件的数的下标存入数组AF2中
NAF2=[]; %将不满足条件的数的下标存入数组NAF2中
for j=2:2047
Y2(j)=Y1(j+1)-Y1(j-1);
if Y2(j)>0.01
AF2=[AF2,j];
else
NAF2=[NAF2,j];
end
end
figure(3)
plot(TIME, qz,'r');
hold on
plot(TIME(AF2),X(AF2),'.');
xlim([TIME(1), TIME(end)]);
%问题：S波的右边检测不出来
%AF3
ataf3=0.02;
for n=2:2047
M(n)=X(n+1)-X(n-1);
end
AF3=[]; %将满足条件的数的下标存入数组AF3中
NAF3=[]; %将不满足条件的数的下标存入数组NAF3中
for i=2:2047
if (M(i)>=ataf3&&M(i)*X(i)>0)
AF3=[AF3,i];
elseif (M(i)<=-ataf3&&M(i)*X(i)<0)
AF3=[AF3,i];
else
NAF3=[NAF3,i];
end
end
figure(4);
plot(TIME, qz,'r');
hold on
plot(TIME(AF3),X(AF3),'.');
xlim([TIME(1), TIME(end)]);
%问题：几乎没有效果
%取顶点_yyp
START=[];
for i=1:length(AF3)-1
if((AF3(i+1)-AF3(i))>2) START=[START,i];
else;
end
end
START=[START 1];
START=sort(START);
PEAK=[];
for i=1:length(START)-1
TEMP=[];
TEMP=AF3(START(i):(START(i+1)-1));
[a b]=max(X(TEMP));
PEAK=[PEAK,TEMP(b)];
end
plot(TIME(PEAK),X(PEAK),'k+');
xlim([TIME(1), TIME(end)]);
%FD1
for n=3:2046
N(n)=-2*X(n-2)-X(n-1)+X(n+1)+2*X(n+2);
end
atfd1=0.07*max(N);
FD1=[]; %将满足条件的数的下标存入数组FD1中
NFD1=[]; %将不满足条件的数的下标存入数组NFD1中
for i=3:2046
if N(i)>atfd1
FD1=[FD1,i];
else
NFD1=[NFD1,i];
end
end
figure(5);
plot(TIME, qz,'r');
hold on
plot(TIME(FD1),X(FD1),'.');
xlim([TIME(1), TIME(end)]);
%问题：Q波的左边检测不出来
%FD2
atfd2=0.03;
for n=2:2047
O(n)=X(n+1)-X(n-1);
end
FD2=[]; %将满足条件的数的下标存入数组FD2中
NFD2=[]; %将不满足条件的数的下标存入数组NFD2中
for i=2:2044
if (O(i)>atfd2&&(O(i+1)>atfd2||O(i+2)>atfd2||O(i+3)>atfd2))
FD2=[FD2,i];
else
NFD2=[NFD2,i];
end
end
figure(6);
plot(TIME, qz,'r');
hold on
plot(TIME(FD2),X(FD2),'.');
xlim([TIME(1), TIME(end)]);
%问题：Q波的左边检测不出来
%FS1
for n=3:2046
P0(n)=abs(X(n+1)-X(n-1));
P1(n)=abs(X(n+2)-2*X(n)+X(n-2));
P2(n)=1.3*P0(n)+1.1*P1(n);
end
FS1=[]; %将满足条件的数的下标存入数组FS1中
NFS1=[]; %将不满足条件的数的下标存入数组NFS1中
for i=3:2046
if P2(i)>=0.06
FS1=[FS1,i];
else
NFS1=[NFS1,i];
end
end
figure(7);
plot(TIME, qz,'r');
hold on
plot(TIME(FS1),X(FS1),'.');
xlim([TIME(1), TIME(end)]);
%效果很好
%FS2
for n=2:2047
Q0(n)=abs(X(n+1)-X(n-1));
end
for i=3:2046
Q1(i)=(Q0(i+1)+2*Q0(i)+Q0(i-1))/4;
end
for j=4:2045
Q2(j)=abs(X(j+2)-2*X(j)+X(j-2));
Q3(j)=Q1(j)+Q2(j);
end
ptfs2=0.1*max(Q3);
stfs2=0.0001*max(Q3);
FS2=[]; %将满足条件的数的下标存入数组FS2中
NFS2=[]; %将不满足条件的数的下标存入数组NFS2中
for w=4:2039
if (Q3(w)>=ptfs2&&Q3(w+1)>=stfs2&&Q3(w+2)>=stfs2&&Q3(w+3)>=stfs2&&Q3(w+4)>=stfs2&&Q3(w+5)>=stfs2&&Q3(w+6)>=stfs2)
FS2=[FS2,w];
else
NFS2=[NFS2,w];
end
end
figure(8);
plot(TIME, qz,'r');
hold on
plot(TIME(FS2),X(FS2),'.');
xlim([TIME(1), TIME(end)]);
%效果比较好
%DF1
for n=5:2048
Z0(n)=X(n)-X(n-4);
end
for n=9:2048
Z1(n)=Z0(n)+4*Z0(n-1)+6*Z0(n-2)+4*Z0(n-3)+Z0(n-4);
end
numa=0;
for i=9:1991
if Z1(i)>0.05
for j=1:57
if Z1(i+j)<-0.05
numa=numa+1;
aa(numa)=i;
end
end
end
end
figure(9);
plot(TIME, qz,'r');
hold on
plot(TIME(aa),X(aa),'.');
xlim([TIME(1), TIME(end)]);
%效果不好
