MATLAB para potenciales evocados

• MATLAB con la Signal Processing Toolbox

  La mayoría de las universidades tienen licencia de campus; consúltala con tu unidad de informática. Si no alcanzas a activarla, MATLAB Online sirve para todo lo que hay aquí. ERPLAB exige la Signal Processing Toolbox. Las de Statistics and Machine Learning y Parallel Computing sólo hacen falta si más adelante haces decoding (MVPA). Los desarrolladores de ERPLAB desaconsejan MATLAB 2025a en adelante; si puedes elegir versión, elige una anterior.
• EEGLAB desde el sitio del SCCN

  Descárgalo de eeglab.org, no como ZIP de GitHub: a ese le faltan plugins. Descomprímelo en una carpeta sin espacios ni tildes en la ruta — es la causa número uno de errores inexplicables en máquinas hispanohablantes.
• ERPLAB 13.00 como extensión de EEGLAB

  Lo más limpio es abrir EEGLAB y usar File → Manage EEGLAB extensions. La alternativa es descomprimir el release dentro de eeglab/plugins/ y reiniciar. ERPLAB 13 trae dos caras: ERPLAB Studio (interfaz nueva) y ERPLAB Classic (el menú clásico dentro de EEGLAB). El curso usa Classic, porque es el que se traduce a script.
• Los datos: ERP CORE, paradigma N400

  40 participantes, siete componentes, con los scripts de procesamiento publicados (erpinfo.org/erp-core, CC BY-SA 4.0). Baja sólo el experimento N400 — son varios GB si bajas todo. Es el mismo set con el que trabaja el libro de Luck, así que cualquier duda tiene respuesta escrita en alguna parte.
• Verificar que todo se ve desde MATLAB

  Abre MATLAB, escribe eeglab y presiona Enter. Debería abrirse la ventana principal y aparecer un menú ERPLAB. Después corre las tres líneas de abajo: si alguna devuelve vacío, algo no quedó en el path.

which eeglab                 % ¿ve MATLAB a EEGLAB?
which pop_basicfilter        % ¿ve a ERPLAB?
ver('signal')                % ¿está la Signal Processing Toolbox?

Empieza escribiendo

Abre MATLAB. La ventana grande del centro es la Command Window: lo que escribes ahí se ejecuta al presionar Enter. Prueba esto, línea por línea:

8 + 5
v = [2 8 3 9]
mean(v)
v(2)
v(2:3)

Tres cosas ya quedaron claras sin que nadie las anunciara. Cuando no asignas el resultado a nada, MATLAB lo guarda en una variable llamada ans. Los corchetes [ ] arman un vector: cuatro números en fila. Y v(2) devolvió 8, el segundo elemento — en MATLAB los índices empiezan en 1, no en 0. Si vienes de Python o de R, esta es la primera y más frecuente fuente de errores por un puesto.

El punto y coma al final de una línea silencia la salida. Con vectores de 500 mil números lo vas a agradecer:

ruido = randn(1, 500000);    % sin salida en pantalla
size(ruido)                  % 1 fila, 500000 columnas

De índices a segundos

Un registro de EEG es una lista de voltajes tomados a intervalos regulares. Si el equipo muestrea a 500 Hz, hay una muestra cada 2 ms. El vector de tiempo se construye con el operador dos puntos, que significa «desde : paso : hasta»:

srate = 500;                      % frecuencia de muestreo, en Hz
t = (0:499) / srate;    % un segundo de tiempo, 500 muestras
numel(t)                          % 500 muestras
t(end)                            % 0.998 s — la última

x = sin(2*pi*10*t);               % un seno de 10 Hz
plot(t, x)
xlabel('Tiempo (s)'), ylabel('Amplitud (\muV)')

El truco que vas a usar todos los días: en vez de contar índices a mano, deja que MATLAB los busque por vos con una comparación. t >= 0.2 devuelve un vector de verdaderos y falsos, y ese vector sirve como índice.

vent = t >= 0.2 & t < 0.4;   % una ventana de 200 ms
sum(vent)                     % 100 muestras caen dentro
mean(x(vent))                 % amplitud media en esa ventana

Eso es, literalmente, lo que hace ERPLAB cuando le pides «amplitud media entre 300 y 500 ms». No hay magia debajo.

Dos dimensiones: canales por tiempo

Un electrodo da un vector. Treinta electrodos dan una matriz: filas para los canales, columnas para el tiempo. Esa es exactamente la forma de EEG.data en EEGLAB, y saber leerla es la mitad del curso.

datos = zeros(4, numel(t));       % 4 canales × 500 muestras
for ch = 1:4
    datos(ch,:) = ch * sin(2*pi*10*t);
end

size(datos)          % 4   500
size(datos, 1)       % 4   → número de canales
datos(3, :)          % todo el canal 3
datos(:, 100)        % los 4 canales en la muestra 100
mean(datos, 2)       % promedio en el tiempo, un valor por canal

El segundo argumento de mean es la dimensión sobre la que promedias: 1 recorre las filas (promedia canales entre sí), 2 recorre las columnas (promedia el tiempo dentro de cada canal). Equivocarse de dimensión produce resultados que parecen razonables, así que compruébalo siempre con size.

El punto que cambia todo

MATLAB nació como lenguaje de álgebra lineal, así que * significa producto matricial. Para multiplicar elemento por elemento — que es lo que casi siempre quieres con señales — se antepone un punto:

a = [1 2 3];
b = [4 5 6];
a .* b        % 4  10  18   ← elemento a elemento
a * b'        % 32          ← producto punto (b' es b transpuesto)
a * b         % ERROR: Inner matrix dimensions must agree

Lo mismo vale para ./ y .^. Cuando veas el error Matrix dimensions must agree, el 90% de las veces falta un punto.

Corrección de línea base, a mano

Restar la línea base es sólo eso: calcular el promedio del período previo al estímulo y restarlo. Escríbelo una vez a mano y nunca más vas a tratar la opción de ERPLAB como una caja negra.

t   = (-100:399) / srate;   % época típica de ERP; enteros / srate: t(101) es 0 exacto
sig = 8 * exp(-((t - 0.400).^2) / (2*0.080^2));   % un pico en 400 ms

base = t < 0;                     % el período pre-estímulo
sig_corr = sig - mean(sig(base)); % restar su media

mean(sig(base))        % antes: distinto de cero
mean(sig_corr(base))   % después: cero (o 1e-16, que es cero)

Del Command Window al script

Todo lo anterior se escribió línea por línea y se perdió al cerrar. Un análisis no se hace así: se escribe en un archivo que se puede volver a correr, corregir y mandar por correo. Ese archivo es un script, un texto con extensión .m.

1. Pestaña Home, botón New → Script.
2. Guárdalo como dia1.m en tu carpeta del curso — sin tildes ni espacios en la ruta, igual que con EEGLAB.
3. Escribe, guarda, y presiona Run (o F5): MATLAB ejecuta el archivo completo, de arriba abajo, como si hubieras tecleado cada línea.

% dia1.m — vectores, tiempo y ventanas de medición
clear                        % borra las variables que quedaron de antes
clc                          % limpia la pantalla

srate = 500;
t = (0:499) / srate;
x = sin(2*pi*10*t);
plot(t, x)

Ese clear de la segunda línea evita el error más difícil de ver de todos: una variable vieja con el mismo nombre sobrevive de la corrida anterior, el script la usa sin quejarse, y los resultados salen mal sin ningún mensaje. Si un script sólo funciona cuando lo corres dos veces, ahí está el problema.

Un atajo que vas a usar toda la semana: dos signos de porcentaje %% parten el archivo en secciones. Con el cursor dentro de una, Ctrl+Enter (Cmd+Enter en Mac) ejecuta sólo esa sección. Así pruebas el bloque que estás escribiendo sin volver a correr los cinco anteriores.

clear; clc;
srate = 250;
t     = (0:499) / srate;   % dos segundos
alfa  = 20 * sin(2*pi*10*t);         % 10 Hz, 20 µV de amplitud

numel(t)
std(alfa)

red   = 3 * sin(2*pi*____*t);    % 50 Hz, 3 µV
senal = alfa + ____;

plot(t, senal)
xlabel('Tiempo (s)'), ylabel('Amplitud (\muV)')

vent = t >= 0.5 & t < ____;      % de 0,5 a 1,5 segundos
sum(vent)                        % ¿cuántas muestras deberían caer adentro?

clear; clc;
srate = 250;
t     = (0:499) / srate;

alfa  = 20 * sin(2*pi*10*t);
red   = 3  * sin(2*pi*50*t);
senal = alfa + red;

plot(t, senal)
xlabel('Tiempo (s)'), ylabel('Amplitud (\muV)')

vent = t >= 0.5 & t < 1.5;

fprintf('muestras en la ventana: %d\n', sum(vent));
fprintf('sd completa: %.4f\n', std(senal));
fprintf('sd ventana : %.4f\n', std(senal(vent)));

% Son casi iguales porque esta señal no cambia de carácter con el tiempo:
% cada segundo se parece a cualquier otro. Un ERP no es así — el componente
% ocurre en un momento y no en los demás — y por eso ahí la ventana sí importa.

mean(datos)
mean(datos, 2)
mean(datos(:))

Primero los números

Simulemos el problema entero. Un componente positivo de 8 µV en 400 ms, ruido con desviación estándar de 25 µV — proporciones realistas para un P3 de un solo ensayo — y 40 ensayos.

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 medición
base = t < 0;                        % pre-estímulo

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));            % cuánto tiembla la línea base
    fprintf('%-6d %8.2f %8.2f %8.2f\n', n, amp, ruido, amp/ruido);
end

Mira la tercera columna. El ruido del período basal cae de 28 a 3,8 al pasar de 1 a 40 ensayos: no dividido por 40, sino por algo cercano a la raíz de 40, que es 6,3. Compruébalo:

ruido_sd / sqrt(40)         % 3.95 — lo que predice la teoría

prom40 = mean(epocas, 1);   % MATLAB no deja indexar el resultado de una
std(prom40(base))           % función, así que van dos pasos:  3.83

Los mismos números, ahora con letras

En la tabla hay tres cantidades y conviene ponerles nombre antes de escribir nada: el ruido de un ensayo suelto, que aquí vale 25 µV; la cantidad de ensayos que entraron al promedio; y el ruido que queda después de promediar, que fue 3,83.

Ahora agrega la amplitud del componente, que el promedio no cambia porque es idéntica en todos los ensayos: 6,30 µV medidos entre 300 y 500 ms. Divide una cosa por la otra y tienes la razón señal-ruido del promedio:

La forma de la fórmula es todo lo que hay que recordar: la n entra bajo una raíz, y entra multiplicando. Nada de esto es una propiedad del EEG —es aritmética de promedios— pero de ahí sale la consecuencia práctica que gobierna cualquier diseño de experimento:

duplicar la calidad de la señal cuesta cuadruplicar los ensayos. Pasar de 40 a 80 ensayos mejora un 41%; pasar de 40 a 160 mejora al doble. Es la razón por la que los paradigmas de ERP son tan largos y tan aburridos, y la razón por la que una condición con pocos ensayos válidos no es comparable con una que los tiene todos.

Las tres estructuras que necesitas

Ya usaste un for dos veces. Estas son las otras dos piezas.

% Condicional
n_malos = sum(max(abs(epocas), [], 2) > 100);   % ensayos que pasan de ±100 µV
if n_malos / ntrials > 0.25
    warning('Se rechazaría más del 25%% de los ensayos');
else
    fprintf('Rechazo: %.1f%%\n', 100 * n_malos / ntrials);
end

% Función anónima: una fórmula guardada en una variable
razon_sn = @(w) mean(w(vent)) / std(w(base));
razon_sn(mean(epocas(1:10,:), 1))

Cuando la operación no cabe en una línea, se escribe una función de verdad. En MATLAB una función puede ir en su propio archivo .m — con el mismo nombre que la función — o al final de un script, nunca al medio ni al principio. Esa restricción no existe en Octave, y es un error que sólo aparece cuando entregas el script a otra persona.

function p = promediar_n(epocas, n)
% PROMEDIAR_N  Promedia los primeros n ensayos de una matriz de épocas.
%   epocas : matriz ensayos × tiempo
%   n      : cuántos ensayos incluir
    if n > size(epocas, 1)
        error('Pediste %d ensayos y sólo hay %d.', n, size(epocas,1));
    end
    p = mean(epocas(1:n, :), 1);
end

Las líneas de comentario que van justo debajo de function son lo que aparece cuando alguien escribe help promediar_n. Escríbelas siempre: es la única documentación que vas a leer en seis meses.

clear; clc; rng(42);

srate = 500;
t     = (-100:399) / srate;         % enteros / srate: t(101) es 0 exacto
senal = 8 * exp(-((t - 0.400).^2) / (2*0.080^2));
vent  = t >= 0.300 & t <= 0.500;
base  = t < 0;
nmax  = 200;

epocas = senal + 25 * randn(nmax, numel(t));   % 200 ensayos de una vez

curva = zeros(1, nmax);
for n = 1:nmax
    p = mean(epocas(1:n,:), 1);
    curva(n) = mean(p(vent)) / std(p(base));
end

plot(1:nmax, curva)
xlabel('Ensayos promediados'), ylabel('SNR')

k = find(curva > ____, 1);      % el 1 significa «sólo el primero»

if isempty(k)
    fprintf('No llega a SNR 2 con %d ensayos\n', ____);
else
    fprintf('Primer n con SNR > 2: %d\n', k);
end

clear; clc; rng(42);

srate = 500;
t     = (-100:399) / srate;         % enteros / srate: t(101) es 0 exacto
senal = 8 * exp(-((t - 0.400).^2) / (2*0.080^2));
vent  = t >= 0.300 & t <= 0.500;
base  = t < 0;
nmax  = 200;

figure; hold on
for sd = [15 25 40]
    epocas = senal + sd * randn(nmax, numel(t));

    curva = zeros(1, nmax);
    for n = 1:nmax
        p = mean(epocas(1:n,:), 1);
        curva(n) = mean(p(vent)) / std(p(base));
    end

    plot(1:nmax, curva, 'LineWidth', 1.2)

    k = find(curva > 2, 1);
    if isempty(k)
        fprintf('sd = %2d uV  ->  no llega a SNR 2 en %d ensayos\n', sd, nmax);
    else
        fprintf('sd = %2d uV  ->  primer n con SNR > 2: %d\n', sd, k);
    end
end

plot([1 nmax], [2 2], '--')      % la meta
legend('15 \muV', '25 \muV', '40 \muV', 'SNR = 2')
xlabel('Ensayos promediados'), ylabel('SNR')

Un dato con nombres

Hasta aquí los datos eran números sueltos. Un registro real trae además la frecuencia de muestreo, los nombres de los electrodos, los eventos, la historia de lo que se le hizo. Todo eso viaja junto en una estructura: una variable con campos nombrados, a los que se llega con un punto.

Vamos a fabricar uno de juguete, con la misma forma que el de verdad. Todo lo que sigue en esta sesión funciona sobre él, así que puedes practicar aunque los datos todavía se estén descargando:

clear; clc;

EEG = struct();
EEG.setname  = 'juguete';
EEG.srate    = 500;
EEG.data     = randn(4, 5000);      % 4 canales, 10 segundos de ruido
EEG.pnts     = size(EEG.data, 2);
EEG.chanlocs = struct('labels', {'Fz','Cz','Pz','Oz'});
EEG.event    = struct('type',    {211, 201, 221, 201, 211}, ...
                      'latency', {500, 900, 1800, 2200, 3500});

fprintf('%s: %d canales, %.1f s a %g Hz\n', ...
    EEG.setname, size(EEG.data,1), EEG.pnts/EEG.srate, EEG.srate);

numel(EEG.chanlocs)     % cuántos electrodos
numel(EEG.event)        % cuántos eventos
isfield(EEG, 'EVENTLIST')   % 0: este dataset todavía no pasó por ERPLAB

Los tres puntos ... continúan una instrucción en la línea siguiente. Los vas a usar mucho: las llamadas a ERPLAB son largas.

Y ojo con lo que hizo esa llave en struct('labels', {'Fz','Cz','Pz','Oz'}). No creó un campo con cuatro nombres adentro: creó cuatro estructuras, una por electrodo, cada una con su propio labels. Ésa es exactamente la forma que tienen chanlocs y event en un dataset real, y es la razón de la sintaxis rara que viene en un momento.

Los campos que importan

Cuando cargas un .set, EEGLAB deja en el workspace una estructura llamada EEG con decenas de campos. Estos son los que vas a tocar:

Recorrer un struct array

Los campos que tienen un elemento por canal o por evento son arreglos de estructuras. Para sacarles una columna completa hay dos sintaxis, y conviene aprenderlas juntas: corchetes cuando el contenido es numérico, llaves cuando es texto.

etiquetas = {EEG.chanlocs.labels};        % celda con todos los nombres
iPz = find(strcmp(etiquetas, 'Pz'));      % ¿en qué posición está Pz?
fprintf('Pz es el canal %d\n', iPz);

codigos = [EEG.event.type];               % vector con todos los códigos
for c = unique(codigos)
    fprintf('  código %3d: %4d veces\n', c, sum(codigos == c));
end

Ese pequeño bucle es lo primero que hay que correr con datos nuevos. Si el paradigma dice que hubo 120 targets y el conteo dice 118, quieres saberlo antes de promediar, no después.

Tres unidades para la misma cosa

Una latencia se puede expresar en muestras, en segundos o en milisegundos, y el objeto EEG usa las tres a la vez. Practica la conversión ahora, con números chicos y a mano, porque equivocarse aquí no produce ningún error: produce una ventana de medición corrida que igual entrega un resultado.

EEG.event(3).latency                        % está en MUESTRAS
EEG.event(3).latency / EEG.srate            % → segundos
1000 * EEG.event(3).latency / EEG.srate     % → milisegundos

% Al revés: dentro de una época que empieza en -200 ms,
% ¿qué muestra corresponde a 300 ms?
round((300 - (-200)) / 1000 * EEG.srate) + 1

El + 1 del final es el mismo de siempre: la primera muestra de la época es la número 1, no la 0. Y el round no es cosmético: sin él queda un índice con decimales y MATLAB contesta Subscript indices must be positive integers.

En la práctica casi nunca vas a hacer esta cuenta. Sobre un dataset cargado de verdad, EEG.times ya trae el eje en milisegundos y una comparación lógica —la misma del día 1— resuelve el problema sin aritmética:

vent = EEG.times >= 300 & EEG.times <= 500;   % en ms, no en segundos

Si escribes 0.3 en vez de 300, vent queda vacío, el promedio da NaN y no aparece ningún mensaje de error. Ése es el precio de que las tres unidades convivan.

Rutas, celdas y nombres de archivo

Los sujetos viven en carpetas. Armar rutas pegando texto con [ ] funciona pero se rompe entre sistemas operativos; fullfile pone el separador correcto solo.

raiz = '~/Documents/ERP_CORE/N400';
sujetos = {'sub-001','sub-002','sub-010'};

for k = 1:numel(sujetos)
    s = sujetos{k};                                  % llaves: sacar de una celda
    archivo = fullfile(raiz, s, [s '_N400.set']);
    fprintf('%s\n', archivo);
end

sprintf('sub-%03d', 7)     % 'sub-007' — ceros a la izquierda

La GUI escribe tu script

Este es el atajo más útil de EEGLAB y casi nadie lo menciona el primer día. Haz una operación con los menús — cargar un archivo, filtrar, lo que sea — y después escribe:

eegh

EEGLAB imprime en la consola el comando exacto, con todos sus argumentos, que equivale a lo que acabas de hacer con el mouse. Copias, pegas en el editor, y ya tienes media línea de pipeline escrita sin memorizar nada. El flujo de trabajo de todo el curso es ese: hacerlo una vez apuntando, pedirle a eegh la traducción, y quedarte con el script.

[ALLEEG, EEG, CURRENTSET, ALLCOM] = eeglab;   % arranca y limpia el workspace

EEG = pop_loadset('filename', 'sub-001_N400.set', ...
                  'filepath', '/ruta/a/N400/sub-001/');
[ALLEEG, EEG, CURRENTSET] = eeg_store(ALLEEG, EEG);
eeglab redraw                                  % refresca la ventana de EEGLAB

ALLEEG es la lista de todos los datasets abiertos y CURRENTSET el índice del activo. Trabajando por script casi siempre basta con EEG; los otros dos importan cuando quieres ver el resultado en la interfaz.

fprintf('\n--- %s ---\n', EEG.setname);
fprintf('canales   : %d\n', size(EEG.data, 1));
fprintf('muestreo  : %g Hz\n', EEG.srate);
fprintf('duración  : %.1f minutos\n', EEG.pnts / EEG.srate / 60);

tipos = {EEG.event.type};              % llaves: funciona en los dos casos

if isnumeric(EEG.event(1).type)
    tipos = cellfun(@num2str, tipos, 'UniformOutput', false);
end

unicos = ____(tipos);                  % los códigos distintos, sin repetir

for k = 1:numel(unicos)
    n = sum(____(tipos, unicos{k}));   % ¿cuántas veces aparece este código?
    fprintf('   %-6s %4d veces\n', unicos{k}, n);
end

function inventario(EEG)
% INVENTARIO  Resumen de un dataset de EEGLAB, antes de tocar nada.
%   inventario(EEG) imprime forma, duración, estado y conteo de eventos.

    fprintf('\n--- %s ---\n', EEG.setname);
    fprintf('canales   : %d\n', size(EEG.data, 1));
    fprintf('muestreo  : %g Hz\n', EEG.srate);

    if isfield(EEG, 'trials') && EEG.trials > 1
        fprintf('estado    : epocado, %d épocas de %d muestras\n', ...
            EEG.trials, EEG.pnts);
    else
        fprintf('estado    : continuo, %.1f minutos\n', EEG.pnts / EEG.srate / 60);
    end

    if ~isfield(EEG, 'event') || isempty(EEG.event)
        fprintf('sin eventos\n');
        return
    end

    tipos = {EEG.event.type};
    if isnumeric(EEG.event(1).type)
        tipos = cellfun(@num2str, tipos, 'UniformOutput', false);
    end

    unicos = unique(tipos);
    fprintf('eventos   : %d en total, %d códigos distintos\n', ...
        numel(tipos), numel(unicos));
    for k = 1:numel(unicos)
        fprintf('   %-6s %4d veces\n', unicos{k}, sum(strcmp(tipos, unicos{k})));
    end
end

EEG.event(1).type = 'S211';
EEG.event(2).type = 'S201';

codigos = [EEG.event.type];
sum(codigos == 211)

Ocho números

Toma esta lista y reemplaza cada valor por el promedio de sí mismo y sus dos vecinos:

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

Los saltos de la serie original —de 2 a 8, de 3 a 9— se aplanaron: la variación rápida se fue, la tendencia lenta quedó. Eso ya es un filtro paso bajo. Todo lo demás es esa misma operación hecha con pesos mejor elegidos.

Ahora los nombres. Los tres pesos [1/3 1/3 1/3] son el núcleo del filtro. Deslizarlos por la señal, multiplicando y sumando en cada posición, se llama convolución. Y el bucle que acabas de escribir, puesto en símbolos, es esto:

Con más pesos la lista se hace larga, así que se escribe con una suma:

Eso es un filtro digital, entero. Butterworth, FIR, paso alto, paso bajo: todos son esta misma suma, y lo único que los distingue es cuántos pesos hay y cuánto valen. MATLAB la calcula en una línea:

conv(x, [1 1 1]/3, 'same')
conv(x, [0.25 0.5 0.25], 'same')        % el centro pesa el doble

Compara con la salida del bucle: en el medio son idénticas, pero en los extremos el bucle dejó ceros y conv devolvió 3,33 y 5,33. Ninguno de los dos tiene razón, porque ahí no hay datos: el bucle se negó a inventarlos y conv supuso que la señal seguía con ceros.

Todo filtro tiene ese problema en los bordes, y en un registro real los bordes son los primeros milisegundos del archivo, los últimos, y cada corte entre bloques. Guárdalo: reaparece más abajo como 'Boundary','boundary'.

Qué frecuencias sobreviven

Un filtro se describe por cuánto deja pasar de cada frecuencia. Esa curva se obtiene transformando el núcleo, y es todo lo que hay detrás de la palabra «respuesta en frecuencia»:

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

for fq = [2 30 50]
    [~, i] = min(abs(f - fq));
    fprintf('ganancia en %5.1f Hz: %.3f\n', f(i), H(i));
end

Una media móvil de 9 puntos a 500 Hz deja pasar casi intacto lo que ocurre a 2 Hz, se come casi la mitad de lo que ocurre a 30 Hz y deja pasar apenas un 11% de los 50 Hz de la red eléctrica. Y no tiene ningún parámetro salvo su largo. Los filtros «de verdad» (Butterworth, FIR) sólo agregan control sobre la frecuencia de corte —dónde empieza a atenuar— y sobre el orden —cuán abrupta es la caída, que se mide en dB por octava.

El filtro que inventa un componente

Aquí está el motivo por el que este día existe. Un paso alto de primer orden cabe en cinco líneas; escribámoslo y pasémosle una onda que es puramente positiva:

srate = 500;
t = (-100:500) / srate;            % enteros / srate: sin residuos de punto flotante
onda = 8 * exp(-((t - 0.400).^2) / (2*0.080^2));

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)));

function y = pasoalto1(x, fc, srate)
    a = 1 / (1 + 2*pi*fc/srate);
    y = zeros(size(x));
    for n = 2:numel(x)
        y(n) = a * (y(n-1) + x(n) - x(n-1));
    end
end

La señal de entrada no baja de cero en ningún momento. La señal filtrada a 2 Hz tiene un valle de −2,13 µV después del pico, y el pico mismo perdió más de la mitad de su amplitud. Grafica onda y y superpuestas y míralo: ese valle no está en los datos, lo fabricó el filtro. Si tu hipótesis es sobre un componente negativo tardío, un paso alto agresivo te lo regala.

Por eso la recomendación de ERPLAB es conservadora: paso alto en 0,1 Hz o menos, paso bajo en 20 Hz o más, con pendiente de 12 dB/octava, salvo que sepas exactamente lo que estás haciendo. Y por eso, en un paper de ERP, los parámetros del filtro son parte de los métodos, no un detalle técnico.

Cómo se hace en la práctica

Dos reglas de orden que vienen de la misma física: el paso alto necesita ver tramos largos para estimar la deriva, así que va sobre el EEG continuo, antes de epocar; el paso bajo sólo suaviza, así que puede ir al final, sobre el ERP ya promediado, donde no afecta la detección de artefactos.

% Paso alto sobre el EEG continuo, canales 1 a 33
EEG = pop_basicfilter(EEG, 1:33, 'Filter', 'highpass', 'Design', 'butter', ...
    'Cutoff', 0.1, 'Order', 2, 'RemoveDC', 'on', 'Boundary', 'boundary');

% Paso bajo, mucho después, sobre el ERP promediado
ERP = pop_filterp(ERP, 1:33, 'Filter', 'lowpass', 'Design', 'butter', ...
    'Cutoff', 20, 'Order', 2);

El argumento 'Boundary','boundary' le dice al filtro que no cruce los cortes del registro: si hubo una pausa entre bloques, filtrar a través de ella produce un transitorio enorme que después se ve como artefacto.

Referencia: aritmética entre canales

El voltaje sólo existe como diferencia entre dos puntos, así que «rereferenciar» es simplemente restar. ERPLAB lo hace con fórmulas escritas en texto plano, una por canal:

% Referencia al promedio de mastoides, y un canal bipolar de EOG
EEG = pop_eegchanoperator(EEG, { ...
    'ch1 = ch1 - (ch30 + ch31)/2 label Fp1', ...
    'ch2 = ch2 - (ch30 + ch31)/2 label Fp2', ...
    'nch34 = ch32 - ch33 label VEOG-bipolar' }, ...
    'ErrorMsg', 'popup', 'KeepChLoc', 'on', 'Warning', 'on');

ch reemplaza un canal existente y nch crea uno nuevo. Para un montaje de 32 canales no vas a escribir 32 líneas a mano: el Reference Assistant de la interfaz las genera y las guarda en un .txt que después le pasas a la función como ruta.

clear; clc;
srate = 500;
t     = (-100:399) / srate;         % enteros / srate: t(101) es 0 exacto
onda  = 8 * exp(-((t - 0.400).^2) / (2*0.080^2));
vent  = t >= 0.300 & t <= 0.500;

y01 = pasoalto1(onda, 0.1, srate);

a0 = mean(onda(vent));
a1 = mean(y01(vent));
fprintf('sin filtrar   : %.3f uV\n', a0);
fprintf('paso alto 0,1 : %.3f uV  (%.1f%% menos)\n', a1, 100*(1 - a1/a0));

y10 = pasoalto1(onda, ____, srate);

a2 = mean(y10(____));
fprintf('paso alto 1,0 : %.3f uV  (%.1f%% menos)\n', a2, 100*(1 - ____/a0));

plot(t, onda, t, y01, t, y10)
legend('sin filtrar', '0,1 Hz', '1,0 Hz')
xlabel('Tiempo (s)'), ylabel('Amplitud (\muV)')

clear; clc; rng(42);
srate = 500;
t     = (-100:399) / srate;         % enteros / srate: t(101) es 0 exacto
onda  = 8 * exp(-((t - 0.400).^2) / (2*0.080^2));
vent  = t >= 0.300 & t <= 0.500;

prom = mean(onda + 25 * randn(40, numel(t)), 1);   % el promedio del día 2

for entrada = {onda, prom}
    x  = entrada{1};
    a0 = mean(x(vent));
    fprintf('\nsin filtrar   : %6.3f uV\n', a0);
    for fc = [0.1 1.0]
        y = pasoalto1(x, fc, srate);
        a = mean(y(vent));
        fprintf('paso alto %3.1f : %6.3f uV  (%.1f%% menos)\n', fc, a, 100*(1 - a/a0));
    end
end

figure; hold on
plot(t, prom)
plot(t, pasoalto1(prom, 0.1, srate))
plot(t, pasoalto1(prom, 1.0, srate))
plot([t(1) t(end)], [0 0], 'k')
legend('sin filtrar', '0,1 Hz', '1,0 Hz')
xlabel('Tiempo (s)'), ylabel('Amplitud (\muV)')

function y = pasoalto1(x, fc, srate)
    a = 1 / (1 + 2*pi*fc/srate);
    y = zeros(size(x));
    for n = 2:numel(x)
        y(n) = a * (y(n-1) + x(n) - x(n-1));
    end
end

Las cinco etapas

ERPLAB impone un orden. No es arbitrario: cada paso necesita lo que dejó el anterior.

Los códigos del paradigma N400

En este experimento cada ensayo tiene tres eventos: aparece una palabra prime, después una palabra target que puede estar semánticamente relacionada con ella o no, y el participante responde si lo están. Los códigos de ERP CORE tienen tres dígitos con un significado cada uno.

El tercer dígito distingue dos listas de palabras que se contrabalancean entre participantes, así que para el análisis se juntan: lo que separa las condiciones son los dos primeros.

El archivo que define tus condiciones

El paso conceptualmente más importante es el segundo, y se controla con un archivo de texto que escribes tú: el bin descriptor file. Cada bloque declara una condición: número de bin, una etiqueta y una secuencia de códigos. El punto marca el evento al que se alinea el tiempo cero.

bin 1
Prime de par relacionado
.{111;112}

bin 2
Prime de par no relacionado
.{121;122}

bin 3
Target relacionado, respuesta correcta
.{211;212}{201}

bin 4
Target no relacionado, respuesta correcta
.{221;222}{201}

El bin 3 se lee así: «alinea el tiempo cero en un 211 o en un 212, y acepta la época sólo si el evento siguiente es un 201». El punto y coma dentro de las llaves significa «cualquiera de estos» —es como se juntan las dos listas de palabras— y el bloque que va después del punto exige un contexto posterior. También se pueden poner llaves antes del punto para exigir un contexto previo: con esa gramática, «rara precedida de frecuente y seguida de respuesta correcta» es una sola línea.

Los primes van en sus propios bins aunque el efecto N400 no esté ahí. Sirven de control: si el prime relacionado y el no relacionado difieren, algo está mal, porque hasta ese momento las dos condiciones son idénticas para el participante.

El pipeline, línea por línea

[ALLEEG, EEG, CURRENTSET, ALLCOM] = eeglab;

raiz = '/ruta/a/ERP_CORE/N400';
s    = 'sub-001';
carp = fullfile(raiz, s);

% --- cargar
EEG = pop_loadset('filename', [s '_N400.set'], 'filepath', carp);

% --- paso alto sobre el continuo
EEG = pop_basicfilter(EEG, 1:33, 'Filter', 'highpass', 'Design', 'butter', ...
    'Cutoff', 0.1, 'Order', 2, 'RemoveDC', 'on', 'Boundary', 'boundary');

% --- 1. EventList
EEG = pop_creabasiceventlist(EEG, 'AlphanumericCleaning', 'on', ...
    'BoundaryNumeric', {-99}, 'BoundaryString', {'boundary'});

% --- 2. asignar a bins
EEG = pop_binlister(EEG, 'BDF', fullfile(raiz, 'bins_n400.txt'), ...
    'IndexEL', 1, 'SendEL2', 'EEG', 'Voutput', 'EEG');

% --- 3. epocar: -200 a 800 ms, línea base = todo el pre-estímulo
EEG = pop_epochbin(EEG, [-200.0 800.0], 'pre');

% --- 4. detectar artefactos (ventana móvil, pico a pico)
EEG = pop_artmwppth(EEG, 'Channel', 1:33, 'Threshold', 100, ...
    'Twindow', [-200 798], 'Windowsize', 200, 'Windowstep', 50, ...
    'Flag', 1, 'Review', 'off');

% --- 5. promediar sólo las épocas buenas
ERP = pop_averager(EEG, 'Criterion', 'good', 'SEM', 'on', 'Compute', 'ERP');

ERP = pop_savemyerp(ERP, 'erpname', [s '_N400'], ...
    'filename', [s '_N400.erp'], 'filepath', carp, 'Warning', 'off');

Dos detalles que valen los cinco minutos que toma entenderlos.

El tercer argumento de pop_epochbin, 'pre', dice que la línea base es todo el período anterior al estímulo. Es el valor por defecto sensato; se puede dar un rango explícito como [-200 0].

Y pop_artmwppth no borra nada. Desliza una ventana de 200 ms por la época, mide la diferencia entre el máximo y el mínimo dentro de esa ventana, y si supera los 100 µV levanta una bandera. Las épocas siguen ahí; es pop_averager con 'Criterion','good' quien decide ignorarlas. Que la marca sea reversible es la razón por la que se puede volver a mirar un rechazo dudoso sin reprocesar todo.

pop_summary_AR_eeg_detection(EEG, '');   % % de rechazo, por bin y total
pop_eegplot(EEG, 1, 1, 1);               % mirar las épocas marcadas
pop_ploterps(ERP, [3 4], 1:33);          % los dos bins de target, todos los canales

iCPz = find(strcmp({ERP.chanlocs.labels}, 'CPz'));   % el sitio donde vive el N400
pop_ploterps(ERP, [3 4], iCPz);

El N400 se mide en CPz —canal 14 en el montaje de ERP CORE, pero búscalo por nombre en vez de confiar en el número— y la ventana estándar es 300 a 500 ms. En el bin 4 (target no relacionado) la onda debería ser más negativa que en el bin 3 en ese tramo. Los bins 1 y 2 deberían quedar prácticamente encima uno del otro.

raiz = '/ruta/a/ERP_CORE/N400';
bdf  = fullfile(raiz, 'bins_n400.txt');

fid = fopen(bdf, 'w');
fprintf(fid, 'bin 1\nPrime de par relacionado\n.{111;112}\n\n');
fprintf(fid, 'bin 2\nPrime de par no relacionado\n.{121;122}\n\n');
fprintf(fid, 'bin 3\nTarget relacionado, respuesta correcta\n.{211;212}{201}\n\n');
fprintf(fid, 'bin 4\nTarget no relacionado, respuesta correcta\n.{221;222}{201}\n');
fclose(fid);

type(bdf)          % míralo antes de usarlo

size(EEG.data)      % antes de epocar: DOS números

EEG = pop_creabasiceventlist(EEG, 'AlphanumericCleaning', 'on', ...
    'BoundaryNumeric', {-99}, 'BoundaryString', {'boundary'});

EEG = pop_binlister(EEG, 'BDF', ____, ...
    'IndexEL', 1, 'SendEL2', 'EEG', 'Voutput', 'EEG');

EEG = pop_epochbin(EEG, [-200.0 ____], 'pre');

size(EEG.data)      % después: TRES números. ¿Qué es cada uno?
EEG.trials

function pipeline_un_sujeto(raiz, s, umbral)
% PIPELINE_UN_SUJETO  De .set continuo a .erp, un participante.
%   raiz   : carpeta del experimento
%   s      : 'sub-001'
%   umbral : criterio de artefactos en µV (100 por defecto)

    if nargin < 3, umbral = 100; end
    carp   = fullfile(raiz, s);
    salida = fullfile(raiz, 'derivados');          % todo lo generado, junto
    if ~exist(salida, 'dir'), mkdir(salida); end

    EEG = pop_loadset('filename', [s '_N400.set'], 'filepath', carp);

    EEG = pop_basicfilter(EEG, 1:33, 'Filter', 'highpass', 'Design', 'butter', ...
        'Cutoff', 0.1, 'Order', 2, 'RemoveDC', 'on', 'Boundary', 'boundary');

    EEG = pop_creabasiceventlist(EEG, 'AlphanumericCleaning', 'on', ...
        'BoundaryNumeric', {-99}, 'BoundaryString', {'boundary'});

    EEG = pop_binlister(EEG, 'BDF', fullfile(raiz, 'bins_n400.txt'), ...
        'IndexEL', 1, 'SendEL2', 'EEG', 'Voutput', 'EEG');

    EEG = pop_epochbin(EEG, [-200.0 800.0], 'pre');

    EEG = pop_artmwppth(EEG, 'Channel', 1:33, 'Threshold', umbral, ...
        'Twindow', [-200 798], 'Windowsize', 200, 'Windowstep', 50, ...
        'Flag', 1, 'Review', 'off');

    pop_summary_AR_eeg_detection(EEG, '');       % queda en la consola

    ERP = pop_averager(EEG, 'Criterion', 'good', 'SEM', 'on', 'Compute', 'ERP');
    ERP = pop_savemyerp(ERP, 'erpname', sprintf('%s_N400_%d', s, umbral), ...
        'filename', sprintf('%s_N400_%d.erp', s, umbral), ...
        'filepath', salida, 'Warning', 'off');
end

El bucle sobre sujetos

Todo lo del día 5 se convierte en una función de un argumento, y el script se reduce a recorrer una lista. La estructura que sigue tiene tres piezas que valen más que el resto: un try/catch para que un sujeto roto no detenga las tres horas de procesamiento, un registro de lo que pasó, y una lista explícita de sujetos en vez de dir('*').

[ALLEEG, EEG, CURRENTSET, ALLCOM] = eeglab;

raiz     = '/ruta/a/ERP_CORE/N400';
sujetos  = {'sub-001','sub-002','sub-003','sub-004','sub-005'};
registro = {};                      % no lo llames log: log() es el logaritmo

for k = 1:numel(sujetos)
    s = sujetos{k};
    fprintf('\n===== %s (%d de %d) =====\n', s, k, numel(sujetos));
    try
        pipeline_un_sujeto(raiz, s);          % la función del día 5
        registro{end+1} = sprintf('%s  OK', s);
    catch err
        registro{end+1} = sprintf('%s  FALLO: %s', s, err.message);
        fprintf(2, 'Falló %s: %s\n', s, err.message);
    end
end

fid = fopen(fullfile(raiz, 'log_procesamiento.txt'), 'w');
fprintf(fid, '%s\n', registro{:});
fclose(fid);

Tres detalles de esas veinte líneas. registro{end+1} = ... agrega un elemento al final de una celda sin saber cuántos van a ser. El 2 de fprintf(2, ...) manda el texto al canal de errores, que MATLAB imprime en rojo: al volver de almorzar, los fallos se ven de lejos. Y err.message es el mensaje que habría aparecido si el error hubiera detenido todo — catch no lo esconde, lo guarda.

Una lista escrita a mano parece un retroceso frente a leer la carpeta automáticamente. No lo es: cuando excluyas un sujeto por calidad de datos, la exclusión queda escrita en el script, con fecha en el control de versiones, en vez de ser una carpeta que alguien movió.

Ondas de diferencia

La resta entre condiciones se hace con la misma gramática de texto que los canales, pero sobre bins. Se ejecuta sobre el objeto ERP, sujeto por sujeto, después de promediar:

ERP = pop_binoperator(ERP, { ...
    'bin5 = bin4 - bin3 label Target no relacionado menos relacionado' });

La diferencia elimina toda la actividad que es idéntica en las dos condiciones y deja sólo lo que las separa. Un aviso de la documentación que conviene tener presente: al crear bins nuevos por operación, ERPLAB pierde las métricas de calidad de datos asociadas, porque no puede estimarlas para una onda derivada.

Gran promedio

pop_gaverager no recibe los ERP: recibe la ruta de un archivo de texto con una ruta .erp por línea. Ese archivo se genera en el mismo script, que es como lo hacen los scripts publicados de ERP CORE.

salida = fullfile(raiz, 'derivados');       % donde quedaron los .erp

lista = fullfile(salida, 'lista_erpsets.txt');
fid = fopen(lista, 'w');
for k = 1:numel(sujetos)
    fprintf(fid, '%s\n', fullfile(salida, sprintf('%s_N400_100.erp', sujetos{k})));
end
fclose(fid);

GA = pop_gaverager(lista, 'ExcludeNullBin', 'on', 'SEM', 'on');
GA = pop_savemyerp(GA, 'erpname', 'GA_N400', ...
    'filename', 'GA_N400.erp', 'filepath', salida);

Qué número medir

Antes de exportar hay que decidir qué se exporta, y la opción intuitiva —la altura del pico— es la peor de las disponibles. Compruébalo con la simulación del día 2: la misma señal de siempre, 40 ensayos, y dos niveles de ruido.

rng(7);
srate = 500;
t     = (-100:399) / srate;         % enteros / srate: t(101) es 0 exacto
senal = 8 * exp(-((t - 0.400).^2) / (2*0.080^2));
vent  = t >= 0.300 & t <= 0.500;

for sd = [25 50]
    p = mean(senal + sd*randn(40, numel(t)), 1);
    fprintf('ruido %2d uV -> pico %5.2f   media %5.2f\n', sd, max(p(vent)), mean(p(vent)));
end
fprintf('verdad       -> pico %5.2f   media %5.2f\n', ...
    max(senal(vent)), mean(senal(vent)));

Tus números van a bailar de corrida en corrida. Promediando dos mil corridas queda así:

El pico medido duplica y triplica el verdadero, y crece cuando crece el ruido. La amplitud media no se mueve. La razón es que buscar un máximo es elegir el punto más alto de la ventana: el ruido que empuja hacia arriba se queda, el que empuja hacia abajo se descarta. Promediar la ventana, en cambio, deja que se cancelen.

La consecuencia es directa y es la razón por la que ERP CORE mide amplitud media: el pico depende de cuánto ruido quedó, y cuánto ruido quedó depende de cuántos ensayos sobrevivieron a la detección de artefactos. Comparar picos entre dos condiciones con distinto número de ensayos aceptados es, en parte, comparar sus niveles de ruido. Con datos reales —filtrados en paso bajo, con ruido menos abrupto que el de esta simulación— la exageración es menor, pero nunca es cero y nunca es simétrica.

Sacar los números

El último paso del análisis en MATLAB es también el más corto: extraer una medida por sujeto, bin y canal, y escribirla a un archivo que abra cualquier programa estadístico.

% ALLERP es la lista de ERPsets que hay en memoria; se llena con pop_loaderp.
% Ventana 300-500 ms; bins 3, 4 y 5; el canal, buscado por nombre.
iCPz = find(strcmp({ALLERP(1).chanlocs.labels}, 'CPz'));

[ALLERP, Amp, Lat] = pop_geterpvalues(ALLERP, [300 500], [3 4 5], iCPz, ...
    'Measure',    'meanbl', ...      % amplitud media, con línea base
    'Baseline',   'pre', ...
    'Erpsets',    1:numel(sujetos), ...
    'Filename',   fullfile(raiz, 'amplitudes_n400.txt'), ...
    'Foutput',    'erpset', ...      % formato largo: una fila por medición
    'Fracreplace','NaN', ...
    'Resolution', 3, ...
    'SendtoWorkspace', 'on', ...
    'Warning',    'off');

Los primeros cuatro argumentos son la ventana en ms, los bins, los canales, y el resto son opciones. Los valores de arriba son los que recomienda ERP CORE para el N400: amplitud media entre 300 y 500 ms en CPz. 'Measure' acepta también 'peaklatbl' (latencia del pico local) y 'areap' (área positiva), entre otras. 'Foutput','erpset' entrega formato largo —una fila por sujeto × bin × canal—, que es el que quiere R; el formato ancho pone todas las medidas de un sujeto en una línea, que es el que quieren las tablas dinámicas.

Lo que evita que te odies en seis meses

• Los datos crudos son de sólo lectura. Cada etapa escribe un archivo nuevo con un sufijo (_elist, _be, _ar). Nunca sobrescribes lo que llegó del equipo.
• El script es el método. Si el pipeline está en un .m versionado, la sección de métodos se escribe leyéndolo, y un revisor puede pedirlo.
• Semilla fija con rng(42) en cualquier cosa que use aleatoriedad, incluida ICA.
• Reporta la calidad de los datos, no sólo el número de ensayos rechazados. ERPLAB calcula el error estándar de medición (SME), que estima directamente cuánto ruido queda en la medida que vas a analizar.
• Guarda el log. Un archivo de texto con qué sujeto se procesó, cuándo, con qué parámetros y con qué porcentaje de rechazo vale más que la memoria.

ERP = pop_binoperator(ERP, { ...
    'bin5 = bin4 - bin3 label Target no relacionado menos relacionado' });

[ALLEEG, EEG, CURRENTSET, ALLCOM] = eeglab;

raiz    = '/ruta/a/ERP_CORE/N400';
salida  = fullfile(raiz, 'derivados');
sujetos = {'sub-001','sub-002','sub-003','sub-004','sub-005'};

archivos = cell(1, numel(sujetos));
for k = 1:numel(sujetos)
    archivos{k} = sprintf('%s_N400_100.erp', sujetos{k});
end

[ERP, ALLERP] = pop_loaderp('filename', archivos, 'filepath', salida);

numel(ALLERP)              % ¿cargaron los cinco?
ALLERP(1).bindescr         % las etiquetas de los bins del primero

iCPz = find(strcmp({ALLERP(1).chanlocs.labels}, '____'));

[ALLERP, Amp, Lat] = pop_geterpvalues(ALLERP, [____ ____], [3 4 5], iCPz, ...
    'Measure',    'meanbl', ...
    'Baseline',   'pre', ...
    'Erpsets',    1:numel(sujetos), ...
    'Filename',   fullfile(salida, 'amplitudes_n400.txt'), ...
    'Foutput',    'erpset', ...
    'Fracreplace','NaN', ...
    'Resolution', 3, ...
    'SendtoWorkspace', 'on', ...
    'Warning',    'off');

size(Amp)                  % ¿cuántas mediciones salieron, y de qué forma?

[ALLEEG, EEG, CURRENTSET, ALLCOM] = eeglab;

raiz     = '/ruta/a/ERP_CORE/N400';
salida   = fullfile(raiz, 'derivados');
sujetos  = {'sub-001','sub-002','sub-003','sub-004','sub-005'};
registro = {};

% --- 1. procesar cada sujeto, sin que uno malo detenga el lote
for k = 1:numel(sujetos)
    s = sujetos{k};
    fprintf('\n===== %s (%d de %d) =====\n', s, k, numel(sujetos));
    try
        pipeline_un_sujeto(raiz, s, 100);
        registro{end+1} = sprintf('%s  OK', s);
    catch err
        registro{end+1} = sprintf('%s  FALLO: %s', s, err.message);
        fprintf(2, 'Falló %s: %s\n', s, err.message);
    end
end

fid = fopen(fullfile(salida, 'log_procesamiento.txt'), 'w');
fprintf(fid, '%s\n', registro{:});
fclose(fid);

% --- 2. cargar los .erp
archivos = cell(1, numel(sujetos));
for k = 1:numel(sujetos)
    archivos{k} = sprintf('%s_N400_100.erp', sujetos{k});
end
[ERP, ALLERP] = pop_loaderp('filename', archivos, 'filepath', salida);

% --- 3. gran promedio, desde una lista en texto
lista = fullfile(salida, 'lista_erpsets.txt');
fid = fopen(lista, 'w');
for k = 1:numel(archivos)
    fprintf(fid, '%s\n', fullfile(salida, archivos{k}));
end
fclose(fid);

GA = pop_gaverager(lista, 'ExcludeNullBin', 'on', 'SEM', 'on');
GA = pop_savemyerp(GA, 'erpname', 'GA_N400', ...
    'filename', 'GA_N400.erp', 'filepath', salida);

% --- 4. medir y exportar
iCPz = find(strcmp({ALLERP(1).chanlocs.labels}, 'CPz'));
pop_ploterps(GA, [3 4 5], iCPz);

[ALLERP, Amp, Lat] = pop_geterpvalues(ALLERP, [300 500], [3 4 5], iCPz, ...
    'Measure', 'meanbl', 'Baseline', 'pre', ...
    'Erpsets', 1:numel(sujetos), ...
    'Filename', fullfile(salida, 'amplitudes_n400.txt'), ...
    'Foutput', 'erpset', 'Fracreplace', 'NaN', 'Resolution', 3, ...
    'SendtoWorkspace', 'on', 'Warning', 'off');

Diagnóstico

Para cada fragmento, responde tres cosas antes de abrir la respuesta: qué está mal, qué resultado produce, y con qué línea de código lo comprobarías en treinta segundos.

1 · El promedio que no promedia lo que crees

% epocas viene de otro script, con forma tiempo × ensayos
size(epocas)                  % 500   40

promedio = mean(epocas, 1);
plot(promedio)

2 · La ventana vacía

vent = ERP.times >= 0.3 & ERP.times <= 0.5;
amp  = mean(ERP.bindata(iCPz, vent, 4));
fprintf('N400: %.2f uV\n', amp);

3 · Un sujeto que son siete

sujetos = 'sub-001';

for k = 1:numel(sujetos)
    s = sujetos(k);
    fprintf('procesando %s\n', s);
    pipeline_un_sujeto(raiz, s);
end

4 · El filtro en el lugar equivocado

EEG = pop_epochbin(EEG, [-200.0 800.0], 'pre');

EEG = pop_basicfilter(EEG, 1:33, 'Filter', 'highpass', 'Design', 'butter', ...
    'Cutoff', 0.1, 'Order', 2, 'RemoveDC', 'on');

5 · La ventana que encontró el efecto

dif = ERP.bindata(iCPz, :, 5);        % la onda de diferencia

[~, i] = max(abs(dif));               % ¿dónde es más grande?
vent   = i-25 : i+25;                 % 100 ms alrededor

amp = mean(dif(vent));
fprintf('efecto: %.2f uV\n', amp);

Dos decisiones

Primera. Un colega mide amplitud de pico y encuentra que su grupo clínico tiene picos más chicos que los controles. En el mismo trabajo reporta que en el grupo clínico rechazó bastantes más ensayos por artefactos. ¿En qué dirección empuja ese desbalance a su resultado?

Segunda. Un revisor pide aplicar un paso bajo de 30 Hz al ERP promediado «para que las figuras se vean más limpias». Tú mides amplitud media entre 300 y 500 ms. ¿Cuánto puede cambiar tu resultado?

Lo que deberías poder hacer sin mirar

Si la semana funcionó, estas cinco cosas se pueden hacer con MATLAB abierto y la página cerrada. Si alguna se atasca, ya sabes a qué día volver.

1. Escribir, de memoria, la corrección de línea base de una época: tres líneas.
2. Explicar, con la fórmula, por qué cuadruplicar los ensayos duplica el SNR — y decir cuántos ensayos harían falta para duplicarlo de nuevo.
3. Nombrar las cinco etapas del pipeline de ERPLAB en orden y decir qué se rompe si intercambias dos.
4. Dar un valor de un parámetro de filtrado que pueda fabricar un componente que no existe, y explicar por dónde aparece.
5. Escribir el párrafo de métodos de tu propio análisis leyendo únicamente tu script.

Y el reflejo general: cuando algo no cuadra, escribe el nombre de la variable sin punto y coma y mírala. size, class, min, max y whos resuelven más problemas que cualquier búsqueda en internet.

MATLAB base

EEGLAB y ERPLAB