%% Función que realiza el filtrado óptimo para analizar la onda P
% Lo que realiza es un filtrado de tipo butterworth paso banda con las
% frecuencias de corte que se le pasen y el orden que se le pase.
%
% Parámetros de entrada:
%   * ecg: registro de ecg de entrada a filtrar
%   * fs: frecuencia de muestreo del ecg
%   * fci: frecuencia de corte inferior del filtro paso banda (Hz)
%   * fcs: frecuencia de corte superior del filtro paso banda (Hz)
%   * orden [opcional]: Orden del filtro 
%   * res_plot [opcional]: Si vale 1 se visualiza la respuesta frecuencial 
%                          del filtro diseñado. Si vale cero 0 no se 
%                          representa la respuesta.
%
% Parámetros de salida: 
%   * ecg_p: registro de ecg de salida al filtrar con línea 
%                 isoeléctrica sin correguir.
%   * ecg_p_corr: línea isoeléctrica corregida mediante la medianda. De 
%                 este modo se va a poder operar correctamente para
%                 cálculos de amplitud de onda P y extracción de
%                 características morfológicas.


function [ecg_p, ecg_p_corr] = filtrado_Paso_Banda(ecg,fs,fci,fcs,orden,res_plot)
%% Configuración de parámetros opcionales
if nargin < 5
    orden = 2;      % Orden por defecto (2º orden)
    res_plot = 0;   % No se representa la respuesta
elseif nargin < 6
    orden = orden;  % Orden introducido por el usuario
    res_plot = 0;   % No se representa la respuesta       
else 
    orden = orden;  % Orden introducido por el usuario
    res_plot = res_plot; % Mostrar la respuesta del filtro según el usuario      
end

%% Filtrado Paso-Banda con respuesta Butterworth
Wi = fci/(fs/2);  % Freciencia de corte inferior normalizada
Ws = fcs/(fs/2);  % Freciencia de corte superior normalizada
Wn = [Wi Ws];     % Frecuencias paso banda normalizadas
[B,A] = butter(orden,Wn); 
ecg_p = filtfilt(B, A, ecg);

%% Correción de la línea isoeléctrica del ecg filtrado (con mediana de señal sin filtrar)
ecg_p_corr = ecg_p - mode(ecg)*ones(1,length(ecg_p))';

%% Representación de la respuesta del filtro diseñado
if res_plot == 1
    figure, freqz(B,A,1e6,fs);   % Visualizar la respuesta en f del filtro
    figure, hold on,             % Visualizar los registros (sin fil, con fil y filt corregido
    plot(ecg), plot(ecg_p,'r'), plot(ecg_p_corr,'g'), plot(zeros(1,length(ecg)),'k')
end
end
