%% Función que extrae caracteristicas de la onda T
% La finalidad de este función es conseguir discernir bien entre normales y
% otros. Analizando las caracteriscias de la Onda T


function [ampli_T mean_ampli_T var_ampli_T t_RT var_t_RT] = detector_onda_T(ecg,fs,picos_R,filtro,debug)
%% Paso 0: Filtrado de la señal 
if filtro == 1
    [b,a] = butter(5,40*2/fs,'low');  % FPB de orden 5 con fc = 40 Hz
    ecg_fil = filtfilt(b,a,ecg);
end

%% Paso 1: Enventanado ondaP y latido completo
v_ant_ms = 275;
v_pos_ms = 425;
[ondaT,indT,indR,iso,latido,lapso] = enventanaT(ecg,fs,picos_R,v_ant_ms,v_pos_ms);

%% Paso 2: Extracción de onda T y latido (llevar a 0 V la línea isoeléctrica)
% Corrección línea isoeléctrica de las ondas T
ondaT_corr = ondaT - iso*ones(1,size(ondaT,2));

% Corrección líena isoeléctrica latidos
latido_corr = [];
iso_latido = [];
for i=1:size(latido,1)
    iso_latido(i,:) = prctile(latido(i,:),25)*ones(1,size(latido,2));
    latido_corr(i,:) = latido(i,:) - iso_latido(i,:);
end

%% Paso 3: Obtención de parámetros (Extracción de características)
%% Paso 3.1: Analisis de los latidos correlados y que en mayor frecuencia se dan.
%% Paso 3.1.1: Clasificación de los latidos con una correlación del 0.9 en 5 clases como mucho
[tipo_latido,patron_latido,num_latido] = ClasificaSenyal(latido_corr,5,0.9); % Clasificación morfológica

if debug == 1
    % Dibujar los patrones
    tp = (1/fs)*[0:1:size(patron_latido,2)-1];
    figure, 
    if size(patron_latido,1) == 1
        plot(tp, patron_latido(1,:)), xlabel('Patron 1'),
    elseif size(patron_latido,1) == 2
        subplot(1,2,1), plot(tp, patron_latido(1,:)), xlabel('Patron 1'),
        subplot(1,2,2), plot(tp, patron_latido(2,:)), xlabel('Patron 2'),    
    elseif size(patron_latido,1) == 3
        subplot(2,2,1), plot(tp, patron_latido(1,:)), xlabel('Patron 1'),
        subplot(2,2,2), plot(tp, patron_latido(2,:)), xlabel('Patron 2'),
        subplot(2,2,3), plot(tp, patron_latido(3,:)), xlabel('Patron 3'),    
    elseif size(patron_latido,1) == 4
        subplot(2,2,1), plot(tp, patron_latido(1,:)), xlabel('Patron 1'),
        subplot(2,2,2), plot(tp, patron_latido(2,:)), xlabel('Patron 2'),
        subplot(2,2,3), plot(tp, patron_latido(3,:)), xlabel('Patron 3'),
        subplot(2,2,4), plot(tp, patron_latido(4,:)), xlabel('Patron 4'),
    elseif size(patron_latido,1) == 5
        subplot(2,3,1), plot(tp, patron_latido(1,:)), xlabel('Patron 1'),
        subplot(2,3,2), plot(tp, patron_latido(2,:)), xlabel('Patron 2'),
        subplot(2,3,3), plot(tp, patron_latido(3,:)), xlabel('Patron 3'),
        subplot(2,3,4), plot(tp, patron_latido(4,:)), xlabel('Patron 4'),
        subplot(2,3,5), plot(tp, patron_latido(5,:)), xlabel('Patron 5'),
    else
        subplot(2,3,1), plot(tp, patron_latido(1,:)), xlabel('Patron 1'),
        subplot(2,3,2), plot(tp, patron_latido(2,:)), xlabel('Patron 2'),
        subplot(2,3,3), plot(tp, patron_latido(3,:)), xlabel('Patron 3'),
        subplot(2,3,4), plot(tp, patron_latido(4,:)), xlabel('Patron 4'),
        subplot(2,3,5), plot(tp, patron_latido(5,:)), xlabel('Patron 5'),
        subplot(2,3,6), plot(tp, patron_latido(6,:)), xlabel('Patron 6'),    
    end
end

%% Paso 3.1.2: Ver cual de la frecuencia con la que se dan los patrones y el % de latidos clasificados en cada partón
% Obtención tabla de frecuencias y porcentajes de cada tipo
% Columnas: 1º tipo latido, 2º frecuencia, 3º porcentaje
tabla_latidos = tabulate(tipo_latido);

% Obtención de patrones relevantes (capaces de explicar al menos un 20%)
patrones_relev = find(tabla_latidos(:,3) >= 20);

% Obtención del numero de latidos pertenecientes a cada patrón relevante
frec_relev = tabla_latidos(patrones_relev,2);
tot_frec_relev = sum(frec_relev); % Nº total de latidos significantes

% Obtención de los pesos de ponderación de cada uno de los patrones relevantes
% Estos pesos con los que se multiplicarán por los parámetros que se han
% obtneido, con la finalidad de buscar la máxima representabilidad
for i=1:1:length(frec_relev)
    w(i) = frec_relev(i)/tot_frec_relev;
end

%% Paso 3.1.3: Análisis onda P una a una de las clases relevantes
for i=1:1:length(patrones_relev)
    %% Obtención de onda P promedio de la misma clase latido y corregir isoeléctrica
    pos_latido_relev = [];
    pos_latido_relev = find(tipo_latido == patrones_relev(i));
    
    if length(pos_latido_relev) > 1
        ondaT_relev(i,:) = mean(ondaT_corr(pos_latido_relev,:));
    else
        ondaT_relev(i,:) = ondaT_corr(pos_latido_relev,:);
    end
    
    iso_relev(i,:) = prctile(ondaT_relev(i,:),25)*ones(1,length(ondaT_relev));
    ondaT_relev_corr(i,:) = ondaT_relev(i,:) - iso_relev(i,:);
     
    %% Obtención de la amplitud media de la onda R
    pos_R_relev = [];
    pos_R_relev = indR(pos_latido_relev);
    ampli_R_mean(i) = mean(ecg(pos_R_relev)); 
    std_R(i) = std(ecg(pos_R_relev));
      
    %% Si R es - se buca mínimo y si en + se busca maximo (AMPLITUD ONDA P)
    if ampli_R_mean(i) < 0  % Si es negativo calcula el minimo ondaT
        % Obtención de la amplitud máxima de la onda T promedio de la clase
        [ampli_T_clase(i),pos_T_clase(i)] = min(ondaT_relev_corr(i,:));
        ampli_T_clase(i) = abs(ampli_T_clase(i)); % Forzar que amplitud sea positiva
        
        % Obtener todas las ondas T pertenecientes a la misma clase
        max_T_clase = [];
        pos_T_relev = [];
        ind_max_T = [];
        for j=1:length(pos_latido_relev)
            [max_T_clase(j), ind_max_T(j)] = min(ondaT_corr(pos_latido_relev(j),:)); 
            pos_T_relev(j,1) = pos_R_relev(j) + lapso +  ind_max_T(j);
        end
        
        % Obtención de la variabilidad tiempo RT
        t_RT_clase(i) = mean((1/fs)*(pos_T_relev - pos_R_relev));
        var_t_RT_clase(i) = var((1/fs)*(pos_T_relev - pos_R_relev));
        
        % Obtención de la amplitud máxima promedio de las ondas T de la
        % clase misma clase, así como su varianza.
        mean_ampli_T_clase(i) = abs(mean(max_T_clase));
        var_ampli_T_clase(i) = var(max_T_clase);
       
    else                    % Si es positivo que calcule el máximo ondaT
        % Obtención de la amplitud máxima de la onda T promedio de la clase
        [ampli_T_clase(i),pos_T_clase(i)] = max(ondaT_relev_corr(i,:));
        ampli_T_clase(i) = abs(ampli_T_clase(i)); % Forzar que amplitud sea positiva
        
        % Obtener todas las ondas T pertenecientes a la misma clase
        max_T_clase = [];
        pos_T_relev = [];
        ind_max_T = [];
        for j=1:length(pos_latido_relev)
            [max_T_clase(j), ind_max_T(j)] = max(ondaT_corr(pos_latido_relev(j),:));
            pos_T_relev(j,1) = pos_R_relev(j) + lapso +  ind_max_T(j);
        end
        
        % Obtención de la variabilidad tiempo RT
        t_RT_clase(i) = mean((1/fs)*(pos_T_relev - pos_R_relev));
        var_t_RT_clase(i) = var((1/fs)*(pos_T_relev - pos_R_relev));
        
        % Obtención de la amplitud máxima promedio de las ondas T de la
        % clase misma clase, así como su varianza.
        mean_ampli_T_clase(i) = abs(mean(max_T_clase));
        var_ampli_T_clase(i) = var(max_T_clase);
    end     
    
    if debug == 1       
        figure, plot(ondaT_relev(i,:)),
        hold on, plot(ondaT_relev_corr(i,:)),
    end
end

%% Paso 3.1.4: Sumatorio ponderado por pesos de cada una de las características
% Obtención de la amplitud de la onda P de los patrones más relevantes
ampli_T = sum(w .* ampli_T_clase);

% Obtención de la amplitud  media de la ondas P pertenecientes a las clases relevantes
mean_ampli_T = sum(w .* mean_ampli_T_clase);

% Obtención de la varianza de la amplitud de la ondas P pertenecientes a las clases relevantes
var_ampli_T = sum(w .* var_ampli_T_clase);

% Obtención de media de intervalo PR de las ondas P pertenecientes a las clases relevantes
t_RT = sum(w .* t_RT_clase);

% Obtención de media de intervalo PR de las ondas P pertenecientes a las clases relevantes
var_t_RT = sum(w .* var_t_RT_clase);

end




