% dia4_filtros.m
% Dia 4 - Filtrar sin romper.
% Usa pasoalto1.m, que vive en su propio archivo en esta misma carpeta.

%% Ocho numeros: una media movil hecha a mano
clear; clc;

x = [2 8 3 9 4 10 5 11];

y = zeros(1, numel(x));
for k = 2:numel(x)-1
    y(k) = (x(k-1) + x(k) + x(k+1)) / 3;
end
y
% y =  0   4.33   6.67   5.33   7.67   6.33   8.67   0
% Los saltos rapidos se aplanaron, la tendencia lenta quedo: eso ya es un
% filtro paso bajo.

%% Lo mismo, con conv
% y[k] = b1*x[k-1] + b2*x[k] + b3*x[k+1]  =  sum_j b_j * x[k-j]
conv(x, [1 1 1]/3, 'same')
conv(x, [0.25 0.5 0.25], 'same')     % el centro pesa el doble

% En el medio son identicas al bucle. En los bordes el bucle dejo ceros y conv
% devolvio 3.33 y 5.33: ninguno tiene razon, porque ahi no hay datos. Todo
% filtro tiene ese problema en los extremos, y en un registro real los extremos
% son el principio del archivo y cada corte entre bloques.

%% Que frecuencias sobreviven
srate = 500;  L = 9;
H = abs(fft(ones(1,L)/L, 512));
f = (0:511) * srate / 512;

for fq = [2 30 50]
    [~, idx] = min(abs(f - fq));
    fprintf('ganancia en %5.1f Hz: %.3f\n', f(idx), H(idx));
end
% Una media movil de 9 puntos a 500 Hz deja pasar casi intacto lo de 2 Hz,
% se come casi la mitad de lo de 30 Hz y deja pasar un 11% de los 50 Hz.

figure;
plot(f(1:256), H(1:256))
xlabel('Frecuencia (Hz)'), ylabel('Ganancia')
title('Respuesta en frecuencia de una media movil de 9 puntos')

%% El filtro que inventa un componente
srate = 500;
t = (-100:500) / srate;            % enteros / srate: sin residuos de punto flotante
onda = 8 * exp(-((t - 0.400).^2) / (2*0.080^2));   % puramente POSITIVA

for fc = [0.1 0.5 2 5]
    y = pasoalto1(onda, fc, srate);
    fprintf('fc = %4.1f Hz  ->  pico %5.2f uV,  posterior %6.2f uV\n', ...
        fc, max(y), min(y(t > 0.6)));
end
fprintf('original      ->  pico %5.2f uV,  posterior %6.2f uV\n', ...
    max(onda), min(onda(t > 0.6)));

figure; hold on
plot(t, onda, 'k', 'LineWidth', 1.5)
plot(t, pasoalto1(onda, 0.5, srate))
plot(t, pasoalto1(onda, 2.0, srate))
plot([t(1) t(end)], [0 0], 'Color', [.6 .6 .6])
legend('sin filtrar', 'paso alto 0,5 Hz', 'paso alto 2 Hz')
xlabel('Tiempo (s)'), ylabel('Amplitud (\muV)')
title('El valle posterior no esta en los datos: lo fabrico el filtro')

%% Como se hace en la practica (necesita EEGLAB + ERPLAB y un dataset cargado)
% Paso alto sobre el EEG CONTINUO, antes de epocar:
%
%   EEG = pop_basicfilter(EEG, 1:33, 'Filter', 'highpass', 'Design', 'butter', ...
%       'Cutoff', 0.1, 'Order', 2, 'RemoveDC', 'on', 'Boundary', 'boundary');
%
% Paso bajo mucho despues, sobre el ERP ya promediado:
%
%   ERP = pop_filterp(ERP, 1:33, 'Filter', 'lowpass', 'Design', 'butter', ...
%       'Cutoff', 20, 'Order', 2);