El criterio MSE y el ecualizador de Wiener

Con ruido no alcanza con anular la ISI. El zero forcing la anula, y para hacerlo amplifica el ruido. El criterio MSE mide las dos cosas con un solo número, la potencia del error que llega al slicer, y elige los taps que lo hacen mínimo. Ese número es una forma cuadrática en los taps: tiene un único mínimo, y el mínimo sale de resolver un sistema lineal, las ecuaciones de Wiener-Hopf. En la cadena canal → AGC → ecualizador → slicer, el criterio dice adónde tiene que llegar el ecualizador. Cómo llega, paso a paso, es cosa del descenso por gradiente y del LMS.

01 · La definición

El ecualizador es un FIR transversal de $N$ coeficientes. La filmina los junta en un vector, y en otro junta las muestras que hay en la línea de retardos, futuras y pasadas:

$$c={[c_{-L},\dots ,c_L]}^T,r_k={[R_{k+L},\dots ,R_k,\dots ,R_{k-L}]}^T$$

con $N=2L+1$. Adentro de $r_k$ está el canal: $R_k$ son los símbolos $A_k$ pasados por el canal discreto equivalente $H(z)$, más el ruido $Z_k$. Todo se supone estacionario en sentido amplio. La salida del ecualizador, que es lo que entra al slicer, es

$$Y_k=c^T{r}_k$$

Si el slicer no se equivoca —durante el entrenamiento, o suponiendo que decide bien—, su salida ${\widehat{A}}_k$ es el símbolo transmitido $A_k$, y el error de la señal es

$$E_k=A_k-Y_k$$

El criterio MSE elige $c$ para que la potencia de ese error sea mínima:

$$\text{MSE}=E[|E_k{|}^2]$$

Mide junto lo que el ZF mira por separado. En $E_k$ hay ISI residual —lo que el ecualizador no terminó de cancelar de los símbolos vecinos— y hay ruido filtrado por $c$, y el MSE suma las dos potencias. La filmina del MMSE lineal lo dice así: desde el punto de vista de la probabilidad de error puede convenir dejar algo de ISI a la entrada del detector, si eso baja el ruido. En esa filmina el error va con el signo cambiado, $e_k=x_k-a_k$. Al cuadrado da lo mismo.

02 · El MSE es cuadrático en los taps

Cómo se obtiene: se desarrolla el módulo al cuadrado y se toma la esperanza.

Paso 1: desarrollar el cuadrado. Con $|x{|}^2=x{x}^{\ast }$ y ${(c^T{r}_k)}^{\ast }=c^H{r}_k^{\ast }$:

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

Los dos términos del medio son conjugados uno del otro, y un número más su conjugado es dos veces su parte real:

$$|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$$

Paso 2: tomar la esperanza. $c$ no es aleatorio: es lo que se elige. Sale afuera de la esperanza, y lo que queda adentro son dos estadísticas de las señales:

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

Paso 3: ponerles nombre.

$$\mathrm{\Phi }=E[r_k^{\ast }{r}_k^T],\alpha =E[A_k{r}_k^{\ast }]$$

$\mathrm{\Phi }$ es la matriz de autocorrelación de las muestras de entrada: su elemento $(i,j)$ dice cuánto se parecen, en promedio, las muestras que están en los taps $i$ y $j$. Como $R_k$ es estacionario, eso depende sólo de la distancia entre los taps: $\mathrm{\Phi }$ es Toeplitz, con ${\phi }_m=E[R_{k+m}{R}_k^{\ast }]$ a lo largo de cada diagonal, y es hermítica. $\alpha$ es la correlación cruzada entre el símbolo que se quiere recuperar y cada muestra de la línea: dice dónde está $A_k$ dentro de la ventana. Con los nombres puestos queda la forma de la filmina:

$$E[|E_k{|}^2]=E[|A_k{|}^2]-2\mathrm{\Re }(c^H\alpha )+c^H\mathrm{\Phi }c$$

Un término que no depende de $c$, uno lineal y uno cuadrático: una forma cuadrática en los taps. La filmina escribe $c^{\prime }$ donde acá dice $c^H$: es la misma traspuesta conjugada, con la comilla de MATLAB. Con símbolos reales, como en PAM-2, los asteriscos no hacen nada y la parte real tampoco.

03 · Un tazón con un solo fondo

Que sea cuadrática no alcanza: una cuadrática también puede ser una silla de montar, o una canaleta sin fondo. Lo que la vuelve un tazón es que $\mathrm{\Phi }$ sea definida positiva, y lo es por construcción:

$$c^H\mathrm{\Phi }c=E[|c^T{r}_k{|}^2]\ge 0$$

Es la potencia de la salida del ecualizador, que no puede ser negativa. Con ruido blanco a la entrada es, además, estrictamente positiva para todo $c\ne 0$: el ruido le suma ${\sigma }_Z^2$ a la diagonal de $\mathrm{\Phi }$, así que $c^H\mathrm{\Phi }c\ge {\sigma }_Z^2{c}^{H}c>0$.

Completar el cuadrado. Es el camino A de la filmina: se suma y se resta ${\alpha }^H{\mathrm{\Phi }}^{-1}\alpha$, y el MSE se reordena en

$$E[|E_k{|}^2]=E[|A_k{|}^2]-{\alpha }^H{\mathrm{\Phi }}^{-1}\alpha +{(c-{\mathrm{\Phi }}^{-1}\alpha )}^H\mathrm{\Phi }(c-{\mathrm{\Phi }}^{-1}\alpha )$$

Para verificarlo alcanza con desarrollar el último término: da $c^H\mathrm{\Phi }c-c^H\alpha -{\alpha }^{H}c+{\alpha }^H{\mathrm{\Phi }}^{-1}\alpha$, porque ${\mathrm{\Phi }}^H=\mathrm{\Phi }$, y el último sumando se cancela con el que se restó. Ahora se lee todo. Los dos primeros términos no dependen de $c$. El tercero es mayor o igual que cero porque $\mathrm{\Phi }$ es definida positiva, y vale cero sólo en $c={\mathrm{\Phi }}^{-1}\alpha$. Ahí está el mínimo, y es uno solo.

Las curvas de nivel son elipses centradas en ese punto, con los ejes sobre los autovectores de $\mathrm{\Phi }$. En la dirección del autovalor ${\lambda }_i$, subir el MSE en $\Delta$ exige alejarse $\sqrt{\Delta /{\lambda }_i}$: donde el autovalor es grande el tazón es empinado, donde es chico es casi plano. En el ejemplo de más abajo, a 10 dB, los autovalores son 10,10 y 0,626, y las elipses son 4,0 veces más largas que anchas. Qué le hace esa forma a un algoritmo que baja por el tazón está en el descenso por gradiente.

04 · Las ecuaciones de Wiener-Hopf

El camino B de la filmina llega al mismo punto sin completar cuadrados: en el fondo del tazón el gradiente es cero.

$${\nabla }_{c}E[|E_k{|}^2]=2\mathrm{\Phi }c-2\alpha =0$$

Son $N$ ecuaciones lineales con $N$ incógnitas, las ecuaciones de Wiener-Hopf. La filmina las llama solución MMSE:

$$\mathrm{\Phi }{c}_{\text{opt}}=\alpha \Rightarrow c_{\text{opt}}={\mathrm{\Phi }}^{-1}\alpha$$

El MSE que queda en el fondo sale de la forma completada, donde el término cuadrático se anula:

$${\xi }_{\text{min}}=E[|A_k{|}^2]-{\alpha }^H{\mathrm{\Phi }}^{-1}\alpha =E[|A_k{|}^2]-{\alpha }^H{c}_{\text{opt}}$$

$c_{\text{opt}}$ es el mejor ecualizador lineal de $N$ taps que existe para esa estadística, y ${\xi }_{\text{min}}$ es un piso: ningún algoritmo adaptativo baja de ahí. El LMS llega cerca de $c_{\text{opt}}$ sin conocer $\mathrm{\Phi }$ ni $\alpha$, y lo que le sobra se mide como un exceso por encima de ${\xi }_{\text{min}}$. Para que haya solución $\mathrm{\Phi }$ tiene que tener inversa, y la filmina avisa cuándo no la tiene: si el canal tiene un cero en frecuencia, las filas dejan de ser independientes. Con ruido blanco no llega a pasar, por el ${\sigma }_Z^2$ de la diagonal, pero cerca de un cero $\mathrm{\Phi }$ está mal condicionada y $c_{\text{opt}}$ tiene taps grandes.

05 · El principio de ortogonalidad

La misma ecuación, leída de otra manera. La correlación entre el error y las muestras de la línea, para un $c$ cualquiera, es

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

y en $c_{\text{opt}}$ vale cero: el error del ecualizador óptimo no se parece a ninguna de las muestras que el ecualizador tiene a mano. Es el principio de ortogonalidad, y vale en los dos sentidos: si el error está descorrelacionado con todas las entradas, $c$ es el óptimo, porque $\alpha -\mathrm{\Phi }c=0$ son las ecuaciones de Wiener-Hopf. La razón se ve sin cuentas: si al error le quedara algo parecido a $R_{k-m}$, se lo podría restar moviendo el tap que multiplica a $R_{k-m}$, y el error bajaría.

Salvo el factor $-2$, esa correlación es el gradiente del MSE. Por eso la filmina la presenta como la primera idea de cómo adaptar: medir cuánto se parece el error a cada entrada y mover los taps hasta que no se parezca a ninguna. Es lo que hace el LMS.

En el ejemplo de la sección siguiente, a 10 dB y sobre medio millón de muestras, la correlación medida en $c_{\text{opt}}$ da $[+\mathrm{0,0004},+\mathrm{0,0003}]$, que es lo que fluctúa un promedio de ese largo. En el ZF da $[-\mathrm{0,0996},+\mathrm{0,0903}]$, y la forma cerrada es ${\sigma }_Z^2[-1,\rho ]=[-\mathrm{0,1},+\mathrm{0,09}]$. Al error del ZF todavía le queda algo de las entradas: no es el óptimo.

06 · Un ejemplo que se hace a mano

Un canal con un solo polo, $H(z)=1/(1-\rho z^{-1})$: cada símbolo deja una cola geométrica, $h_n={\rho }^n$, y la ISI no termina nunca. Con $\rho =\mathrm{0,9}$, PAM-2 de $\pm 1$ (así $E[|A_k{|}^2]=1$), ruido blanco de potencia ${\sigma }_Z^2$, y un ecualizador de dos taps. Con $N=2$ no hay tap central, así que la línea es causal: $r_k={[R_k,R_{k-1}]}^T$. La SNR es la potencia de los símbolos sobre la del ruido, $1/{\sigma }_Z^2$.

Dos taps invierten este canal exacto. El ZF es $1-\rho z^{-1}$, que pone un cero justo sobre el polo del canal —la cuenta de polos y ceros—, así que $c_{\text{ZF}}={[1,-\rho ]}^T$ es un punto del plano y no una aproximación truncada. No deja ISI: todo su error es ruido, que pasa por $1-\rho z^{-1}$ y sale multiplicado por $1+{\rho }^2=\mathrm{1,81}$.

La autocorrelación de $R_k$ sale de sumar la serie geométrica, ${\phi }_m={\rho }^{|m|}/(1-{\rho }^2)+{\sigma }_Z^2{\delta }_m$, y $\alpha$ se lee del canal: $R_k$ contiene a $A_k$ con peso $h_0=1$, y $R_{k-1}$ sólo tiene símbolos anteriores, independientes de $A_k$:

$$\mathrm{\Phi }=\left[\begin{array}{cc}{\phi }_0& {\phi }_1\\ {\phi }_1& {\phi }_0\end{array}\right]=\left[\begin{array}{cc}\mathrm{5,263}+{\sigma }_Z^2& \mathrm{4,737}\\ \mathrm{4,737}& \mathrm{5,263}+{\sigma }_Z^2\end{array}\right],\alpha =\left[\begin{array}{c}1\\ 0\end{array}\right]$$

El sistema de 2 × 2 se resuelve a mano:

$$c_{\text{opt}}=\frac{1}{{\phi }_0^2-{\phi }_1^2}\left[\begin{array}{c}{\phi }_0\\ -{\phi }_1\end{array}\right],{\xi }_{\text{min}}=1-\frac{{\phi }_0}{{\phi }_0^2-{\phi }_1^2}$$

Sin ruido, ${\phi }_0^2-{\phi }_1^2$ vale $1/(1-{\rho }^2)={\phi }_0$, y $c_{\text{opt}}$ da ${[1,-\rho ]}^T$: el óptimo es el ZF, con ${\xi }_{\text{min}}=0$. A 10 dB, con ${\sigma }_Z^2=\mathrm{0,1}$, $c_{\text{opt}}$ vale ${[\mathrm{0,8478},-\mathrm{0,7488}]}^T$ y ${\xi }_{\text{min}}=\mathrm{0,1522}$ (−8,18 dB), contra ${\sigma }_Z^2(1+{\rho }^2)=\mathrm{0,1810}$ (−7,42 dB) del ZF. Simulado, sobre medio millón de muestras, el MSE medido da 0,1517 y 0,1805.

A 30 dB los dos puntos se tapan. Al bajar la SNR el mínimo se despega del ZF y viaja hacia el origen por la curva punteada: Wiener achica los dos taps. El tazón entero sube, porque su fondo ${\xi }_{\text{min}}$ sube con el ruido, y se ensancha. Y hay una cosa que el dibujo no muestra y las lecturas sí: si multiplicás $c$ por un número, con los botones o arrastrándolo sobre una recta que pase por el origen, el MSE cambia y la BER no se mueve. Esa es la sección 08.

07 · MMSE contra ZF, según la SNR

Para ver qué negocia Wiener hay que partir su error. El pulso total, canal por ecualizador, es $p=h\ast c$; con este canal vale $p_0=c_0$ y $p_n={\rho }^{n-1}(\rho c_0+c_1)$ para $n\ge 1$. Con símbolos independientes y ruido independiente de los símbolos, el MSE se parte en tres potencias:

$$E[|E_k{|}^2]=\underset{\text{cursor}}{\underbrace{{(1-p_0)}^2}}+\underset{\text{ISI}}{\underbrace{\sum _{n\ge 1}{p}_n^2}}+\underset{\text{ruido}}{\underbrace{{\sigma }_Z^2(c_0^2+c_1^2)}}$$

El cursor que no llega a 1, los vecinos que no se cancelaron y el ruido que pasó por los taps. El ZF tiene las dos primeras en cero por construcción, y todo su MSE es ruido, ${\sigma }_Z^2(1+{\rho }^2)$. Las tres de Wiener, calculadas en MATLAB:

A 30 dB son el mismo ecualizador: los taps difieren en la tercera cifra decimal, y el MSE, en 0,008 dB. A 0 dB son dos cosas distintas. El ZF amplifica el ruido hasta 1,81, más que la potencia de los símbolos: al slicer le llega más error que señal. Wiener no lo hace. Achica los dos taps, acepta un cursor de 0,373 en vez de 1, deja pasar 0,0151 de ISI, y a cambio el ruido le queda en 0,2188. Su MSE es 0,627, 4,6 dB menos que el del ZF; a 10 dB la diferencia es de 0,75 dB, y a 30 dB, de 0,008 dB. Con este canal lo que más cede no es ISI sino cursor: lo barato es achicar la señal entera, con el ruido adentro.

Es el dibujo de la filmina del MMSE lineal, donde el canal se llama $M(z)$. Tres curvas del ecualizador en frecuencia: con SNR → ∞ tiende a $1/M(z)$, el ZF; con SNR → 0 deja de ecualizar; en el medio levanta lo que el canal atenuó, pero menos. Y la misma filmina avisa lo que el ZF no puede evitar: el ruido sale amplificado, y ya no es blanco.

08 · Mínimo MSE no es mínima BER

La filmina lo advierte: minimizar el error total no es lo mismo que minimizar la probabilidad de error, porque la ISI residual no tiene una distribución gaussiana. Es una aproximación razonable. El ejemplo muestra cuánto se separan las dos cosas. A 0 dB, con la BER calculada exacta —enumerando los vecinos que pesan— y medida sobre un millón de símbolos:

Wiener y el ZF están a 4,6 dB en MSE y a 0,008 en BER. Y $3c_{\text{opt}}$ da vuelta el orden: tiene más MSE que el ZF y la misma BER que Wiener. Su MSE sale de la forma cuadrática, $1-2a{c}_0+a^2{c}_0$ para $a{c}_{\text{opt}}$, porque $\mathrm{\Phi }{c}_{\text{opt}}=\alpha$. Son dos razones.

La escala. En PAM-2 el umbral está en cero, y multiplicar $c$ por un número positivo no le cambia el signo a ninguna $Y_k$: no cambia ninguna decisión. El MSE sí cambia, porque mide la distancia a $A_k$, no de qué lado del umbral cayó $Y_k$.

La forma del error. Si el error fuera gaussiano, la BER se leería del MSE: $Q(1/\sqrt{{\xi }_{\text{min}}})=\mathrm{0,1033}$. La verdadera es el doble. Lo que falla es el cursor: la distancia al umbral no es 1 sino $p_0=\mathrm{0,373}$, y con esa distancia, $Q(p_0/\sqrt{\text{ISI}+\text{ruido}})=\mathrm{0,2202}$ ya pega. El MSE junta en un número el sesgo y la dispersión, y a la BER le importa la dispersión comparada con el cursor.

En este canal, igual, Wiener apunta bien. Como la BER de PAM-2 depende sólo de la dirección de $c$, se puede buscar la dirección de mínima BER con fminbnd: a 0 dB está a 4·10⁻⁴ rad de la de Wiener, y la BER difiere en la séptima cifra. Lo que el MSE mide mal no es hacia dónde apuntar, sino cuánto vale lo que se logró. Con 16-QAM, cuyos umbrales no están en cero, el cursor corto sí mueve decisiones.

09 · El análisis en MATLAB

El ejemplo entero: $\mathrm{\Phi }$ y $\alpha$ analíticas y estimadas, Wiener-Hopf con la barra invertida, el ZF, la ortogonalidad, el barrido de SNR y la BER. Lo que figura como comentario es la salida real de correrlo en R2026a.

%% 1. El canal, los símbolos y el ruido
rho = 0.9;                          % canal H(z) = 1/(1 - rho z^-1): h_n = rho^n
SNRdB = 10;                         % E|A_k|^2 / sigma_Z^2, con PAM-2 de +-1
s2 = 10^(-SNRdB/10);                % potencia del ruido

%% 2. Phi y alpha, analíticas (dos taps causales: r_k = [R_k; R_{k-1}])
phi = @(m) rho.^abs(m) / (1 - rho^2);        % autocorrelación de la señal sin ruido
Phi = toeplitz([phi(0) phi(1)]) + s2*eye(2);
alpha = [1; 0];                              % E[A_k R_k] = h_0, E[A_k R_{k-1}] = 0
fprintf('Phi = [%.4f %.4f; %.4f %.4f], autovalores %.4f y %.4f\n', Phi, eig(Phi))
% Phi = [5.3632 4.7368; 4.7368 5.3632], autovalores 0.6263 y 10.1000

%% 3. Wiener-Hopf, con la barra invertida
c_opt  = Phi \ alpha;
xi_min = 1 - alpha' * c_opt;
fprintf('c_opt = [%.4f %.4f], xi_min = %.4f (%.2f dB)\n', c_opt, xi_min, 10*log10(xi_min))
% c_opt = [0.8478 -0.7488], xi_min = 0.1522 (-8.18 dB)

%% 4. El ZF: con este canal, dos taps lo invierten exacto
c_zf = [1; -rho];
mse  = @(c) 1 - 2*alpha'*c + c'*Phi*c;       % la forma cuadrática
fprintf('MSE ZF = %.4f (%.2f dB)\n', mse(c_zf), 10*log10(mse(c_zf)))
% MSE ZF = 0.1810 (-7.42 dB)

%% 5. Phi y alpha estimadas con muestras
N = 1e6;  rng(1)
A = 2*randi([0 1], N, 1) - 1;
R = filter(1, [1 -rho], A) + sqrt(s2)*randn(N, 1);
k = (101:N)';                                % fuera el transitorio del canal
r = [R(k) R(k-1)];                           % cada fila es r_k.'
ent = 1:floor(numel(k)/2);  pru = ent(end)+1:numel(k);
Phi_est   = r(ent,:)' * r(ent,:) / numel(ent);
alpha_est = r(ent,:)' * A(k(ent)) / numel(ent);
c_est = Phi_est \ alpha_est;
fprintf('c_est = [%.4f %.4f]\n', c_est)
% c_est = [0.8472 -0.7483]

%% 6. El MSE medido, sobre muestras que no se usaron para estimar
Ek = @(c) A(k(pru)) - r(pru,:)*c;            % E_k = A_k - Y_k
fprintf('MSE medido: Wiener %.4f, estimado %.4f, ZF %.4f\n', ...
    mean(Ek(c_opt).^2), mean(Ek(c_est).^2), mean(Ek(c_zf).^2))
% MSE medido: Wiener 0.1517, estimado 0.1517, ZF 0.1805

%% 7. El principio de ortogonalidad: E[E_k r_k] = alpha - Phi*c
fprintf('E[E r] Wiener: [%+.4f %+.4f]   cerrado [%+.4f %+.4f]\n', mean(Ek(c_opt).*r(pru,:)), alpha - Phi*c_opt)
% E[E r] Wiener: [+0.0004 +0.0003]   cerrado [+0.0000 +0.0000]
fprintf('E[E r] ZF:     [%+.4f %+.4f]   cerrado [%+.4f %+.4f]\n', mean(Ek(c_zf).*r(pru,:)), alpha - Phi*c_zf)
% E[E r] ZF:     [-0.0996 +0.0903]   cerrado [-0.1000 +0.0900]
Ee = A(k(ent)) - r(ent,:)*c_est;             % el estimado, sobre sus propias muestras
fprintf('E[E r] estimado, mismas muestras: [%.1e %.1e]\n', mean(Ee.*r(ent,:)))
% E[E r] estimado, mismas muestras: [6.7e-15 -5.8e-15]

%% 8. El barrido de SNR: cursor, ISI, ruido y BER
fprintf('\n SNR   c0      c1      cursor  ISI     ruido   MSE     MSE_ZF  BER       BER_ZF\n')
for SNRdB = [30 20 10 5 0]
    s2  = 10^(-SNRdB/10);
    Phi = toeplitz([phi(0) phi(1)]) + s2*eye(2);
    c   = Phi \ alpha;
    p0  = c(1);                              % pulso total p = h*c: p_0 = c_0
    kp  = rho*c(1) + c(2);                   % p_n = rho^(n-1)*kp para n >= 1
    isi = kp^2 / (1 - rho^2);
    ruido = s2 * (c'*c);
    fprintf('%3d  %6.4f %7.4f  %6.4f  %6.4f  %6.4f  %6.4f  %6.4f  %8.2e  %8.2e\n', SNRdB, c, ...
        (1-p0)^2, isi, ruido, (1-p0)^2 + isi + ruido, s2*(1+rho^2), ...
        ber_pam2(p0, kp, ruido, rho), qf(1/sqrt(s2*(1+rho^2))))
end
%  SNR   c0      c1      cursor  ISI     ruido   MSE     MSE_ZF  BER       BER_ZF
%  30  0.9982 -0.8982  0.0000  0.0000  0.0018  0.0018  0.0018  1.77e-122  1.81e-122
%  20  0.9822 -0.8823  0.0003  0.0000  0.0174  0.0178  0.0181  5.19e-14  5.31e-14
%  10  0.8478 -0.7488  0.0232  0.0011  0.1280  0.1522  0.1810  9.13e-03  9.37e-03
%   5  0.6419 -0.5450  0.1282  0.0056  0.2242  0.3581  0.5724  9.03e-02  9.31e-02
%   0  0.3730 -0.2821  0.3931  0.0151  0.2188  0.6270  1.8100  2.20e-01  2.29e-01

%% 9. La BER medida, con 1e6 símbolos, a 0 dB
s2 = 1;  Phi = toeplitz([phi(0) phi(1)]) + s2*eye(2);  c_opt = Phi \ alpha;
R = filter(1, [1 -rho], A) + sqrt(s2)*randn(N, 1);  r = [R(k) R(k-1)];
fprintf('\nBER medida a 0 dB: Wiener %.4f, ZF %.4f, 3*c_opt %.4f\n', ...
    mean(sign(r*c_opt) ~= A(k)), mean(sign(r*c_zf) ~= A(k)), mean(sign(r*3*c_opt) ~= A(k)))
% BER medida a 0 dB: Wiener 0.2206, ZF 0.2292, 3*c_opt 0.2206
fprintf('MSE medido a 0 dB: Wiener %.4f, ZF %.4f, 3*c_opt %.4f\n', ...
    mean((A(k) - r*c_opt).^2), mean((A(k) - r*c_zf).^2), mean((A(k) - r*3*c_opt).^2))
% MSE medido a 0 dB: Wiener 0.6277, ZF 1.8112, 3*c_opt 2.1212
fprintf('Si el error fuera gaussiano: Q(1/sqrt(xi_min)) = %.4f\n', qf(1/sqrt(1 - alpha'*c_opt)))
% Si el error fuera gaussiano: Q(1/sqrt(xi_min)) = 0.1033

%% Funciones
function p = qf(x)
    p = 0.5*erfc(x/sqrt(2));
end

function p = ber_pam2(p0, kp, ruido, rho)
    % BER exacta de PAM-2 con ISI p_n = rho^(n-1)*kp: se enumeran los 18
    % primeros vecinos; la cola, de varianza despreciable, se suma al ruido.
    K = 18;
    s = 1 - 2*mod(floor((0:2^K-1)' ./ 2.^(0:K-1)), 2);       % todas las secuencias de +-1
    isi = s * (kp * rho.^(0:K-1)');
    cola = kp^2 * rho^(2*K) / (1 - rho^2);
    p = mean(qf((p0 + isi) / sqrt(ruido + cola)));
end

Los detalles que muerden:

• Phi \ alpha, no inv(Phi)*alpha. La barra resuelve el sistema sin invertir la matriz: hace menos cuentas y pierde menos precisión cuando $\mathrm{\Phi }$ está mal condicionada.
• La comilla de MATLAB conjuga. Con una fila $r_k^T$ por muestra, r'*r suma $r_k^{\ast }{r}_k^T$, que es justo $\mathrm{\Phi }$, y r'*A suma $A_k{r}_k^{\ast }$, que es $\alpha$. Con QAM eso importa; trasponer sin conjugar es .'.
• filter arranca con el canal descargado. En las primeras muestras la cola de ${\rho }^n$ todavía no se armó, así que no son estacionarias y sesgan la estimación. El bloque 5 tira las cien primeras.
• La ortogonalidad no se chequea con las mismas muestras. Si $c$ se estima con un bloque y la correlación se mide sobre ese mismo bloque, da cero —6,7·10⁻¹⁵— por construcción: es la ecuación que se acaba de resolver, no una verificación. El bloque 7 mide sobre la otra mitad.
• La estimación mejora como $1/\sqrt{N}$. Con 200 repeticiones a 10 dB, la distancia de $c$ estimado a $c_{\text{opt}}$, en valor cuadrático medio, es 0,049 con 100 muestras, 0,015 con mil, 0,0048 con diez mil y 0,0014 con cien mil. Cada década de muestras la divide por algo más de tres.

Para tener a mano

1. $E_k=A_k-Y_k$, con $Y_k=c^T{r}_k$. El criterio minimiza $E[|E_k{|}^2]$, que suma la potencia de la ISI residual y la del ruido.
2. Es una forma cuadrática en los taps:$$E[|E_k{|}^2]=E[|A_k{|}^2]-2\mathrm{\Re }(c^H\alpha )+c^H\mathrm{\Phi }c$$con $\mathrm{\Phi }=E[r_k^{\ast }{r}_k^T]$ y $\alpha =E[A_k{r}_k^{\ast }]$.
3. $\mathrm{\Phi }$ es Toeplitz, hermítica y definida positiva: el MSE es un tazón con un solo fondo, de curvas de nivel elípticas con los ejes sobre los autovectores de $\mathrm{\Phi }$.
4. Wiener-Hopf: $\mathrm{\Phi }{c}_{\text{opt}}=\alpha$, y ${\xi }_{\text{min}}=E[|A_k{|}^2]-{\alpha }^H{c}_{\text{opt}}$. En MATLAB, c_opt = Phi \ alpha.
5. Ortogonalidad: $E[E_k{r}_k^{\ast }]=\alpha -\mathrm{\Phi }c$, que se anula sólo en $c_{\text{opt}}$. Es el gradiente, salvo un $-2$.
6. Con SNR alta Wiener es el ZF. Con SNR baja achica los taps, deja cursor corto y algo de ISI, y no amplifica el ruido: a 0 dB, 0,627 de MSE contra 1,810.
7. Mínimo MSE no es mínima BER: la BER de PAM-2 no ve la escala de $c$, y el error no es gaussiano. A 0 dB, 4,6 dB de MSE son 0,2203 contra 0,2287 de BER.