基于matlab編程的心電信號特征提取以及分析，應用小波閾值去燥法，小波熵演算法，哪位大神幫幫忙編程出來，萬分感謝！！
uj5u.com熱心網友回復：
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)]);
%效果不好
uj5u.com熱心網友回復：
<br />樓主我想問一下你，基于小波變換的非平穩信號特征提取的代碼如何寫？轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/99485.html
標籤：軟件測試
下一篇：鴻蒙自定義組件之鴻蒙畫板