Ancho de banda, adquisición y tracking

Una vez elegido el filtro PI, al lazo le queda una sola decisión de fondo: cuánto ancho de banda darle. Todo lo que hace el PLL depende de ese número. Con más ancho de banda engancha antes y desde más lejos, porque la filmina da un tiempo de enganche en frecuencia que cae como $1/B_n^3$. Pero deja pasar más ruido a la fase estimada, porque la varianza del error crece como $B_n$. El ancho que importa para el ruido no es el de −3 dB sino el equivalente de ruido $B_n$, y con él se arman también $K_p$ y $K_i$.

01 · El ancho de −3 dB

Con el filtro PI, $L(s)=K_p+K_i/s$, el PLL linealizado es de segundo orden (el filtro de lazo):

$$F(s)=\frac{2\zeta {\omega }_{n}s+{\omega }_n^2}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2}$$ $$\zeta =\frac{K_p}{2}\sqrt{\frac{K_V{K}_{\mathrm{PD}}}{K_i}},{\omega }_n=\sqrt{K_V{K}_{\mathrm{PD}}{K}_i}$$

Es un pasabajos de fase: copia lo que $\theta$ hace despacio y deja pasar como error lo que hace rápido. Su primer ancho de banda es el de siempre, donde la ganancia de potencia cae a la mitad.

Paso 1: el módulo al cuadrado. Con $s=j\omega$, el numerador es ${\omega }_n^2+j2\zeta {\omega }_n\omega$ y el denominador $({\omega }_n^2-{\omega }^2)+j2\zeta {\omega }_n\omega$:

$${|F(j\omega )|}^2=\frac{{\omega }_n^4+4{\zeta }^2{\omega }_n^2{\omega }^2}{{({\omega }_n^2-{\omega }^2)}^2+4{\zeta }^2{\omega }_n^2{\omega }^2}$$

Paso 2: igualarlo a 1/2. Pasando el denominador multiplicando, con $x={\omega }^2/{\omega }_n^2$ y dividiendo todo por ${\omega }_n^4$, queda $2+8{\zeta }^{2}x={(1-x)}^2+4{\zeta }^{2}x$, o sea una cuadrática en $x$:

$$x^2-2(1+2{\zeta }^2)x-1=0$$

Paso 3: la raíz positiva. El producto de las raíces es $-1$, así que una sola es positiva, y es la que sirve porque $x$ es un cuadrado:

$${\omega }_{3\mathrm{dB}}={\omega }_n\sqrt{1+2{\zeta }^2+\sqrt{{(1+2{\zeta }^2)}^2+1}}$$

Es la fórmula de la filmina. Con las constantes de phase_system.m, $K_p=\mathrm{0,005}$, $K_i=K_p/100$ y $K_{\mathrm{PD}}{K}_V=1$ por muestra, sale $\zeta =\mathrm{0,3536}$ y ${\omega }_n=\mathrm{0,00707}$ rad por muestra, o sea $f_n=\mathrm{18,01}$ MHz a $f_s=16$ GHz. La fórmula da $f_{3\mathrm{dB}}=\mathrm{30,40}$ MHz y bandwidth(F), sobre el lazo discreto del script, 30,34 MHz.

02 · El ancho de banda equivalente de ruido, paso a paso

Para el ruido, el ancho de −3 dB no dice cuánto pasa: el ruido blanco pasa por todas las frecuencias, y lo que llega a la salida depende del área entera bajo ${\left|F\right|}^2$, pico y faldas incluidos. El ancho de banda equivalente de ruido resume esa área en un número: es el ancho del pasabajos rectangular ideal que deja pasar la misma potencia de ruido blanco que el PLL. Es el mismo concepto de la Unidad 3, con una diferencia de normalización: el rectángulo de la filmina tiene la altura de la ganancia en continua, ${|F(0)|}^2=1$, no la del pico, y va de $-B_n$ a $B_n$. Como ${\left|F\right|}^2$ es par, igualar las áreas es

$$B_n={\int }_0^{\infty }{|F(j2\pi f)|}^{2}df$$

Con un pico, la normalización importa. El lazo de la clase tiene ${\left|F\right|}^2$ de hasta 3,14 en su pico; normalizando por el pico, el mismo lazo tendría un ancho equivalente 3,14 veces más chico y la misma potencia de ruido.

Paso 1: separar el numerador. El módulo al cuadrado de la sección anterior tiene un denominador $D(\omega )={({\omega }_n^2-{\omega }^2)}^2+4{\zeta }^2{\omega }_n^2{\omega }^2$ y un numerador con dos términos, así que la integral en $\omega$ se parte en dos:

$${\int }_0^{\infty }{\left|F\right|}^{2}d\omega ={\omega }_n^4\underset{I_0}{\underbrace{{\int }_0^{\infty }\frac{d\omega }{D}}}+4{\zeta }^2{\omega }_n^2\underset{I_2}{\underbrace{{\int }_0^{\infty }\frac{{\omega }^{2}d\omega }{D}}}$$

Paso 2: las dos integrales de un segundo orden. $I_0$ es el área bajo el módulo al cuadrado de $1/(s^2+2\zeta {\omega }_{n}s+{\omega }_n^2)$ e $I_2$ la de $s/(s^2+2\zeta {\omega }_{n}s+{\omega }_n^2)$. Por residuos, cerrando en el semiplano superior, dan

$$I_0=\frac{\pi }{4\zeta {\omega }_n^3},I_2=\frac{\pi }{4\zeta {\omega }_n}$$

Paso 3: sumar. ${\omega }_n^4{I}_0=\pi {\omega }_n/(4\zeta )$ y $4{\zeta }^2{\omega }_n^2{I}_2=\pi \zeta {\omega }_n$, así que ${\int }_0^{\infty }{\left|F\right|}^{2}d\omega =\pi {\omega }_n(\zeta +1/4\zeta )$.

Paso 4: pasar a hertz. La definición integra en $f$, y $d\omega =2\pi df$. Dividir por $2\pi$ da la fórmula de la filmina:

$$B_n=\frac{{\omega }_n}{2}(\zeta +\frac{1}{4\zeta })$$

El paso 4 es el que se olvida. ${\omega }_n$ entra en radianes por segundo y $B_n$ sale en hertz: el $2\pi$ ya está adentro. En el lazo discreto, con ${\omega }_n$ en radianes por muestra, sale en ciclos por muestra, que es $B_n{T}_s$. Para la clase, $B_n{T}_s=\mathrm{3,75}·{10}^{-3}$, o sea $B_n=\mathrm{60,00}$ MHz. La integral numérica con integral da el mismo $\mathrm{3,750000}·{10}^{-3}$ que la cerrada. Y en el lazo discreto del script, donde el ancho se mide por Parseval como la mitad de la energía de la respuesta al impulso, $B_n{T}_s=\frac{1}{2}{\sum }_{n}f{[n]}^2$, da 59,55 MHz: un 0,74 % menos, por ser discreto.

Contra el de −3 dB, en las mismas unidades. Los dos son proporcionales a ${\omega }_n$, así que su cociente depende sólo de $\zeta$, y sale de las dos formas cerradas:

$B_n$ es siempre más ancho que el de −3 dB, y la diferencia crece cuando baja $\zeta$, porque el pico suma área. Es el gráfico de la filmina con $f_n=100$ kHz, normalizado. Para un ${\omega }_n$ dado, $B_n$ tiene su mínimo donde se anula la derivada de $\zeta +1/4\zeta$, que es $1-1/4{\zeta }^2$: en $\zeta =\mathrm{0,5}$, con $B_n={\omega }_n/2$. Y la expresión no cambia si se reemplaza $\zeta$ por $1/4\zeta$: el 0,354 de la clase y el 0,707 de los libros son ese par, y con el mismo ${\omega }_n$ dejan pasar exactamente el mismo ruido.

03 · La receta de las constantes

El diseño arranca al revés: se eligen $B_n$ y $\zeta$, y de ahí salen las constantes del filtro.

Paso 1: igualar denominadores. El PI cerrado da $s^2+K_V{K}_{\mathrm{PD}}{K}_{p}s+K_V{K}_{\mathrm{PD}}{K}_i$ abajo; comparado con la forma de segundo orden, $K_p{K}_V{K}_{\mathrm{PD}}=2\zeta {\omega }_n$ y $K_i{K}_V{K}_{\mathrm{PD}}={\omega }_n^2$.

Paso 2: despejar ${\omega }_n$ de $B_n$. ${\omega }_n=2B_n/(\zeta +1/4\zeta )$.

Paso 3: reemplazar.

$$K_p{K}_V{K}_{\mathrm{PD}}=\frac{4\zeta B_n}{\zeta +\frac{1}{4\zeta }},K_i{K}_V{K}_{\mathrm{PD}}=\frac{4B_n^2}{{(\zeta +\frac{1}{4\zeta })}^2}$$

La filmina escribe $2\zeta B_n$ en el numerador de la primera, pero la cuenta da 4: $2\zeta ·2B_n/(\dots )$. La segunda está bien. Con el 2, $K_p$ sale a la mitad, y como $K_i$ no cambia, el lazo queda con la mitad del $\zeta$ pedido.

Aplicada a la clase, con $B_n{T}_s=\mathrm{3,75}·{10}^{-3}$ y $\zeta =\mathrm{0,3536}$, la receta devuelve las constantes del script: $K_p=\mathrm{0,0050}$ y $K_i=\mathrm{5,00}·{10}^{-5}$. Con el 2 de la filmina, $K_p=\mathrm{0,0025}$: un lazo con $\zeta =\mathrm{0,177}$ y 90,0 MHz de ancho de ruido en lugar de 60. Hay un solo $\zeta$ en el que el error no cambia $B_n$: el 0,707, porque su mitad es su espejo $1/4\zeta$. Y la receta de la filmina con $\zeta =\mathrm{0,707}$ da justo $K_p=\mathrm{0,005}$ y $K_i=5·{10}^{-5}$.

La receta es continua y el lazo de la clase es discreto. Mientras $B_n{T}_s$ sea chico, el ancho que se pide es el que se obtiene. Medido como $\frac{1}{2}{\sum }_{n}f{[n]}^2$ sobre el lazo de phase_system.m armado con la receta y $\zeta =\mathrm{0,707}$:

El discreto siempre queda un poco por debajo, y la diferencia crece con el ancho: por debajo del 1 % hasta $B_n{T}_s\approx \mathrm{0,01}$, y del 12 % con 0,1. Para anchos así hace falta la receta discreta, que es otro tema (el PLL en tiempo discreto).

04 · Tracking: la varianza de la fase

Enganchado, el lazo sigue a la entrada, pero la entrada trae ruido y una parte de ese ruido llega a la fase estimada. La filmina mide la calidad del seguimiento con la varianza del error, ${\sigma }_{{\theta }_e}^2=\mathrm{E}\{{|\theta -\widehat{\theta }|}^2\}$, y la da para un tono de potencia $P_{\mathrm{in}}$ en ruido blanco de densidad $N_0/2$.

Paso 1: el ruido como fase. El tono es $A\cos ({\omega }_{v}t+\theta )$, con $P_{\mathrm{in}}=A^2/2$. El ruido alrededor de la portadora se separa en una componente en fase con el tono y otra en cuadratura, $n_s$, cada una con densidad $N_0$. La que está en fase sólo le cambia el largo; la que está en cuadratura lo hace girar, y con ruido chico lo gira $n_s/A$ radianes. Para el lazo, entonces, el ruido es una fase extra de densidad $N_0/A^2=N_0/(2P_{\mathrm{in}})$, sumada a $\theta$.

Paso 2: lo filtra $F$. Esa fase extra entra al lazo por el mismo lugar que $\theta$, así que llega a $\widehat{\theta }$ multiplicada por ${\left|F\right|}^2$, y su potencia es la densidad por el área, que por la sección 02 vale $2B_n$:

$${\sigma }_{{\theta }_e}^2=\frac{N_0}{2P_{\mathrm{in}}}{\int }_{-\infty }^{\infty }{|F(j2\pi f)|}^{2}df=\frac{N_0{B}_n}{P_{\mathrm{in}}}$$

Con $\theta$ constante el error es sólo ese ruido filtrado, cambiado de signo, y la fórmula vale para el error. Es el motivo para preferir $B_n$ a $f_{3\mathrm{dB}}$: la varianza de la fase sale en una línea, y depende de la forma del lazo sólo a través de $B_n$. Más potencia de señal la baja; más ruido o más ancho de banda la suben, y la suben en proporción.

En banda base y en muestras, con un tono de potencia 1 y ruido complejo de varianza ${\sigma }_w^2$ por muestra, $N_0={\sigma }_w^2{T}_s$ y la fórmula queda ${\sigma }_{{\theta }_e}^2=B_n{T}_s·{\sigma }_w^2$. Simulado con el lazo de la clase, el detector de producto y dos millones de muestras:

La varianza medida es la de la fórmula dentro de un 3,2 %, en tres décadas de ruido. Dos cosas la separan. El lazo discreto deja pasar un 0,74 % menos que el $B_n$ continuo. Y la propia medición se mueve, porque el error está correlacionado a lo largo de cientos de muestras: el mismo lazo con 10 dB, simulado con ocho semillas de 800 000 muestras, da entre 0,943 y 1,032 veces la fórmula. Con 10 dB por muestra el lazo de la clase deja la fase con 1,11° de desvío. Una muestra sola tiene una varianza de fase de ${\sigma }_w^2/2$, y el lazo la divide por $1/(2B_n{T}_s)=133$: promedia el ruido de más de cien muestras.

05 · Adquisición: primero la frecuencia, después la fase

Todo lo anterior supone el lazo enganchado. Cuando se prende, el oscilador está corrido $\mathrm{\Delta }f$ de la entrada, el error de fase crece sin parar y el detector recorre su curva S entera, una vuelta tras otra. Ahí el modelo lineal no vale. La filmina parte la adquisición en dos etapas.

Enganche en frecuencia. El VCO se ajusta hasta igualar la frecuencia de la entrada. Mientras el error da vueltas, el detector promedia casi cero, pero no cero: el error gira más despacio en la mitad de la curva que lo acerca, y el integrador del PI acumula ese desbalance vuelta tras vuelta. Por eso es lento, y por eso necesita el integrador. La filmina lo estima con

$$T_{\mathrm{FL}}\approx 4\frac{{(\mathrm{\Delta }f)}^2}{B_n^3}$$

Enganche en fase. Cuando la frecuencia ya coincide queda un error de fase fijo, que para el lazo es un escalón: lo corrige en la escala de tiempo del lazo lineal, $T_{\mathrm{PL}}\approx \mathrm{1,3}/B_n$, y el total es $T_{\mathrm{LOCK}}\approx T_{\mathrm{FL}}+T_{\mathrm{PL}}$. Las unidades tienen que ser coherentes: con $\mathrm{\Delta }f$ y $B_n$ en hertz los tiempos salen en segundos; con los dos en ciclos por muestra, en muestras.

Las dos fórmulas son aproximaciones sin derivación en la filmina, así que se cruzan con simulación. El lazo es el de phase_system.m en el tiempo, sin retardo, con un detector no lineal, y la entrada salta $\mathrm{\Delta }f=200$ MHz, que con $B_n=60$ MHz pide $T_{\mathrm{FL}}\approx$ 11 852 muestras, 741 ns. $T_{\mathrm{FL}}$ medido es la última vuelta del error; el enganche, la última vez que el error pasa de 0,1 rad. Con el mismo $B_n$ y dos $\zeta$:

Con el detector de producto, que tiene una curva S senoidal, y $\zeta =\mathrm{0,707}$, la fórmula acierta: 11 268 muestras contra 11 852, un 5 % menos, y el enganche total, 12 287 contra 12 199. Pero el 4 no es una constante del PLL. La forma clásica de la que sale, que la filmina no muestra, es $T_{\mathrm{FL}}\approx {(2\pi \mathrm{\Delta }f)}^2/(2\zeta {\omega }_n^3)$. Con ${\omega }_n=2B_n/(\zeta +1/4\zeta )$ el coeficiente da 4,16 para $\zeta =\mathrm{0,707}$ y 8,33 para el 0,354 de la clase: con el mismo $B_n$ los dos tienen el mismo ${\omega }_n$, y el tiempo va como $1/\zeta$. La simulación da 1,94 veces la fórmula. Con el mismo ancho de ruido, el lazo de la clase tarda el doble en engancharse en frecuencia.

El detector también cuenta. Con el diente de sierra, el detector ideal que dibuja la filmina, el mismo salto engancha 4,1 veces antes con 0,354 y 11,5 veces antes con 0,707: el 4 es de un detector senoidal. Y por debajo de cierto salto no hay enganche en frecuencia que medir. Buscado por bisección, el mayor salto que el lazo de la clase alcanza sin dar ninguna vuelta es de 45,4 MHz con el producto y de 89,3 MHz con el diente de sierra. Ahí $T_{\mathrm{FL}}$ es cero, aunque la fórmula dé un número.

El 1,3 de $T_{\mathrm{PL}}$ también es de un $\zeta$. En el escalón de fase lineal, $\mathrm{1,3}/B_n$ cae casi justo en el mínimo del error de un lazo con $\zeta =1$, −13,5 % del escalón, y de ahí en adelante ese error sólo decae. Con otro $\zeta$ el mismo instante deja otro resto: −4,5 % con $\zeta =2$, −20,3 % con 0,707, −29,4 % con 0,5 y −39,7 % con el 0,354 de la clase, que recién queda por debajo del 13,5 % en $\mathrm{2,15}/B_n$.

El rango de captura. La filmina da hasta dónde el lazo puede engancharse solo:

$${(\mathrm{\Delta }f)}_{\text{pull-in}}\approx \left(2\pi \sqrt{2}\zeta \right)B_n$$

Para la clase da $2\pi \sqrt{2}·\mathrm{0,354}·60$ MHz = 188 MHz. Con el detector de producto la simulación no encuentra ese borde: engancha un salto de 1,6 GHz después de 106 267 vueltas y 99 µs, y uno de 6,4 GHz después de 10,75 millones de vueltas y 2,3 ms. Con un integrador que no pierde, el empuje de cada vuelta es chico pero no cero, y se acumula sin límite: lo que crece con el salto es el precio, como ${(\mathrm{\Delta }f)}^2$. El único techo que aparece es el del muestreo: a 7,84 GHz, cerca de $f_s/2$, no enganchó en 109 millones de muestras. Con el diente de sierra sí hay un borde, y es abrupto: con 0,354 engancha 592 MHz en 33 446 muestras, y a menos de 0,1 MHz de ahí ya no engancha en dos millones. Con 0,707 el borde está en 699 MHz. Son 3,1 y 1,9 veces lo que da la fórmula, 188 y 377 MHz. El rango de captura de la filmina es un orden de magnitud, no una frontera.

El compromiso. Juntas, las dos mitades de la filmina dicen lo mismo desde lados opuestos. Más $B_n$ acorta $T_{\mathrm{FL}}$ como $1/B_n^3$, $T_{\mathrm{PL}}$ como $1/B_n$ y agranda el rango de captura en proporción; y agranda ${\sigma }_{{\theta }_e}^2$ en proporción. Se elige $B_n$ con dos datos: cuánto puede estar corrido el oscilador y cuánta varianza de fase tolera la constelación.

Arriba, mover $\zeta$ cambia la forma de ${\left|F\right|}^2$ y no su área. Con 0,354 el pico llega a 3,14 y el −3 dB cae a la mitad de $B_n$; con 2 casi no hay pico, pero la falda se estira. Mover $B_n$ no cambia nada arriba, porque el eje está en unidades de $B_n$, y cambia todo abajo. De 60 a 30 MHz, el salto de 200 MHz pasa de 196 vueltas a 1 623 y $T_{\mathrm{FL}}$ de 1,44 a 12,1 µs, 8,4 veces más; el desvío de la fórmula baja de 1,11° a 0,79°, un factor $\sqrt{2}$. De 60 a 120 MHz quedan 21 vueltas y 146 ns. Desde 320 MHz el error ya no da ninguna vuelta: $T_{\mathrm{FL}}$ medido es cero y la fórmula sigue dando 4,9 ns. Con el mismo $B_n$, pasar de 0,354 a 0,707 parte $T_{\mathrm{FL}}$ a la mitad sin tocar la varianza. Y con el diente de sierra el lazo de la clase engancha el mismo salto en 27 vueltas y 182 ns.

06 · El análisis en MATLAB

El lazo de la clase de punta a punta: sus dos anchos de banda, $B_n$ por integración y por Parseval, la receta, una simulación de tracking con ruido y una de adquisición. Lo que figura como comentario es la salida real de correrlo en R2026a. La función del final es el lazo de phase_system.m en el tiempo.

%% 1. El lazo de la clase, como lo arma phase_system.m
fs = 16e9; Ts = 1/fs;
z = tf('z', Ts);
Kp = 0.005; Ki = Kp/100;
L = Kp + Ki*z/(z-1);
NCO = z/(z-1);
F = feedback(L*NCO, 1);

%% 2. El equivalente continuo: zeta, wn, Bn y el ancho de -3 dB
zeta = Kp/(2*sqrt(Ki));               % con K_PD*K_V = 1 por muestra
wn = sqrt(Ki);                        % rad/muestra
Bn = wn/2*(zeta + 1/(4*zeta));        % ciclos/muestra: es Bn*Ts
w3 = wn*sqrt(1 + 2*zeta^2 + sqrt((1 + 2*zeta^2)^2 + 1));
fprintf('zeta = %.4f   fn = %.2f MHz   Bn*Ts = %.3e   Bn = %.2f MHz\n', ...
        zeta, wn/(2*pi)*fs/1e6, Bn, Bn*fs/1e6)
fprintf('f3dB cerrada = %.2f MHz   bandwidth(F) = %.2f MHz\n', ...
        w3/(2*pi)*fs/1e6, bandwidth(F)/(2*pi)/1e6)
% zeta = 0.3536   fn = 18.01 MHz   Bn*Ts = 3.750e-03   Bn = 60.00 MHz
% f3dB cerrada = 30.40 MHz   bandwidth(F) = 30.34 MHz

%% 3. Bn por integración numérica de |F|^2, en ciclos por muestra
s = @(f) 2j*pi*f;
F2 = @(f) abs((2*zeta*wn*s(f) + wn^2)./(s(f).^2 + 2*zeta*wn*s(f) + wn^2)).^2;
Bn_int = integral(F2, 0, Inf, 'RelTol', 1e-12, 'AbsTol', 1e-15);
fprintf('integral = %.6e   cerrada = %.6e   diferencia = %.1e\n', Bn_int, Bn, abs(Bn_int - Bn))
% integral = 3.750000e-03   cerrada = 3.750000e-03   diferencia = 0.0e+00

%% 4. Bn del lazo discreto: Parseval, la mitad de la energía de f[n]
[b, a] = tfdata(F, 'v');
f = filter(b, a, [1 zeros(1, 2e5)]);
Bn_d = sum(f.^2)/(2*sum(f)^2);        % sum(f) = F(1) = 1
fprintf('discreto: Bn*Ts = %.4e   Bn = %.2f MHz\n', Bn_d, Bn_d*fs/1e6)
% discreto: Bn*Ts = 3.7221e-03   Bn = 59.55 MHz

%% 5. La receta: de (Bn, zeta) a (Kp, Ki)
kp = 4*zeta*Bn/(zeta + 1/(4*zeta));
ki = 4*Bn^2/(zeta + 1/(4*zeta))^2;
fprintf('receta(Bn, %.4f): Kp = %.4f   Ki = %.2e\n', zeta, kp, ki)
% receta(Bn, 0.3536): Kp = 0.0050   Ki = 5.00e-05
kp2 = 2*zeta*Bn/(zeta + 1/(4*zeta));  % con el 2 de la filmina
z2 = kp2/(2*sqrt(ki));
fprintf('con 2*zeta*Bn: Kp = %.4f -> zeta = %.4f, Bn = %.1f MHz\n', ...
        kp2, z2, sqrt(ki)/2*(z2 + 1/(4*z2))*fs/1e6)
% con 2*zeta*Bn: Kp = 0.0025 -> zeta = 0.1768, Bn = 90.0 MHz

%% 6. Cuánto se separa el discreto de la receta continua (zeta = 0.7071)
for BT = [0.001 0.00375 0.01 0.03 0.1]
    ze = 1/sqrt(2);
    k1 = 4*ze*BT/(ze + 1/(4*ze)); k2 = 4*BT^2/(ze + 1/(4*ze))^2;
    Fd = feedback((k1 + k2*z/(z-1))*z/(z-1), 1);
    [bd, ad] = tfdata(Fd, 'v');
    fd = filter(bd, ad, [1 zeros(1, 4e5)]);
    fprintf('Bn*Ts pedido %.5f   medido %.5f   cociente %.4f\n', BT, sum(fd.^2)/2, sum(fd.^2)/2/BT)
end
% Bn*Ts pedido 0.00100   medido 0.00100   cociente 0.9987
% Bn*Ts pedido 0.00375   medido 0.00373   cociente 0.9950
% Bn*Ts pedido 0.01000   medido 0.00987   cociente 0.9868
% Bn*Ts pedido 0.03000   medido 0.02884   cociente 0.9613
% Bn*Ts pedido 0.10000   medido 0.08801   cociente 0.8801

%% 7. Tracking: un tono de potencia 1 en ruido blanco, detector de producto
rng(1)
for snr = [10 20 30]                  % dB por muestra
    sw2 = 10^(-snr/10);               % varianza del ruido complejo: N0*fs
    e = pll(Kp, Ki, 0, sw2, 2e6, 'producto');
    v = mean(e(5e4:end).^2);          % sin el transitorio
    fprintf('SNR %2d dB: var = %.3e   N0*Bn/Pin = %.3e   cociente %.3f\n', snr, v, Bn*sw2, v/(Bn*sw2))
end
% SNR 10 dB: var = 3.727e-04   N0*Bn/Pin = 3.750e-04   cociente 0.994
% SNR 20 dB: var = 3.630e-05   N0*Bn/Pin = 3.750e-05   cociente 0.968
% SNR 30 dB: var = 3.670e-06   N0*Bn/Pin = 3.750e-06   cociente 0.979

%% 8. Adquisición: un salto de 200 MHz, sin ruido
df = 200e6/fs;
for ze = [zeta 1/sqrt(2)]
    for det = ["producto" "argumento"]
        k1 = 4*ze*Bn/(ze + 1/(4*ze)); k2 = 4*Bn^2/(ze + 1/(4*ze))^2;
        e = pll(k1, k2, df, 0, 3e5, det);
        k = round(e(end)/(2*pi)); e = e - 2*pi*k;          % la vuelta en la que quedó
        TFL = max([0 find(abs(e) >= pi, 1, 'last')]);       % la última vuelta
        TL = max([0 find(abs(e) > 0.1, 1, 'last')]);        % enganchado a 0,1 rad
        fprintf('zeta %.3f %-9s: %3d vueltas  TFL = %5d (4*df^2/Bn^3 = %5.0f)  TLOCK = %5d (+ 1.3/Bn: %5.0f)\n', ...
                ze, det, abs(k), TFL, 4*df^2/Bn^3, TL, 4*df^2/Bn^3 + 1.3/Bn)
    end
end
% zeta 0.354 producto : 196 vueltas  TFL = 23028 (4*df^2/Bn^3 = 11852)  TLOCK = 24738 (+ 1.3/Bn: 12199)
% zeta 0.354 argumento:  27 vueltas  TFL =  2907 (4*df^2/Bn^3 = 11852)  TLOCK =  4324 (+ 1.3/Bn: 12199)
% zeta 0.707 producto :  96 vueltas  TFL = 11268 (4*df^2/Bn^3 = 11852)  TLOCK = 12287 (+ 1.3/Bn: 12199)
% zeta 0.707 argumento:  10 vueltas  TFL =  1033 (4*df^2/Bn^3 = 11852)  TLOCK =  1995 (+ 1.3/Bn: 12199)

function eu = pll(kp, ki, df, sw2, N, det)
    % El lazo de phase_system.m en el tiempo: L = Kp + Ki*z/(z-1) y
    % NCO = z/(z-1), sin retardo. Como el NCO ya suma la corrección de esta
    % muestra, el error sale de resolver e = W(d - K*e), con K = Kp + Ki.
    K = kp + ki; ph = 0; ac = 0; eu = zeros(1, N);
    w = sqrt(sw2/2)*(randn(1, N) + 1j*randn(1, N));
    for n = 1:N
        th = 2*pi*df*(n-1);
        r = exp(1j*th) + w(n);
        d = angle(r) - ph - ac;
        if det == "argumento"         % W = diente de sierra
            e = (mod(d + pi, 2*pi) - pi)/(1 + K);
        else                          % W = |r| sen(.), parte imaginaria
            A = abs(r); e = 0;
            for it = 1:2
                g = e - A*sin(d - K*e);
                e = e - g/(1 + A*K*cos(d - K*e));
            end
        end
        ac = ac + ki*e;
        ph = ph + kp*e + ac;
        eu(n) = th - ph;              % theta - phi, sin envolver
    end
end

Tres detalles que muerden:

• $B_n$ sale en ciclos, no en radianes. Con ${\omega }_n$ en rad/muestra, wn/2*(zeta + 1/(4*zeta)) es $B_n{T}_s$: se multiplica por fs, no por fs/(2*pi). Dividir por $2\pi$ da 9,55 MHz, un lazo $2\pi$ veces más angosto que el verdadero. El bloque 3 lo comprueba integrando en $f$.
• El lazo del script no tiene retardo. L y NCO suman la muestra actual, así que la fase de salida depende del error de esa misma muestra. Para feedback no es problema; en el tiempo hay que despejar el error, y por eso la función resuelve $e=W(d-Ke)$: en una línea con el diente de sierra, con dos pasos de Newton con el seno.
• La varianza se mide enganchado. Las primeras muestras son el transitorio del filtro; el bloque 7 las tira antes de promediar. Y se mide $\theta -\widehat{\theta }$, no la salida del detector, que trae el ruido entero sin filtrar.

Para tener a mano

1. ${\omega }_{3\mathrm{dB}}={\omega }_n\sqrt{1+2{\zeta }^2+\sqrt{{(1+2{\zeta }^2)}^2+1}}$, de igualar ${\left|F\right|}^2$ a 1/2. Para la clase, 30,40 MHz.
2. $B_n={\int }_0^{\infty }{|F(j2\pi f)|}^{2}df=\frac{{\omega }_n}{2}(\zeta +\frac{1}{4\zeta })$ en hertz con ${\omega }_n$ en rad/s. Para la clase, $B_n{T}_s=\mathrm{3,75}·{10}^{-3}$ ciclos por muestra: 60 MHz, y 59,55 en el lazo discreto.
3. $B_n$ siempre es más ancho que el de −3 dB: 1,97 veces en la clase. Mínimo en $\zeta =\mathrm{0,5}$ para un ${\omega }_n$ dado, e igual para $\zeta$ y $1/4\zeta$.
4. Receta: $K_p{K}_V{K}_{\mathrm{PD}}=4\zeta B_n/(\zeta +1/4\zeta )$ —con 4, no con el 2 de la filmina— y $K_i{K}_V{K}_{\mathrm{PD}}=4B_n^2/{(\zeta +1/4\zeta )}^2$. Con 60 MHz y 0,354 devuelve 0,005 y $5·{10}^{-5}$.
5. ${\sigma }_{{\theta }_e}^2=N_0{B}_n/P_{\mathrm{in}}$: la densidad de la fase de ruido por el área de ${\left|F\right|}^2$. En muestras, $B_n{T}_s·{\sigma }_w^2$. Con 10 dB por muestra, el lazo de la clase deja 1,11° de desvío.
6. $T_{\mathrm{FL}}\approx 4{(\mathrm{\Delta }f)}^2/B_n^3$ acierta con detector senoidal y $\zeta =\mathrm{0,707}$, dentro del 5 %. El 4 es de ese caso: con el 0,354 de la clase tarda 1,94 veces más, con diente de sierra entre 4 y 11 veces menos, y por debajo de 45 MHz el lazo de la clase engancha sin dar vueltas.
7. $T_{\mathrm{PL}}\approx \mathrm{1,3}/B_n$ es el de $\zeta =1$; con el 0,354 de la clase el error todavía es el 40 % del escalón en ese instante.
8. El rango de captura, $(2\pi \sqrt{2}\zeta )B_n$, da 188 MHz para la clase y no es un borde: con detector senoidal el lazo engancha 6,4 GHz, pagando 2,3 ms; con diente de sierra el borde existe, pero en 592 MHz.