-
Notifications
You must be signed in to change notification settings - Fork 0
/
ensayo1.m
194 lines (176 loc) · 6.5 KB
/
ensayo1.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
close all;
clear all;
addpath('C:\Program Files\MATLAB\R2018a\prtools')
addpath('C:\MATLAB2018\MATLAB\mcode\BioMedicine\Dream Abnea\BiomedicasFinal\apnea-ecg');
% m represents the size
m = 6000;
% NUMBER OF DATASETS
N = 7;
% STRINGS MEMORY BANK
names=cell(8,1);
for i = 1:4
dbNames = strcat({'a0'},int2str(i),{'erm.mat'});
names{i}=dbNames;
end
for i=1:3
string=strcat({'c0'},int2str(i),{'erm.mat'});
names{i+4}=string;
end
% SAVE IN PHYSICAL UNITS
for i = 1:N
signal = load(char(names{i}));
signal = signal.val;
ECG(:,i) = (signal(1,:)+30)./200;
RespA(:,i) = (signal(2,:)+1903)./20000;
RespC(:,i) = (signal(3,:)+1903)./20000;
RespN(:,i) = (signal(4,:)+1903)./20000;
SpO2(:,i) = signal(5,:);
end
% NORMALIZATION MIN MAX
FullECG = (ECG-min(ECG))./(max(ECG)-min(ECG));
FullRespC = (RespC-min(RespC))./(max(RespC)-min(RespC));
FullRespA = (RespA-min(RespA))./(max(RespA)-min(RespA));
FullRespN = (RespN-min(RespN))./(max(RespN)-min(RespN));
FullSpO2 = (SpO2-min(SpO2))./(max(SpO2)-min(SpO2));
% Sample Frequency
Fs = 100;
% Squared for easier peak detection
for k = 1:N
ECGSquared(:,k) = abs(FullECG((1:m),k)).^2;
end
ECGSquared = ECGSquared';
% PEAKS DETECTION
[R1,TR1]=findpeaks(ECGSquared(1,(1:m)),Fs,'MinPeakHeight',0.35);
[R2,TR2]=findpeaks(ECGSquared(2,(1:m)),Fs,'MinPeakHeight',0.35);
[R3,TR3]=findpeaks(ECGSquared(3,(1:m)),Fs,'MinPeakHeight',0.35);
[R4,TR4]=findpeaks(ECGSquared(4,(1:m)),Fs,'MinPeakHeight',0.35);
[R5,TR5]=findpeaks(ECGSquared(5,(1:m)),Fs,'MinPeakHeight',0.35);
[R6,TR6]=findpeaks(ECGSquared(6,(1:m)),Fs,'MinPeakHeight',0.35);
[R7,TR7]=findpeaks(ECGSquared(7,(1:m)),Fs,'MinPeakHeight',0.35);
%% PARA PROPOSITOS DE VISUALIZACION
t = [0:m-1]/Fs;
t1m = (1:6000)/Fs;
t30s = [0:3000-1]./Fs;
t5min = (1:30000)/Fs;
t30m = (1:180000)/Fs;
t1h = (1:360000)/Fs;
figure
plot(t,ECGSquared(1,:)),grid on, axis tight
hold on
plot(TR1,R1,'^r'),title('Peaks in ECG'),ylabel('Magnitude'), xlabel('Time (s)'),grid on, axis tight
legend('ECG', 'R')
%% INTERVALS RR
figure
title('RR intervals')
RRinterval1 = diff(TR1);
RRinterval2 = diff(TR2);
RRinterval3 = diff(TR3);
RRinterval4 = diff(TR4);
RRinterval5 = diff(TR5);
RRinterval6 = diff(TR6);
RRinterval7 = diff(TR7);
RRintMean1 = mean(RRinterval1);
RRintMean2= mean(RRinterval2);
RRintMean3= mean(RRinterval3);
RRintMean4= mean(RRinterval4);
RRintMean5= mean(RRinterval5);
RRintMean6= mean(RRinterval6);
RRintMean7= mean(RRinterval7);
disp(['1. El intervalo RR promedio es ',num2str(RRintMean1),' segundos'])
disp(['2. El intervalo RR promedio es ',num2str(RRintMean2),' segundos'])
disp(['3. El intervalo RR promedio es ',num2str(RRintMean3),' segundos'])
disp(['4. El intervalo RR promedio es ',num2str(RRintMean4),' segundos'])
disp(['6. El intervalo RR promedio es ',num2str(RRintMean5),' segundos'])
disp(['5. El intervalo RR promedio es ',num2str(RRintMean6),' segundos'])
disp(['7. El intervalo RR promedio es ',num2str(RRintMean7),' segundos'])
plot((1:length(TR1)-1)/60,RRinterval1)
xlabel('Tiempo (min)');
ylabel('RR (sec)');
h_RR = get(gca);
%% PULSACIONES
title('C�lculo de las pulsaciones por minuto')
plot((1:length(TR1)-1)/60,1./RRinterval1*60,...
(1:length(TR1)-1)/60,mean(1./RRinterval1*60),'ro');
ylabel('Pulsaciones por minuto')
%% Peaks detecion with wvs
% There are two keys for using wavelets as general feature detectors:
% The wavelet transform separates signal components into different frequency
% bands enabling a sparser representation of the signal. You can often find
% a wavelet which resembles the feature you are trying to detect. The 'sym4'
% wavelet resembles the QRS complex, which makes it a good choice for QRS
% detection. To illustrate this more clearly, extract a QRS complex and plot
% the result with a dilated and translated 'sym4' wavelet for comparison.
%[mpdict,~,~,longs] = wmpdictionary(numel(ECGa(1:m)),'lstcpt',{{'sym4',3}});
% Remove drift baseline
DetrendedECG = Detrending(ECGSquared,10);
Medians = ECGSquared-DetrendedECG;
figure
plot(Medians(1:m)),title('Drift baseline'),axis tight,grid on,xlabel('samples'),ylabel('amplitude')
%% Respiratory signal
BreathA = FullRespA(1:end,:).^2;
[RRespA,TRRespA]=findpeaks(BreathA(1:end,1),Fs,'MinPeakHeight',0.4,'MinPeakDistance',35);
findpeaks(BreathA(1:end,1),Fs,'MinPeakHeight',0.4,'MinPeakDistance',35);
RRintervalRespA = diff(TRRespA); % 203 x 1, 203 intervalos almacenados
% 10 segundos equivalen a 100 hz * 10 seg = 1000 muestras
meanRRintervalRespA = mean(RRintervalRespA);
aux = 0;
for p = 1:length(RRintervalRespA)
if RRintervalRespA(p) > 10
ApneaRespA = 1;
aux = [aux ApneaRespA];
else
ApneaRespA = 0;
aux = [aux ApneaRespA];
end
end
ApneaEventsWithRespA = aux;
plot(ApneaEventsWithRespA),title('Apnea events')
ApneaEventsA = length(ApneaEventsWithRespA(ApneaEventsWithRespA==1));
disp(['Apnea events are: ',num2str(ApneaEventsA),' in this dataset'])
% anotaciones
%[ana01er]=rdann('apnea-ecg/a01er','apn');
% figure
% plot(BreathA((1:end),1))
% hold on
% %plot(TRRespA,RRespA,'v')
% hold on
% plot(ana01er,BreathA(ana01er),'k*')
%% 3. EXTRACCION DE CARACTERISTICAS
% 3.1 Extracci�n del rms de la envolvente de la se�al y cruces por cero
mc1=[];
mc2=[];
%% EXTRACCIÓN DE CARACTERÍSTICAS
% CARACTERISTICAS:
%1. Valor RMS de la señal
%2. Coeficientes de auto-regresión
%3. Valor absoluto medio
%4. Cruces por cero
%5. Baja en la oxigenacion de la sangre SP02
%Declaramos el tama?o de la ventana
ECGm = FullECG';
%Declaramos el tama�o de la ventana
win_size = 256;
win_inc = 128; % El solapamiento de la ventana es del 50% en entrenamiento
%% MAIN FEATURES EXTRACTION
% Register
Registro = 1;
k = 1;
% We design a pace for reading the signal each minute
pace = 6000;
%for i=1:7 PARA ECG 1
for j=1:pace:996000
temp=ECGm(Registro,(j:j+6000-1));
[feature1,feature2,feature3,feature4] = extract_feature(temp',win_size,win_inc);
mc1(k) = mean(feature1);
mc2(k,:)= mean(feature2);
mc3(k) = mean(feature3);
mc4(k) = mean(feature4);
k = k + 1;
end
%end
mc5 = SPO2Detector(FullSpO2((1:966000),Registro));
mc5 = [mc5 zeros(1,length(mc1)-length(mc5))];
%% ANNOTATIONS
class_training=getclass(Registro);
feature_training = [mc1' mc2(:,1) mc2(:,2) mc2(:,3) mc2(:,4) mc3' mc4' mc5'];
[Data_training,PC_training,Ws_training,W_training,Ap_training] = entrenamiento1(feature_training,class_training);