El PLL en tiempo discreto

El PLL de la recuperación fina de portadora no tiene VCO ni capacitores: es aritmética que corre una vez por símbolo. Los dos integradores del lazo de la filmina pasan a ser dos acumuladores, Laplace pasa a ser la transformada Z, y aparece un bloque que en continuo no existe: la latencia. El script de la clase, phase_system.m, arma ese lazo en seis líneas. Hecha a mano, la cuenta da una transferencia de segundo orden con polos que se escriben en una línea, y resulta ser exactamente el lazo continuo pasado por Euler hacia atrás. Por eso valen el mismo $\zeta$ y el mismo ${\omega }_n$ de la filmina, mientras el lazo sea lento frente a la tasa de símbolo.

01 · La definición

Es el modelo equivalente de fase de la filmina, muestreado. Entra la fase $\theta [n]$ de la portadora recibida y sale la fase $\phi [n]$ del NCO, una muestra por símbolo. Linealizado alrededor del enganche, el lazo son tres ecuaciones:

$$\begin{array}{l}e[n]=\theta [n]-\phi [n-\mathrm{Lat}]\\ c[n]=K_{p}e[n]+K_i\sum _{k\le n}e[k]\\ \phi [n]=\phi [n-1]+c[n]\end{array}$$

La primera es el detector de fase: compara la entrada con la fase del NCO, que le llega Lat muestras tarde. La segunda es el filtro de lazo PI de la filmina, $L(s)=K_p+K_i/s$, con la integral cambiada por una suma: la rama proporcional mira el error de ahora y la integral, todos los errores hasta ahora. La tercera es el NCO, que acumula el control igual que el VCO lo integraba. En diagrama, θ[n] → (+) → L(z) → NCO → φ[n], con el retorno por z^−Lat.

En continuo el detector tenía ganancia $K_{\mathrm{PD}}$ y el VCO, $K_V$. Acá las dos están normalizadas, $K_{\mathrm{PD}}{K}_V=1$ por muestra, y lo que quede de ellas vive dentro de $K_p$ y $K_i$. Es lo que permite usar las fórmulas de la filmina con los números del script sin arrastrar unidades.

02 · El script, pieza por pieza

• fs = 16e9, con el comentario «fs = BR para el FCR». El lazo de recuperación fina de portadora corre una vez por símbolo, así que su frecuencia de muestreo es la tasa de símbolo, 16 GBd, y Ts son 62,5 ps. Ts sólo sirve para que los ejes salgan en segundos y en hertz: el lazo cuenta muestras.
• z = tf('z', Ts) crea $z$ como una transferencia discreta de período Ts. Desde ahí, toda expresión con z es otra tf, y MATLAB hace el álgebra de polinomios. Todas tienen que compartir el mismo Ts; mezclar períodos da error.
• L = Kp + Ki*z/(z-1) es el PI, con Kp = 0.005 y Ki = Kp/100. El cociente $z/(z-1)$ es el acumulador de la sección siguiente.
• NCO = z/(z-1) es otro acumulador, el mismo.
• G = L*NCO es el lazo abierto, de $e$ a $\phi$. Con dos acumuladores en serie su denominador es ${(z-1)}^2$: es un lazo de tipo 2.
• H = z^-Lat, con Lat = 0, es la latencia. Va en la realimentación; con cero queda H = 1.
• F = feedback(G, H) cierra el lazo: calcula $G/(1+GH)$, con realimentación negativa por omisión. Es la $F=\Phi /\Theta$ de la filmina, en $z$.

Y las cuatro figuras:

• step(F) es un escalón de fase de 1 rad: la portadora recibida salta de golpe. La fase del NCO lo alcanza con un sobrepico de 39,72 % y entra en la banda del 2 % a las 1477 muestras, 92,3 ns.
• step(F*z/(z-1)) junto a step(z/(z-1)) es la rampa. No hay una función de rampa para una tf, y el truco es que un escalón que pasa por un acumulador es una rampa: $\theta [n]=n+1$, un radián por muestra. Una rampa de fase es un corrimiento de frecuencia, y las dos curvas terminan una encima de la otra: el PI lo sigue sin error.
• pzmap(F) muestra dos polos complejos en 0,997463 ± 0,006581j, adentro del círculo y pegados a $z=1$, y dos ceros: uno en el origen y otro en $K_p/(K_p+K_i)=\mathrm{0,990099}$.
• bode(F) es el pasabajos de fase de la filmina, dibujado hasta $f_s/2$ = 8 GHz: un pico de 4,97 dB en 16,35 MHz, cerca de $f_n$, y el corte a −3 dB en 30,34 MHz, contra los 30,40 MHz de la fórmula continua de la filmina. Lo que pasa arriba de eso no lo copia el NCO: queda como error.

Dos restos del script: clear all borra todo y close all está comentado, así que las figuras se acumulan entre corridas; y hold all es el nombre viejo de hold on.

03 · Los dos acumuladores

La rama integral. Llamando $i[n]$ a la suma de la segunda ecuación, cada muestra le suma el error nuevo a la anterior:

$$i[n]=i[n-1]+K_{i}e[n]⟹I(z)(1-z^{-1})=K_{i}E(z)⟹\frac{I}{E}=\frac{K_i}{1-z^{-1}}=K_i\frac{z}{z-1}$$

Es la propiedad de retardo de la transformada Z, $x[n-1]\leftrightarrow z^{-1}X(z)$, y nada más. La suma incluye la muestra actual, así que el acumulador no retarda: el error de hoy ya está en la salida de hoy.

El NCO. La tercera ecuación es la misma cuenta con $c[n]$ en lugar de $K_{i}e[n]$: $\Phi /C=z/(z-1)$. Con un control constante, la fase crece una cantidad fija por muestra: el NCO cambió de frecuencia, igual que el VCO de la filmina con tensión constante. Un NCO que sumara sólo hasta la muestra anterior tendría un $z^{-1}$ más; en el script, eso es Lat = 1.

El polo en $z=1$. $z/(z-1)$ vale infinito en $z=e^{j0}=1$, que es frecuencia cero. Es el $L(0)\to \infty$ que la tabla de filtros de la filmina pide para seguir una rampa de fase sin error, escrito en $z$.

04 · La transferencia, paso a paso

Con Lat = 0, que es como viene el script.

Paso 1: el PI con un solo denominador. Se multiplica $K_p$ por $(z-1)/(z-1)$ y se suma:

$$L(z)=K_p+K_i\frac{z}{z-1}=\frac{(K_p+K_i)z-K_p}{z-1}$$

Paso 2: el lazo abierto. El NCO multiplica por otro $z/(z-1)$:

$$G(z)=L(z)\frac{z}{z-1}=\frac{z[(K_p+K_i)z-K_p]}{{(z-1)}^2}=\frac{N}{D}$$

Paso 3: cerrar el lazo. Con $H=1$, se multiplica arriba y abajo por $D$:

$$F=\frac{G}{1+G}=\frac{N/D}{1+N/D}=\frac{N}{D+N}$$

Paso 4: desarrollar el denominador. $D+N=z^2-2z+1+(K_p+K_i)z^2-K_{p}z$, y agrupando por potencias de $z$:

$$F(z)=\frac{\Phi (z)}{\Theta (z)}=\frac{z[(K_p+K_i)z-K_p]}{(1+K_p+K_i)z^2-(2+K_p)z+1}$$

Con los números de la clase, $(\mathrm{0,00505}{z}^2-\mathrm{0,005}z)/(\mathrm{1,00505}{z}^2-\mathrm{2,005}z+1)$, que es lo que devuelve tfdata(F, 'v'). Contra feedback coinciden exactamente: la diferencia es 0 en los coeficientes normalizados y 0 en la respuesta en frecuencia, sobre 4000 puntos de 0 a $\pi$. feedback llega a los mismos coeficientes por su propio camino, sumando los polinomios de $G$.

Paso 5: la comprobación en $z=1$. El numerador da $K_p+K_i-K_p=K_i$, y el denominador, $1+K_p+K_i-2-K_p+1=K_i$. $F(1)=1$: una fase constante se copia entera.

Paso 6: el error. La transferencia del error es lo que le falta a $F$:

$$\frac{E(z)}{\Theta (z)}=1-F=\frac{D}{D+N}=\frac{{(z-1)}^2}{(1+K_p+K_i)z^2-(2+K_p)z+1}$$

Doble cero en $z=1$, uno por acumulador. Con el teorema del valor final discreto, $e_{\mathrm{ss}}={lim}_{z\to 1}(z-1)E(z)$, el escalón $\Theta =z/(z-1)$ se lleva un cero y sobra otro, y la rampa $\Theta ={\omega }_{0}z/{(z-1)}^2$ se lleva los dos y deja una $(z-1)$ suelta: los dos errores son cero. Medido con lsim sobre una rampa de ${10}^{-3}$ rad por muestra —a 16 GBd es un corrimiento de 2,546 MHz—, el error a las 20 000 muestras es de 8·10⁻¹² rad: redondeo sobre una fase que ya vale 20 rad. Sin la rama integral el error es $(z-1)/((1+K_p)z-1)$, la rampa deja ${\omega }_0/K_p$ —la misma fila de la tabla de la filmina— y la medición da 0,2000 rad, 11,5°.

05 · Los polos

Son las raíces del denominador de $F$. El discriminante es ${(2+K_p)}^2-4(1+K_p+K_i)=K_p^2-4K_i$, porque el 4 y el $4K_p$ se cancelan. Queda:

$$z_{1,2}=\frac{(2+K_p)\pm j\sqrt{4K_i-K_p^2}}{2(1+K_p+K_i)},|z|=\frac{1}{\sqrt{1+K_p+K_i}}$$

El módulo sale del producto de las raíces, que es el coeficiente independiente sobre el principal: $1/(1+K_p+K_i)$. Tres lecturas:

• Son complejos si y sólo si $K_p^2<4K_i$, o sea $K_p/(2\sqrt{K_i})<1$. Es exactamente el $\zeta <1$ de la filmina, sin aproximar nada: el lazo discreto oscila cuando el continuo oscila.
• Con Lat = 0 siempre son estables: para cualquier $K_p,K_i>0$ el módulo es menor que 1. Ninguna ganancia saca este lazo del círculo. La latencia sí (sección 07).
• Con los números de la clase quedan en 0,997463 ± 0,006581j, de módulo 0,997485, a 0,00705 de $z=1$. Polos tan pegados a 1 quieren decir que en una muestra casi no pasa nada: el lazo tarda cientos de muestras. Contra pole(F), la fórmula difiere en 1,1·10⁻¹⁴.

06 · De s a z

Con $K_{\mathrm{PD}}{K}_V=1$ y el tiempo contado en muestras, la $F(s)$ de la filmina con el PI es

$$F(s)=\frac{K_{p}s+K_i}{s^2+K_{p}s+K_i}=\frac{2\zeta {\omega }_{n}s+{\omega }_n^2}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2},{\omega }_n=\sqrt{K_i},\zeta =\frac{K_p}{2\sqrt{K_i}}$$

Con los números de la clase, $\zeta =\mathrm{0,3536}$ y ${\omega }_n=\mathrm{0,00707}$ rad por muestra; multiplicando por $f_s/2\pi$, $f_n=\mathrm{18,01}$ MHz. El ancho de banda equivalente de ruido de la filmina, $B_n=\frac{{\omega }_n}{2}(\zeta +\frac{1}{4\zeta })$, da $B_n{T}_s=\mathrm{3,75}·{10}^{-3}$: 60 MHz contra 16 GBd. Ese producto es el que decide todo lo que sigue.

Paso 1: reemplazar cada integrador por un acumulador. El acumulador de la sección 03 es $1/(1-z^{-1})$ donde había $1/s$. O sea, $s\to 1-z^{-1}=(z-1)/z$: la derivada aproximada por la diferencia hacia atrás, que es Euler hacia atrás con $T_s=1$.

Paso 2: sustituir en $F(s)$. El numerador queda $K_p(z-1)/z+K_i=[(K_p+K_i)z-K_p]/z$, y el denominador,

$$\frac{{(z-1)}^2}{z^2}+K_p\frac{z-1}{z}+K_i=\frac{(1+K_p+K_i)z^2-(2+K_p)z+1}{z^2}$$

El cociente es, letra por letra, la $F(z)$ de la sección 04. El lazo del script es el lazo continuo de la filmina discretizado por Euler hacia atrás, no una aproximación de él.

Paso 3: dónde caen los polos. Si $s=1-z^{-1}$, entonces $z=1/(1-s)$, y cada polo continuo $s_c=-\zeta {\omega }_n\pm j{\omega }_n\sqrt{1-{\zeta }^2}$ va a parar a $1/(1-s_c)$; contra pole(F) la diferencia es 1,1·10⁻¹⁴. Esa función manda todo el semiplano izquierdo adentro del círculo unidad, y por eso el lazo con Lat = 0 no se desestabiliza con ninguna ganancia. Un muestreo ideal del lazo continuo, en cambio, pondría el polo en $e^{s_c}$. Las dos funciones coinciden al primer orden y se separan al segundo:

$$\frac{1}{1-s}=1+s+s^2+\dots e^s=1+s+\frac{s^2}{2}+\dots$$

La diferencia es del orden de $|s_c{|}^2/2={\omega }_n^2/2$, y ${\omega }_n$ en radianes por muestra es proporcional a $B_n{T}_s$. Con la clase son 2,5·10⁻⁵ entre $e^{s_c}$ y el polo de pole(F): invisible. Es el sentido exacto de «el lazo es lento frente a la tasa de símbolo».

Cuándo deja de valer. Con $\zeta$ fijo en el 0,3536 de la clase ($K_i=2K_p^2$) y $K_p$ subiendo, el modelo continuo tiene siempre el mismo sobrepico, 40,08 %, porque con el tiempo medido en $1/{\omega }_n$ la forma de la curva depende sólo de $\zeta$. El discreto, medido con stepinfo:

El sobrepico continuo sale de la forma cerrada de la tabla de la filmina; stepinfo, que la muestrea, da 40,07 %. La diferencia crece en proporción a $B_n{T}_s$: −0,36 puntos con la clase, el doble con el doble. Hasta $B_n{T}_s=\mathrm{0,015}$ el modelo continuo se equivoca en menos de 1,5 puntos; con 0,075 ya son 6,5, y con 0,375 el discreto tiene menos de la mitad del sobrepico que predice $\zeta$. El sentido es siempre el mismo: Euler hacia atrás amortigua. Los polos discretos, leídos de vuelta como $s=\ln z$, tienen un $\zeta$ más alto que el diseñado: 0,357 con la clase, 0,412 con 0,075, 0,692 con 0,75.

07 · La latencia

Con H = z^-Lat en la realimentación, el detector compara $\theta [n]$ con $\phi [n-\mathrm{Lat}]$: corrige con una fase vieja. En continuo no hay nada parecido, y en el hardware es inevitable por dos razones.

La primera: Lat = 0 no se puede construir. Con Lat = 0, $\phi [n]$ depende de $e[n]$ y $e[n]$ depende de $\phi [n]$. Es un lazo algebraico: MATLAB lo resuelve despejando,

$$\phi [n]=\frac{\phi [n-1]+i[n-1]+(K_p+K_i)\theta [n]}{1+K_p+K_i}$$

y por eso ante el escalón la primera muestra ya se mueve: $\phi [0]=(K_p+K_i)/(1+K_p+K_i)=\mathrm{0,005025}$, el mismo valor que da step. Un circuito no despeja ecuaciones: entre medir el error y mover el NCO hay por lo menos un registro, que es Lat = 1.

La segunda: a 16 GBd nada corre a 16 GHz. El receptor procesa muchos símbolos en paralelo en cada ciclo de reloj, y cada etapa del pipeline —el detector, el filtro, el NCO, la rotación de la constelación— agrega ciclos. Contada en símbolos, que es lo que cuenta Lat, esa demora se multiplica por el paralelismo.

Lo que hace con los polos. Con $H=z^{-\mathrm{Lat}}$, $F=G/(1+G{z}^{-\mathrm{Lat}})$, y multiplicando por $z^{\mathrm{Lat}}{(z-1)}^2$ el denominador pasa a ser

$$z^{\mathrm{Lat}}{(z-1)}^2+(K_p+K_i)z^2-K_{p}z$$

de grado $\mathrm{Lat}+2$. Con $\mathrm{Lat}\ge 1$ tiene una raíz en el origen que se cancela con un cero del numerador, así que el lazo tiene $\mathrm{Lat}+1$ polos: el par de siempre, cerca de $z=1$, y $\mathrm{Lat}-1$ nuevos. Con Lat = 1 el par casi no se mueve: 0,997475 ± 0,006605j. Con Lat = 50 está en 0,998590 ± 0,007756j, y los 49 nuevos forman un anillo de módulos entre 0,886 y 0,950, lejos de $z=1$. El par de siempre se va corriendo hacia afuera, y el sobrepico, medido con stepinfo, sube con él:

Con 100 símbolos de latencia el par de siempre sale del círculo y el lazo de la clase, que con Lat = 0 no se desestabilizaba con ninguna ganancia, es inestable. El techo no lo pone la ganancia: lo pone el retardo.

Con los valores de la clase el discreto y el continuo son la misma curva, y los polos azules tapan a los rombos grises. Subí $K_p$: la ventana de la derecha se agranda, los polos se alejan de $z=1$ y se despegan de sus $e^s$, y el escalón discreto se queda por debajo del punteado. Con $K_i$ fijo, subir $K_p$ también sube $\zeta$, y las dos curvas cambian de forma; lo que mide la diferencia entre ellas es $B_n{T}_s$. Después volvé a la clase y subí la latencia: aparece el anillo de polos nuevos, el par de siempre se arrima al círculo, el escalón oscila cada vez más, y antes de los 100 símbolos sale. En toda esa excursión el modelo continuo no se mueve: no tiene cómo ver la latencia.

08 · El análisis en MATLAB

El script de la clase tal cual, y abajo la verificación de lo derivado a mano. Lo que figura como comentario es la salida real de correrlo en R2026a.

clear all
%close all

fs = 16e9; %GHz (fs = BR para el FCR)
Ts = 1/fs;
z = tf('z', Ts);

%%
% System modelling
Kp = 0.005;
Ki = Kp/100;

L = Kp + Ki*z/(z-1);
NCO = z/(z-1);
G = L*NCO;
Lat = 0;
H = z^-Lat;
F = feedback(G, H); %Funcion del PLL

%% Step response
figure
step(F)

%% Ramp response
figure
step(F*z/(z-1))
hold all
step(z/(z-1))

% %%
figure
pzmap(F)

%%
figure
bode(F)
hold all

%% Verificación 1: la transferencia derivada a mano contra feedback
Fm = tf([Kp+Ki, -Kp, 0], [1+Kp+Ki, -(2+Kp), 1], Ts);
[nF, dF] = tfdata(F, 'v');
[nM, dM] = tfdata(Fm, 'v');
fprintf('coeficientes: %.1e\n', max(abs([nF/dF(1) - nM/dM(1), dF/dF(1) - dM/dM(1)])))
% coeficientes: 0.0e+00
w = linspace(1e-4, pi, 4000) / Ts;
fprintf('respuesta en frecuencia: %.1e\n', ...
        max(abs(squeeze(freqresp(F, w)) - squeeze(freqresp(Fm, w)))))
% respuesta en frecuencia: 0.0e+00

%% Verificación 2: los polos, por la fórmula cerrada
p  = pole(F);
pm = ((2+Kp) + [1; -1]*1j*sqrt(4*Ki - Kp^2)) / (2*(1+Kp+Ki));
fprintf('polos: %.6f +- %.6fj, |p| = %.6f, formula - pole = %.1e\n', ...
        real(pm(1)), imag(pm(1)), 1/sqrt(1+Kp+Ki), max(abs(sort(pm) - sort(p))))
% polos: 0.997463 +- 0.006581j, |p| = 0.997485, formula - pole = 1.1e-14

%% Verificación 3: el equivalente continuo
zeta = Kp/(2*sqrt(Ki));  wn = sqrt(Ki);
fprintf('zeta = %.4f, wn = %.5f rad/muestra, fn = %.2f MHz\n', zeta, wn, wn*fs/(2*pi)/1e6)
% zeta = 0.3536, wn = 0.00707 rad/muestra, fn = 18.01 MHz
sc = roots([1 Kp Ki]);                     % polos continuos, en rad/muestra
fprintf('1/(1-s) - pole = %.1e;  exp(s) - pole = %.1e\n', ...
        max(abs(sort(1./(1-sc)) - sort(p))), max(abs(sort(exp(sc)) - sort(p))))
% 1/(1-s) - pole = 1.1e-14;  exp(s) - pole = 2.5e-05

%% Lo que dicen las figuras
S = stepinfo(F);
fprintf('sobrepico %.2f %%, asentamiento %.0f muestras, -3 dB en %.1f MHz\n', ...
        S.Overshoot, S.SettlingTime/Ts, bandwidth(F)/(2*pi)/1e6)
% sobrepico 39.72 %, asentamiento 1477 muestras, -3 dB en 30.3 MHz

Tres detalles que muerden:

• El eje de tiempo está en segundos. Como z se creó con Ts, step y stepinfo devuelven tiempos en segundos: el asentamiento sale 9,23·10⁻⁸ y hay que dividir por Ts para leer 1477 muestras. Lo mismo bandwidth, que da rad/s.
• Restar dos tf no es comparar. F - Fm arma una transferencia nueva con denominador producto, y sus coeficientes no dicen si las dos eran iguales. Lo que compara es evaluar las dos por separado, con freqresp, o los coeficientes normalizados por el principal: feedback no garantiza que el denominador empiece en 1.
• Con latencia, pzmap muestra un polo y un cero de más en el origen. feedback no cancela factores comunes: con Lat = 5, pole(F) devuelve siete raíces y una es $z=0$, tapada por un cero en el mismo lugar: en el gráfico no se distinguen.

Para tener a mano

1. El lazo de FCR corre una muestra por símbolo, $f_s=\mathrm{BR}$, con $K_{\mathrm{PD}}{K}_V=1$ por muestra adentro de $K_p$ y $K_i$.
2. El PI es $K_p+K_{i}z/(z-1)$ y el NCO es $z/(z-1)$: dos acumuladores sin retardo, ${(z-1)}^2$ en el lazo abierto, tipo 2.
3. Lazo cerrado:$$F(z)=\frac{z[(K_p+K_i)z-K_p]}{(1+K_p+K_i)z^2-(2+K_p)z+1}$$y el error, $E/\Theta$, tiene un cero doble en $z=1$: ni escalón ni rampa dejan error.
4. Polos complejos si y sólo si $\zeta <1$; módulo $1/\sqrt{1+K_p+K_i}$: sin latencia, estable con cualquier ganancia.
5. Es el lazo continuo por Euler hacia atrás, $s\to 1-z^{-1}$: ${\omega }_n=\sqrt{K_i}$ y $\zeta =K_p/(2\sqrt{K_i})$ por muestra. Vale mientras $B_n{T}_s$ sea chico: con $\zeta =\mathrm{0,3536}$, el sobrepico discreto se aparta 0,36 puntos con 3,75·10⁻³ y 6,5 con 0,075, siempre hacia menos.
6. La clase: $\zeta =\mathrm{0,3536}$, $f_n=\mathrm{18,01}$ MHz, $B_n{T}_s=\mathrm{3,75}·{10}^{-3}$; sobrepico 39,72 %, asentamiento 1477 muestras.
7. La latencia es $z^{-\mathrm{Lat}}$ en la realimentación: $\mathrm{Lat}+1$ polos. Lat = 0 es un lazo algebraico que el hardware no puede hacer. Con el lazo de la clase, 100 símbolos lo vuelven inestable.