El algoritmo LMS

El MSEG baja por la superficie del MSE con el gradiente exacto, $2\mathrm{\Phi }c-2\alpha$, y para eso necesita $\mathrm{\Phi }$ y $\alpha$: las estadísticas del canal, que es justo lo que el receptor no conoce. El LMS —la cátedra lo llama también SG, gradiente estocástico— hace la jugada barata: en cada símbolo usa el error de ese símbolo como si fuera el promedio. Cada paso apunta un poco mal, en promedio apunta bien, y el precio es que los taps nunca se quedan quietos. En la cadena canal → AGC → ecualizador → slicer, el LMS es lo que mueve los taps del ecualizador con el error que vuelve del slicer.

01 · La definición

El ecualizador es un FIR transversal de $N$ taps. En cada símbolo $k$, con $r_k$ las $N$ muestras de $R_k$ que hay en la línea de retardos y $A_k$ el símbolo de referencia:

$$Y_k=c_k^T{r}_k,E_k=A_k-Y_k,c_{k+1}=c_k+\beta E_k{r}_k^{\ast }$$

$\beta$ es el paso de adaptación. $A_k$ es el símbolo conocido mientras dura el entrenamiento, o la decisión del slicer después; de dónde sale es tema de entrenamiento y modo dirigido por decisiones, y al algoritmo le da igual. Con señales reales, como una PAM, el conjugado no hace nada.

Mirado tap por tap, que es como lo dibuja la filmina: el coeficiente ${[c_k]}_j$, el que multiplica a la muestra $R_{k-j}$, tiene su propio acumulador, que suma $\beta E_k{R}_{k-j}^{\ast }$ símbolo tras símbolo. Integra la correlación entre el error y lo que ese tap está viendo. Si todavía hay correlación, al tap le queda información sin usar y se mueve; en el óptimo esa correlación es cero en promedio, que es el principio de ortogonalidad.

El script de la clase, SG.m, escribe el error al revés, error = y - y_ref, y por eso resta: c = c - step*error*buffer. Son dos signos cambiados a la vez y es el mismo algoritmo, y step es exactamente el $\beta$ de la filmina: el gradiente de $|E_k{|}^2$ trae un 2 que se cancela con el $\beta /2$ de la recursión (sección 02). No pasa lo mismo en MSEG.m, cuyo beta multiplica al gradiente entero y vale la mitad del $\beta$ de la filmina. Lo que no se puede es mezclar el error de una convención con el signo de la otra: los taps suben por la superficie en lugar de bajar.

02 · Del MSEG al LMS, paso a paso

El MSEG de la filmina corrige los taps en contra del gradiente del MSE:

$$c_{j+1}=c_j-\frac{\beta }{2}{\nabla }_{c}E[|E_k{|}^2]=(I-\beta \mathrm{\Phi })c_j+\beta \alpha$$

con $\mathrm{\Phi }=E[r_k^{\ast }{r}_k^T]$ y $\alpha =E[A_k{r}_k^{\ast }]$. La filmina del LMS arranca de ahí: en un canal real $\mathrm{\Phi }$ no se conoce de antemano, así que la recursión no se puede evaluar. La esperanza no se puede calcular, pero el módulo al cuadrado de un error ya observado sí.

Paso 1: sacar la esperanza. Es el mismo desarrollo del cuadrado que da el MSE, sin el $E[·]$ de afuera:

$$|E_k{|}^2=|A_k{|}^2-2\mathrm{\Re }(c^H{A}_k{r}_k^{\ast })+c^H{r}_k^{\ast }{r}_k^{T}c$$

Donde el MSE tenía $\alpha$ ahora está $A_k{r}_k^{\ast }$, y donde tenía $\mathrm{\Phi }$ está $r_k^{\ast }{r}_k^T$: las dos cosas que había que promediar, tomadas en un solo símbolo.

Paso 2: derivar término a término. Con las mismas reglas que dieron $2\mathrm{\Phi }c-2\alpha$, cada término se deriva igual que su pariente con esperanza. El primero no depende de $c$:

$${\nabla }_c|E_k{|}^2=-2A_k{r}_k^{\ast }+2r_k^{\ast }{r}_k^{T}c$$

Paso 3: sacar factor común. Los dos términos llevan $2r_k^{\ast }$, y lo que queda adentro es el error, porque $r_k^{T}c=c^T{r}_k=Y_k$:

$${\nabla }_c|E_k{|}^2=-2r_k^{\ast }(A_k-r_k^{T}c)=-2E_k{r}_k^{\ast }$$

Paso 4: meterlo en la recursión. Donde el MSEG tenía el gradiente del MSE va el de $|E_k{|}^2$, y el $\beta /2$ se come el 2:

$$c_{k+1}=c_k-\frac{\beta }{2}{\nabla }_c|E_k{|}^2=c_k+\beta E_k{r}_k^{\ast }$$

El índice cambió de nombre y de significado. En el MSEG, $j$ contaba iteraciones y cada una necesitaba el MSE, o sea un promedio sobre muchos símbolos; en el LMS, $k$ es el tiempo, y hay una actualización por símbolo. La filmina lo resume así: se sustituye el promedio del ensamble por el promedio del tiempo, y es sólo una aproximación, el precio que se paga por no conocer el canal.

03 · Lo que cuesta

Por símbolo, $2N+1$ multiplicaciones y $2N$ sumas, según la filmina: $N$ multiplicaciones para filtrar, una para $\beta E_k$ y $N$ para multiplicar ese escalar por el vector. No estima $\mathrm{\Phi }$, no invierte nada, no resuelve ningún sistema. Para comparar, el script MSEG.m estima el gradiente con lotes: por cada actualización filtra 1000 muestras con el $c$ actual y otras cinco veces con un tap corrido, y promedia seis veces los cuadrados. Con 5 taps son $6·(5000+1000)=36000$ multiplicaciones por actualización. El LMS con los mismos 5 taps hace 11.

04 · En promedio, el MSEG

El gradiente de un símbolo apunta mal, pero no está sesgado. Con $c$ fijo, su esperanza es exactamente el gradiente del MSEG:

$$E[-2E_k{r}_k^{\ast }]=-2E[A_k{r}_k^{\ast }]+2E[r_k^{\ast }{r}_k^T]c=2\mathrm{\Phi }c-2\alpha$$

La filmina lo marca como «estimador sin sesgo del MSEG», y se puede hacer porque el error es cuadrático y el gradiente y la esperanza son lineales: se pueden intercambiar. Tomando esperanza en la actualización:

$$E[c_{k+1}]=E[(I-\beta r_k^{\ast }{r}_k^T)c_k]+\beta E[A_k{r}_k^{\ast }]\approx (I-\beta \mathrm{\Phi })E[c_k]+\beta \alpha$$

El $\approx$ es donde está la trampa. $c_k$ se armó con las muestras de antes, y $r_k$ comparte con ellas todas sus muestras menos una, así que no son independientes y la esperanza del producto no se parte. La filmina lo resuelve con el paso chico: si $\beta$ es chico, $c_k$ cambia despacio comparado con la entrada y se lo puede tratar como constante. Dentro de esa aproximación, el promedio de las trayectorias del LMS es la trayectoria del MSEG: los mismos modos, el mismo punto de llegada $c_{\text{opt}}$ y la misma condición para que la media converja, $0<\beta <2/{\lambda }_{\text{max}}$.

Para medirlo hace falta un ecualizador de verdad. El de toda esta nota tiene $N=11$ taps y ecualiza una PAM-2 que pasó por el canal $p=[\mathrm{0,2}1\mathrm{0,3}]$ —un precursor, el cursor y un poscursor— con ruido de varianza 0,01, buscando el símbolo con un retardo de 6. La potencia de entrada es ${\phi }_0=\mathrm{1,14}$, el mínimo de Wiener es ${\xi }_{\text{min}}=-\mathrm{18,26}$ dB, y la entrada está correlacionada: la dispersión de autovalores de $\mathrm{\Phi }$ es 7,9. Arrancando de $c=0$, sobre 2000 realizaciones de 3000 símbolos, la mayor distancia entre el promedio de los $c_k$ y el $c_k$ del MSEG con el mismo $\beta$:

Distancia euclídea en el espacio de los 11 taps, contra $\Vert c_{\text{opt}}\Vert =\mathrm{1,208}$. Con paso chico el promedio se pega al MSEG, y cualquier realización suelta anda lejos: «en promedio» quiere decir literalmente en promedio. Con el paso grande se separan los dos, porque la aproximación que los hacía iguales deja de valer.

05 · Nunca se queda quieto

En $c_{\text{opt}}$ el gradiente vale cero en promedio, pero el de un símbolo no. El error del óptimo, el que ningún ecualizador de $N$ taps saca, sigue ahí, y en cada símbolo empuja a los taps un poco para algún lado. La filmina lo dice así: el LMS no converge al vector óptimo sino a su entorno, y cuanto más grande es $\beta$, más grande es la fluctuación. Para medir ese entorno define el error de coeficientes

$$q_k=c_k-c_{\text{opt}},E[|E_k{|}^2]={\xi }_{\text{min}}+E[q_k^H\mathrm{\Phi }{q}_k]$$

El segundo término es el exceso de MSE: lo que el error supera al mínimo porque los taps no están en el óptimo. Para poder seguir, la filmina supone entrada blanca de media nula, $\mathrm{\Phi }={\phi }_{0}I$ con ${\phi }_0=E[|R_k{|}^2]$, y el exceso queda ${\phi }_{0}E\Vert q_k{\Vert }^2$.

Paso 1: restar el óptimo. Restando $c_{\text{opt}}$ a los dos lados de la actualización, y escribiendo el error como el del óptimo menos lo que agrega $q_k$, $E_k=E_k^{\text{opt}}-r_k^T{q}_k$:

$$\Vert q_{k+1}{\Vert }^2=\Vert q_k{\Vert }^2+2\beta \mathrm{\Re }(E_k{q}_k^H{r}_k^{\ast })+{\beta }^2|E_k{|}^2\Vert r_k{\Vert }^2$$

Paso 2: el término lineal tira hacia el óptimo. Con $q_k$ independiente de $r_k$, la parte de $E_k^{\text{opt}}$ muere por ortogonalidad y queda $-2\beta {\phi }_{0}E\Vert q_k{\Vert }^2$. Es lo que haría el gradiente exacto.

Paso 3: el término cuadrático empuja al azar. Aproximando la esperanza del producto por el producto de esperanzas, ${\beta }^2({\xi }_{\text{min}}+{\phi }_{0}E\Vert q_k{\Vert }^2)·N{\phi }_0$. Es el ruido del gradiente, y tiene una parte que no depende de $q_k$ y no se va nunca.

Paso 4: juntar. Es la ecuación en diferencias de la filmina:

$$E\Vert q_{k+1}{\Vert }^2=\gamma E\Vert q_k{\Vert }^2+{\beta }^{2}N{\phi }_0{\xi }_{\text{min}},\gamma =1-2\beta {\phi }_0+{\beta }^{2}N{\phi }_0^2$$

La filmina la escribe $E\Vert q_{k+1}\Vert =\gamma E\Vert q_k\Vert +\beta N{\phi }_0{\xi }_{\text{min}}$, sin los cuadrados y con $\beta$ en lugar de ${\beta }^2$. Así no puede ser: con ese término el exceso final sería $N{\phi }_0{\xi }_{\text{min}}/(2-\beta N{\phi }_0)$, que no se achica al achicar $\beta$, y con el paso óptimo el MSE terminaría en $(1+N{\phi }_0){\xi }_{\text{min}}$ y no en el $2{\xi }_{\text{min}}$ que dice la misma filmina dos páginas después. La simulación de la sección 07 da el $2{\xi }_{\text{min}}$.

Paso 5: el piso. La recursión es de primer orden: $E\Vert q_k{\Vert }^2$ cae como ${\gamma }^k$ hacia el valor que la deja quieta, ${\beta }^{2}N{\phi }_0{\xi }_{\text{min}}/(1-\gamma )$. Multiplicado por ${\phi }_0$ es el exceso que queda para siempre, y dividido por ${\xi }_{\text{min}}$ es el desajuste:

$$M=\frac{\beta N{\phi }_0}{2-\beta N{\phi }_0},E[|E_{\infty }{|}^2]={\xi }_{\text{min}}(1+M)=\frac{2{\xi }_{\text{min}}}{2-\beta N{\phi }_0}$$

Crece con $\beta$ y con $N$: pasos más grandes, más temblor, y más taps temblando. Con paso chico, $M\approx \beta N{\phi }_0/2$.

06 · Cuánto puede valer el paso

Todo lo que importa de la recursión está en $\gamma$, que es una parábola en $\beta$. El término $-2\beta {\phi }_0$ acelera; el ${\beta }^{2}N{\phi }_0^2$, que es el castigo por usar un gradiente ruidoso, frena, y a partir de cierto punto gana.

La estabilidad. La norma del error de coeficientes decrece si $|\gamma |<1$. Como el mínimo de $\gamma$ es positivo, alcanza con $\gamma <1$:

$$0<\beta <\frac{2}{N{\phi }_0}=2{\beta }_{\text{opt}}$$

Es mucho más estricta que la condición de la media, $\beta <2/{\lambda }_{\text{max}}$: como $N{\phi }_0$ es la traza de $\mathrm{\Phi }$, la suma de los $N$ autovalores, nunca es menor que ${\lambda }_{\text{max}}$, y con entrada blanca es $N$ veces más grande. La media de los taps puede ir tranquila al óptimo mientras cada realización se dispara: cada actualización mete $N$ empujones al azar, y con paso grande el ruido le gana a la corrección.

El paso óptimo. Derivando, $d\gamma /d\beta =-2{\phi }_0+2\beta N{\phi }_0^2$, que se anula en

$${\beta }_{\text{opt}}=\frac{1}{N{\phi }_0},\gamma ({\beta }_{\text{opt}})=1-\frac{1}{N}$$

Es el paso que más rápido achica el error de coeficientes, y la parábola es simétrica alrededor de él: vuelve a valer 1 justo en $2{\beta }_{\text{opt}}$. Conviene medir el paso en esa unidad, $\beta /{\beta }_{\text{opt}}=\beta N{\phi }_0$, que es justo lo que aparece en el desajuste: $M=f/(2-f)$ con $f=\beta /{\beta }_{\text{opt}}$, y con el paso óptimo $M=1$: el MSE termina en el doble del mínimo.

En las tres fórmulas $\beta$ aparece multiplicado por ${\phi }_0$, la potencia de entrada. Si la señal llega más fuerte, el mismo $\beta$ se comporta como un paso más grande y el borde de la estabilidad se acerca; por eso la filmina termina pidiendo un control automático de ganancia delante del ecualizador. Y el borde también baja con $N$: un paso que andaba con 11 taps puede no andar con 41.

La constante de tiempo. El exceso cae como ${\gamma }^k$; definiendo $\tau$ por ${\gamma }^{\tau }=1/e$,

$$\tau =\frac{-1}{\ln \gamma }\approx \frac{1}{\beta {\phi }_0(2-\beta N{\phi }_0)}\approx \frac{1}{2\beta {\phi }_0}$$

La última es la de la filmina y vale para paso chico. La filmina la usa también con ${\beta }_{\text{opt}}$ y obtiene $\tau \approx N/2$, pero ahí el paso ya no es chico: el término cuadrático es la mitad del lineal. Con su propio $\gamma =1-1/N$, $\tau =-1/\ln (1-1/N)\approx N-1/2$. Medido con entrada blanca y 11 taps, sobre 4000 realizaciones, el exceso cae por debajo de $1/e$ en la muestra 11, contra 10,49 de $\gamma$ y 5,5 de $N/2$; con la mitad del paso, en la 15, contra 14,16 y 11. La conclusión de la filmina queda en pie: el mejor $\tau$ crece con el número de taps, sólo que es del orden de $N$ y no de $N/2$.

El borde, medido. Todo esto supone entrada blanca. Con el ecualizador de la sección 04, cuya entrada no lo es, sobre 500 realizaciones de 3000 símbolos: el MSE promedio entre los símbolos 1501 y 3000, cuántas realizaciones tuvieron en algún momento un tap por encima de 10 —el óptimo entero mide $\Vert c_{\text{opt}}\Vert =\mathrm{1,208}$— y cuántas terminaron con alguno por encima de 100:

El borde de la cuenta es $2{\beta }_{\text{opt}}$; con este canal el MSE promedio ya se dispara en $\mathrm{1,3}{\beta }_{\text{opt}}$. Y la tabla muestra cómo se rompe el LMS cerca del borde: no diverge de golpe. Una realización se escapa lejos del óptimo y casi siempre vuelve —con $\mathrm{1,5}{\beta }_{\text{opt}}$, de 263 que se escaparon sólo 2 terminan arriba de 100—, pero cada escapada pesa en el promedio. Contar las realizaciones que terminan mal esconde el problema; el MSE promedio, no. Cerca del borde no conviene andar.

07 · El compromiso, medido

Un solo número, $\beta$, decide dos cosas que tiran para lados opuestos: qué tan rápido llega el ecualizador y qué tan cerca del mínimo se queda. Con el ecualizador de 11 taps, sobre 1000 realizaciones de 6000 símbolos, el MSE asintótico promediado entre los símbolos 3001 y 6000, y la primera muestra en que la curva promediada llega a 1 dB de ese piso:

La columna de la velocidad va como $1/\beta$: veinte veces más paso, 23 veces menos muestras. El piso sube: con el paso óptimo el MSE queda 4,55 dB arriba de ${\xi }_{\text{min}}$, y con $\mathrm{0,05}{\beta }_{\text{opt}}$, a 0,16 dB. Es la última filmina del tema dicha en números: si lo que importa es el piso, el paso más chico posible; si lo que importa es la velocidad, el óptimo, pagando con el MSE.

La fórmula del desajuste acierta en dB con paso chico y es optimista con paso grande: con ${\beta }_{\text{opt}}$ predice que el MSE se duplica, y casi se triplica. No es que la fórmula esté mal. Sale de dos hipótesis —entrada blanca, y un error que no depende de lo que hay en la línea de retardos— y este canal no cumple ninguna: la entrada tiene una dispersión de autovalores de 7,9, y la ISI residual que queda en el error está hecha de los mismos símbolos que llenan $r_k$. Con el canal ideal, $p=1$, la entrada es blanca y la fórmula da en el clavo: 0,334 medido contra 0,333 con $\mathrm{0,5}{\beta }_{\text{opt}}$, y 0,998 contra 1 con ${\beta }_{\text{opt}}$.

Con $\mathrm{0,02}{\beta }_{\text{opt}}$ la realización azul casi no se despega de la línea cortada, y el promedio de las 400, si se lo prende, cae encima: el LMS con paso chico es el MSEG. El precio es el tiempo, más de doscientas muestras para llegar a 1 dB del piso, y el premio es que el piso queda pegado a ${\xi }_{\text{min}}$. Subiendo el paso, la llegada se acorta y la realización se vuelve un garabato alrededor del óptimo, cada vez más ancho: el piso de la curva de abajo se despega de ${\xi }_{\text{min}}$. Pasado $\mathrm{0,5}{\beta }_{\text{opt}}$ el piso medido se separa también de la línea de la filmina. Alrededor del paso óptimo empiezan a escaparse realizaciones —algún tap pasa de 10—, y con $\mathrm{1,2}{\beta }_{\text{opt}}$ se escapa una de cada cuatro: el MSE promedio pasa de 0 dB, que es lo que da no ecualizar nada ($c=0$). Todo eso, bastante antes del $2{\beta }_{\text{opt}}$ de la cuenta.

08 · La curva de aprendizaje

Es el MSE en dB contra $k$, y la primera figura que hay que mirar después de correr cualquier cosa adaptativa: si baja y se aplana, convergió; si se aplana muy arriba, el paso es grande o faltan taps; si todavía baja al final, el paso es chico o la simulación es corta; si sube en recta, diverge. En la teoría, la esperanza es sobre realizaciones.

Una sola corrida es ruido. Con el ecualizador de 11 taps y $\mathrm{0,5}{\beta }_{\text{opt}}$, en régimen, $|E_k{|}^2$ de una corrida anda entre −33,9 dB y −12,0 dB (los percentiles 10 y 90), con la mediana en −19,7 dB. El MSE, promediando 1000 corridas, es −16,48 dB: ni siquiera la mediana de una corrida se le parece, porque la distribución del cuadrado es asimétrica. Por eso la curva se promedia sobre muchas realizaciones, o, ya en régimen, sobre muchas muestras de una.

Promediar en el tiempo tiene su trampa. SG.m dibuja el promedio móvil de $|E_k{|}^2$ sobre 2560 muestras. Para leer el piso sirve. Para ver la convergencia, no: con el paso del script la constante de tiempo medida es de 8 muestras, y durante las primeras 2560 la ventana arrastra el error grande del arranque. La curva queda alta hasta la 2560 y recién ahí baja, como si el LMS hubiera tardado eso en converger.

09 · El LMS normalizado, en una línea

Como lo que manda es $\beta N{\phi }_0$, se puede dividir el paso por la energía que hay en la línea de retardos, $\Vert r_k{\Vert }^2\approx N{\phi }_0$: $c_{k+1}=c_k+(\mu /\Vert r_k{\Vert }^2)E_k{r}_k^{\ast }$. Es el NLMS: $\mu$ hace de $\beta /{\beta }_{\text{opt}}$ y el paso deja de depender del nivel de la señal. No está en las filminas: la cátedra resuelve el mismo problema del otro lado, fijando ${\phi }_0$ con el AGC.

10 · El análisis en MATLAB

SG.m, el script de la clase, no ecualiza: identifica. La misma entrada blanca gaussiana pasa por un FIR «desconocido», h_ref = fir1(4,0.3), de 5 taps, y por el FIR adaptativo, y el error entre las dos salidas —con un ruido de desvío 0,0025 sumado a la referencia— maneja el LMS. Con ${\phi }_0=1$, ${\beta }_{\text{opt}}=\mathrm{0,2}$, y el step = 0.1 del script es la mitad. Corre tal cual en R2026a y los taps terminan pegados a los de referencia. La filmina predice un piso 1,25 dB arriba del ruido ($M=1/3$); promediando 2000 realizaciones da 1,45 dB ($M=\mathrm{0,40}$), porque la entrada gaussiana tampoco cumple del todo la hipótesis del paso 3.

El análisis del ecualizador de esta nota. Corre las realizaciones en paralelo: cada fila de las matrices es una corrida, y el bucle en $k$ es el de SG.m. Lo que figura como comentario es la salida real de correrlo en R2026a.

%% 1. El canal, la potencia de entrada y el paso óptimo de la filmina
rng(1)
p  = [0.2 1 0.3];            % canal discreto equivalente: precursor, cursor, poscursor
s2 = 0.01;                   % varianza del ruido Z_k
N  = 11;  d = 6;             % taps del ecualizador y retardo del símbolo buscado
phi0     = sum(p.^2) + s2;   % potencia de entrada, E|R_k|^2
beta_opt = 1/(N*phi0);
% Wiener sólo para tener contra qué medir: el LMS no usa nada de esto
fi   = [conv(p, fliplr(p)), zeros(1, N)];
fi   = fi(numel(p):numel(p)+N-1);  fi(1) = fi(1) + s2;
Phi  = toeplitz(fi);
alfa = zeros(N, 1);  alfa(d+1-(0:numel(p)-1)) = p;
c_opt  = Phi \ alfa;
xi_min = 1 - alfa' * c_opt;
fprintf('phi0 = %.4f   beta_opt = %.5f   xi_min = %.2f dB\n', phi0, beta_opt, 10*log10(xi_min))
% phi0 = 1.1400   beta_opt = 0.07974   xi_min = -18.26 dB

%% 2. Una corrida contra el promedio de mil
[mse, ~, e1] = lms(p, s2, N, d, 0.5*beta_opt, 1000, 6000);
v = 10*log10(e1(3001:end).^2);
fprintf('una corrida, |E_k|^2 en régimen: 10 %% = %.1f dB, mediana = %.1f dB, 90 %% = %.1f dB\n', ...
        prctile(v, 10), prctile(v, 50), prctile(v, 90))
fprintf('promedio de 1000 corridas, en régimen: %.2f dB\n', 10*log10(mean(mse(3001:end))))
figure; semilogx(10*log10(e1.^2), 'Color', [.7 .7 .7]); hold on
semilogx(10*log10(mse), 'LineWidth', 1.5); yline(10*log10(xi_min), '--');
xlabel('k'); ylabel('|E_k|^2 [dB]'); legend('una corrida', 'promedio de 1000', '\xi_{min}')
% una corrida, |E_k|^2 en régimen: 10 % = -33.9 dB, mediana = -19.7 dB, 90 % = -12.0 dB
% promedio de 1000 corridas, en régimen: -16.48 dB

%% 3. El compromiso: barrido de beta
fprintf('  beta/beta_opt  k(1 dB)  MSE medido  predicho   M medido  M predicho\n')
for f = [0.05 0.1 0.25 0.5 1]
    mse  = lms(p, s2, N, d, f*beta_opt, 1000, 6000);
    piso = mean(mse(3001:end));
    k1   = find(mse <= piso*10^0.1, 1);
    fprintf('  %5.2f  %8d  %8.2f dB  %6.2f dB  %7.3f  %7.3f\n', f, k1, ...
            10*log10(piso), 10*log10(xi_min*2/(2-f)), piso/xi_min - 1, f/(2-f))
end
%   beta/beta_opt  k(1 dB)  MSE medido  predicho   M medido  M predicho
%    0.05      1517    -18.10 dB  -18.15 dB    0.038    0.026
%    0.10       781    -17.94 dB  -18.04 dB    0.075    0.053
%    0.25       313    -17.45 dB  -17.68 dB    0.205    0.143
%    0.50       142    -16.48 dB  -17.01 dB    0.505    0.333
%    1.00        67    -13.71 dB  -15.25 dB    1.851    1.000

%% 4. La estabilidad: el MSE promedio y las realizaciones que se escapan
for f = [1 1.1 1.2 1.3 1.4 1.5 1.6]
    [mse, ~, ~, c, ~, pico] = lms(p, s2, N, d, f*beta_opt, 500, 3000);
    fprintf('beta = %.1f beta_opt: MSE en régimen %6.2f dB, %3d de 500 pasan de 10, %3d terminan arriba de 100\n', ...
            f, 10*log10(mean(mse(1501:end))), sum(pico > 10), sum(~all(abs(c) < 1e2, 2)))
end
% beta = 1.0 beta_opt: MSE en régimen -13.66 dB,   0 de 500 pasan de 10,   0 terminan arriba de 100
% beta = 1.1 beta_opt: MSE en régimen -12.68 dB,   0 de 500 pasan de 10,   0 terminan arriba de 100
% beta = 1.2 beta_opt: MSE en régimen -10.95 dB,   0 de 500 pasan de 10,   0 terminan arriba de 100
% beta = 1.3 beta_opt: MSE en régimen  -3.36 dB,  11 de 500 pasan de 10,   0 terminan arriba de 100
% beta = 1.4 beta_opt: MSE en régimen  34.24 dB,  72 de 500 pasan de 10,   1 terminan arriba de 100
% beta = 1.5 beta_opt: MSE en régimen  43.87 dB, 263 de 500 pasan de 10,   2 terminan arriba de 100
% beta = 1.6 beta_opt: MSE en régimen  95.50 dB, 461 de 500 pasan de 10,   4 terminan arriba de 100

%% 5. En promedio, el MSEG
for f = [0.1 0.5 1]
    beta = f*beta_opt;  K = 3000;
    [~, c_media, ~, ~, c_una] = lms(p, s2, N, d, beta, 2000, K);
    c_mseg = zeros(K, N);  c = zeros(N, 1);
    for k = 1:K, c_mseg(k,:) = c.';  c = c - beta*(Phi*c - alfa); end
    fprintf('beta = %.1f beta_opt: max||media - MSEG|| = %.3f   max||una - MSEG|| = %.3f   ||c_opt|| = %.3f\n', ...
            f, max(vecnorm(c_media - c_mseg, 2, 2)), max(vecnorm(c_una - c_mseg, 2, 2)), norm(c_opt))
end
% beta = 0.1 beta_opt: max||media - MSEG|| = 0.013   max||una - MSEG|| = 0.113   ||c_opt|| = 1.208
% beta = 0.5 beta_opt: max||media - MSEG|| = 0.074   max||una - MSEG|| = 0.311   ||c_opt|| = 1.208
% beta = 1.0 beta_opt: max||media - MSEG|| = 0.192   max||una - MSEG|| = 0.553   ||c_opt|| = 1.208

%% 6. Entrada blanca (canal ideal): donde valen las hipótesis de la filmina
for f = [0.5 1]
    mse  = lms(1, s2, N, 5, f/(N*(1 + s2)), 4000, 400);
    xb   = s2/(1 + s2);  piso = mean(mse(201:end));
    ex   = (mse - piso)/(mse(1) - piso);
    fprintf('blanca, beta = %.1f beta_opt: M medido = %.3f (predicho %.3f), tau medido = %d, -1/ln(gamma) = %.2f, 1/(2 beta phi0) = %.2f\n', ...
            f, piso/xb - 1, f/(2-f), find(ex < exp(-1), 1) - 1, -1/log(1 - 2*f/N + f^2/N), N/(2*f))
end
% blanca, beta = 0.5 beta_opt: M medido = 0.334 (predicho 0.333), tau medido = 15, -1/ln(gamma) = 14.16, 1/(2 beta phi0) = 11.00
% blanca, beta = 1.0 beta_opt: M medido = 0.998 (predicho 1.000), tau medido = 11, -1/ln(gamma) = 10.49, 1/(2 beta phi0) = 5.50

function [mse, c_media, e1, c, c_una, pico] = lms(p, s2, N, d, beta, R, K)
% LMS entrenado, R realizaciones en paralelo: cada fila es una corrida.
L  = numel(p);
A  = sign(randn(R, K+N+L));                          % símbolos PAM-2 conocidos
Rx = filter(p, 1, A, [], 2) + sqrt(s2)*randn(R, K+N+L);
c  = zeros(R, N);  mse = zeros(1, K);  c_media = zeros(K, N);
e1 = zeros(1, K);  c_una = zeros(K, N);  pico = zeros(R, 1);
for k = 1:K
    n = k + N + L - 1;
    r = Rx(:, n:-1:n-N+1);                           % r_k de cada realización
    E = A(:, n-d) - sum(r.*c, 2);                    % E_k = A_k - Y_k
    mse(k) = mean(E.^2);  c_media(k,:) = mean(c, 1);
    e1(k) = E(1);  c_una(k,:) = c(1,:);
    c = c + beta*E.*r;                               % c_{k+1} = c_k + beta E_k r_k^*
    pico = max(pico, max(abs(c), [], 2));            % el tap más grande que tuvo cada corrida
end
end

Detalles que muerden si no se los conoce:

• rand('seed',1) no fija randn. Es la primera línea de SG.m, y pone a rand en modo legado; randn sigue con su propio estado, así que dos randn después de dos rand('seed',1) dan números distintos y cada corrida del script da otro piso. Peor: en la misma sesión, un rng(1) posterior tira error, porque rng no acepta el generador legado. Hay que llamar antes a rng('default').
• El piso se mide en régimen, y el transitorio con el ensamble. El mean(mse(3001:end)) del bloque 3 descarta la primera mitad; con el paso más chico la curva llega a 1 dB del piso recién en la muestra 1517, así que una ventana que empiece antes mezcla transitorio con piso.
• Contar divergencias al final no alcanza. c al terminar sólo dice cómo quedó cada corrida; el bloque 4 guarda además pico, el tap más grande que tuvo cada una en todo el camino, porque cerca del borde las realizaciones se escapan y vuelven. Y una sola fila con Inf hace que mean(E.^2) dé Inf o NaN.
• SG.m sólo sirve con señales reales. sum(c.*buffer) y c - step*error*buffer no conjugan nada. Con QAM el filtro es c.'*buffer y la actualización lleva conj(buffer). Y la variable error tapa la función error de MATLAB mientras exista.

Para tener a mano

1. El LMS es el MSEG con la esperanza sacada: ${\nabla }_c|E_k{|}^2=-2E_k{r}_k^{\ast }$, y $c_{k+1}=c_k+\beta E_k{r}_k^{\ast }$, una vez por símbolo. Cuesta $2N+1$ multiplicaciones.
2. El gradiente de un símbolo no tiene sesgo, y con paso chico el promedio de las trayectorias es la trayectoria del MSEG. Cada realización, no.
3. No se queda quieto en $c_{\text{opt}}$: con entrada blanca, $E\Vert q_{k+1}{\Vert }^2=\gamma E\Vert q_k{\Vert }^2+{\beta }^{2}N{\phi }_0{\xi }_{\text{min}}$, con $\gamma =1-2\beta {\phi }_0+{\beta }^{2}N{\phi }_0^2$ (la filmina tiene $\beta$ donde va ${\beta }^2$).
4. ${\beta }_{\text{opt}}=1/(N{\phi }_0)$, estable si $0<\beta <2{\beta }_{\text{opt}}$ con entrada blanca; con el canal de la nota el MSE promedio ya se dispara en $\mathrm{1,3}{\beta }_{\text{opt}}$. Todo depende de $\beta {\phi }_0$: de ahí el AGC.
5. $\tau \approx 1/(2\beta {\phi }_0)$ con paso chico; con ${\beta }_{\text{opt}}$, del orden de $N$ (medido: 11 muestras con 11 taps), no $N/2$.
6. Desajuste $M=f/(2-f)$ con $f=\beta /{\beta }_{\text{opt}}$: con el óptimo el MSE se duplica. Exacta con entrada blanca, optimista con un canal correlacionado (1,85 medido contra 1).
7. Más paso, llega antes y se queda más lejos. La curva de aprendizaje se promedia sobre realizaciones: en régimen, el error de una corrida sola barre 22 dB entre sus percentiles 10 y 90.