% dia2_promediado.m
% Dia 2 - Por que se promedia.
% Simula el problema completo y demuestra la ley de la raiz de N.

%% Simular 40 ensayos
clear; clc;
rng(42);                    % semilla: los resultados se pueden repetir

srate  = 500;
t      = (-100:399) / srate;         % enteros / srate: t(101) es 0 exacto
ntimes = numel(t);

senal    = 8 * exp(-((t - 0.400).^2) / (2*0.080^2));
ruido_sd = 25;
ntrials  = 40;

epocas = zeros(ntrials, ntimes);     % preasignar: filas = ensayos
for k = 1:ntrials
    epocas(k,:) = senal + ruido_sd * randn(1, ntimes);
end

vent = t >= 0.300 & t <= 0.500;      % ventana de medicion
base = t < 0;                        % pre-estimulo

fprintf('%-6s %8s %8s %8s\n', 'n', 'pico', 'ruido', 'SNR');
for n = [1 2 10 40]
    p     = mean(epocas(1:n,:), 1);  % promedio de los primeros n ensayos
    amp   = mean(p(vent));           % amplitud media en la ventana
    ruido = std(p(base));            % cuanto tiembla la linea base
    fprintf('%-6d %8.2f %8.2f %8.2f\n', n, amp, ruido, amp/ruido);
end

%% La ley de la raiz de N
ruido_sd / sqrt(40)         % lo que predice la teoria

prom40 = mean(epocas, 1);   % MATLAB no deja indexar el resultado de una
std(prom40(base))           % funcion, asi que van dos pasos

% sigma_promedio = sigma_ensayo / sqrt(n)
% SNR(n)         = A / sigma_promedio = (A / sigma_ensayo) * sqrt(n)
%
% Duplicar la calidad de la senal cuesta CUADRUPLICAR los ensayos.

figure;
plot(t, epocas(1,:), 'Color', [.75 .75 .75]); hold on
plot(t, mean(epocas(1:10,:), 1), 'LineWidth', 1.2)
plot(t, prom40, 'LineWidth', 1.8)
plot(t, senal, 'k--')
legend('1 ensayo', 'promedio de 10', 'promedio de 40', 'componente latente')
xlabel('Tiempo (s)'), ylabel('Amplitud (\muV)')

%% Condicional y funcion anonima
n_malos = sum(max(abs(epocas), [], 2) > 100);   % ensayos que pasan de +-100 uV
if n_malos / ntrials > 0.25
    warning('Se rechazaria mas del 25%% de los ensayos');
else
    fprintf('Rechazo: %.1f%%\n', 100 * n_malos / ntrials);
end

razon_sn = @(w) mean(w(vent)) / std(w(base));   % una formula en una variable
razon_sn(mean(epocas(1:10,:), 1))

% No le pongas a una variable el nombre de una funcion: snr, log, mean, max e i
% ya existen en MATLAB. Una variable con ese nombre las tapa en silencio.
% "which log" muestra cual esta ganando; "clear log" arregla el desastre.

%% La funcion promediar_n vive en su propio archivo
p10 = promediar_n(epocas, 10);
fprintf('SNR del promedio de 10: %.2f\n', razon_sn(p10));