%% Detector de la presencia de onda P o no en un registro de ECG

%% Tiene pinta de que estos van a ser muy prometedores para discernir entre N y A. Para otros 
%% se podría analizar la ventana de complejo QRST.
%% Parámetros para Otros podría ser numeros de T detectado o amplitud de Complejo T, tiempo RT, varianza RT...
% Esta función se encarga de extraer acracterísticas que gusrdan relación
% con la onda P.

function [Q ampli_P mean_ampli_P var_ampli_P t_PR var_t_PR] = detector_onda_P(ecg,fs,QRS,filtro,debug)
%% Paso 0: Filtrado de la señal 
if filtro == 1
    [b,a] = butter(4,35*2/fs,'low');  % FPB de orden 4 con fc = 35 Hz
    ecg_fil = filtfilt(b,a,ecg);
end

%% Paso 1: Enventanado ondaP y latido completo
v_ant_ms = 275;
v_pos_ms = 425;
[ondaP,indP,indR,iso,latido] = enventanaP(ecg_fil,fs,QRS,v_ant_ms,v_pos_ms);
v_ant = round((v_ant_ms/10^3)*fs); % Pasa a nº de muestras importante para después

%% Paso 2: Corrección onda P y latido (llevar a 0 V la línea isoeléctrica)
% Corrección línea isoeléctrica de las ondas P
ondaP_corr = ondaP - iso*ones(1,size(ondaP,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

% Obtención del medidior de la calidad de la fiabilidad de los parámetros
Q = tot_frec_relev/sum(tabla_latidos(:,2));

% En el caso de que se detecten más de 3 patrones relevantes
if (length(patrones_relev) > 3)
    Q = 0.333;
end

% En el caso de que se detecten 5 patrones y todos con = frecuencia
if (tot_frec_relev == sum(tabla_latidos(:,2))) && (length(num_latido) == 5)
    Q = 0.2;  % Coeficiente de calidad bastante malo.
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
        ondaP_relev(i,:) = mean(ondaP_corr(pos_latido_relev,:));
    else
        ondaP_relev(i,:) = ondaP_corr(pos_latido_relev,:);
    end
    
    iso_relev(i,:) = prctile(ondaP_relev(i,:),25)*ones(1,length(ondaP_relev));
    ondaP_relev_corr(i,:) = ondaP_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 ondaP
        % Obtención de la amplitud máxima de la onda P promedio de la clase
        [ampli_P_clase(i),pos_P_clase(i)] = min(ondaP_relev_corr(i,:));
        ampli_P_clase(i) = abs(ampli_P_clase(i)); % Forzar que amplitud sea positiva
        
        % Obtener todas las ondas P pertenecientes a la misma clase
        max_P_clase = [];
        pos_P_relev = [];
        ind_max_P = [];
        for j=1:length(pos_latido_relev)
            [max_P_clase(j), ind_max_P(j)] = min(ondaP_corr(pos_latido_relev(j),:));
            pos_P_relev(j,1) = pos_R_relev(j) - v_ant + ind_max_P(j);
        end
        
        % Obtención de la variabilidad tiempo RT
        t_PR_clase(i) = mean((1/fs)*(pos_R_relev - pos_P_relev));
        var_t_PR_clase(i) = var((1/fs)*(pos_R_relev - pos_P_relev));
        
        % Obtención de la amplitud máxima promedio de las ondas P de la
        % clase misma clase, así como su varianza.
        mean_ampli_P_clase(i) = abs(mean(max_P_clase));
        var_ampli_P_clase(i) = var(max_P_clase);
       
    else                    % Si es positivo que calcule el máximo ondaP
        % Obtención de la amplitud máxima de la onda P promedio de la clase
        [ampli_P_clase(i),pos_P_clase(i)] = max(ondaP_relev_corr(i,:));
        ampli_P_clase(i) = abs(ampli_P_clase(i)); % Forzar que amplitud sea positiva
        
        % Obtener todas las ondas P pertenecientes a la misma clase
        max_P_clase = [];
        pos_P_relev = [];
        ind_max_P = [];
        for j=1:length(pos_latido_relev)
            [max_P_clase(j), ind_max_P(j)] = max(ondaP_corr(pos_latido_relev(j),:));
            pos_P_relev(j,1) = pos_R_relev(j) - v_ant + ind_max_P(j);
        end
        
        % Obtención de la variabilidad tiempo RT
        t_PR_clase(i) = mean((1/fs)*(pos_R_relev - pos_P_relev));
        var_t_PR_clase(i) = var((1/fs)*(pos_R_relev - pos_P_relev));
        
        % Obtención de la amplitud máxima promedio de las ondas P de la
        % clase misma clase, así como su varianza.
        mean_ampli_P_clase(i) = abs(mean(max_P_clase));
        var_ampli_P_clase(i) = var(max_P_clase);
    end     
    
    if debug == 1       
        figure, plot(ondaP_relev(i,:)),
        hold on, plot(ondaP_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_P = sum(w .* ampli_P_clase);

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

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

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

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

end

