% Script que procesa un EEG contaminado
% Clase práctica.
% 1. Cargar el fichero de datos del EEG contaminado. Muestras tomadas con una
% frecuencia de muestreo de 200 Hz.
% 2. Representar 1 segundo de dicho EEG  con el plot
% 3. Determinar "a ojo" la frecuencia del ruido añadido sobre la señal.

load eeg_cont;
Fm=200;
Tm=1/Fm;
t=0:Tm:1;
y=sn(1:length(t));
plot(t,y)
pause
% Para contar mejor abrimos una ventana ampliada de 0.2 segundos
plot(t,y)
axis([0 0.2 -0.03 0.03])
pause

% Devolvemos el axis al modo automatico (octave)
axis;

% El resultado debe ser aproximadamente sobre 75 Hz.

% Ahora queremos filtrar este ruido.
% Primero calculamos el espectro de la señal (usamos solo 1 segundo para no sobrecargar
% la memoria)


w=-2*pi*80:2*pi*80; %Vector de frecuencias donde calcularemos el espectro
% Obsérvese que tomamos unos 20 Hz más haya de donde creemos
% que acaba el espectro de la señal.

sp=fourier(y,t,w);
% Representamos la densidad espectral de energia.
plot(w/(2*pi),abs(sp)/(2*pi));
title('Espectro.');
pause

% Hacemos un zoom para ver mas claramente algunas componentes del espectro
% del EEG de un ser humano.
plot(w/(2*pi),abs(sp)/(2*pi));
axis([0 20 0 1.5e-3]);
title('Zona relevante 0-20 Hz.');
pause

% Ahora vemos la zona del ruido 

% Hacemos un zoom para ver mas claramente algunas componentes del espectro
% del EEG de un ser humano.
plot(w/(2*pi),abs(sp)/(2*pi));
axis([60 80 0 1.5e-3]);
title('Pico de ruido 75  Hz aprox.');
pause

% Ahora vamos a hacer el filtrado.
% primero buscamos los indices del vector del espectro situados en 50 y -50 Hz resp.
ind1=find(w>=-50*2*pi);
ind1=min(ind1);
ind2=find(w>=50*2*pi);
ind2=min(ind2);

% Resintetizamos la señal solo tomando en consideracion esta parte del espectro.

yr=sintetf(sp(ind1:ind2),w(ind1:ind2),t);

% Pintamos la señal "limpia".
yr=real(yr); % Para eliminar la "basura" introducida por la imprecision del 
% calculo numerico en la parte imaginaria de la integral.
subplot(211)
plot(t,y);
axis([0 1 -0.04 0.02]);
title('Con ruido.')
subplot(212)
plot(t,yr);
axis([0 1 -0.04 0.02]);
title('Sin ruido.')

