Más muestras no arreglan el periodograma

El apunte 5 llega a la densidad espectral de potencia por el camino corto: si el proceso es estacionario en sentido amplio, su autocorrelación depende de un solo retardo, y la transformada de esa autocorrelación es la PSD. Estimarla parece igual de corto —tomar $N$ muestras, transformarlas, quedarse con el módulo al cuadrado— y ese estimador tiene nombre: el periodograma. Lo que sigue es el reflejo que todo el mundo tiene delante de una medición ruidosa: juntá más datos. Con una media muestral funciona, y el error baja como $1/\sqrt{N}$. Con el periodograma no baja nada. Cien veces más muestras dejan el error relativo exactamente donde estaba, y eso no es un defecto de la simulación: es una propiedad del estimador.

Medido sobre ruido blanco, 400 repeticiones, mirando siempre el mismo bin: con $N=1000$ el desvío del valor estimado vale 1,050 veces su media; con $N=10000$, 0,985; con $N=100000$, 0,919. El estimador es tan grande como su propia incertidumbre, y lo sigue siendo con cien veces más datos. Lo que sí la baja es partir la serie en $K$ tramos y promediar sus periodogramas, y la baja exactamente como $1/\sqrt{K}$: 0,496 contra 0,500 con cuatro tramos, 0,256 contra 0,250 con dieciséis, 0,128 contra 0,125 con sesenta y cuatro. Lo que sigue es por qué una de las dos cosas no hace nada y la otra sí.

01 · De la autocorrelación a la PSD

Un proceso estacionario en sentido amplio —WSS— pide dos cosas: que su media no dependa del instante, y que su autocorrelación dependa sólo de la distancia entre los dos instantes, no de dónde caen. La segunda es la que compra todo lo demás: $R_{XX}(\tau )=E[X[n]X[n+\tau ]]$ es una función de una sola variable, y no de dos. Con eso el proceso deja de ser una familia infinita de variables aleatorias y pasa a ser algo que se puede dibujar.

La PSD es la transformada de esa autocorrelación —el teorema de Wiener-Khinchin—, y es la forma en que la unidad mide «cuánta potencia hay en cada frecuencia»:

$$S_{XX}(f)=\sum _{\tau =-\infty }^{\infty }{R}_{XX}(\tau )e^{-j2\pi f\tau }$$

El ruido blanco es el caso más simple y el que usa todo lo que viene: no se acuerda de nada, así que su autocorrelación es un solo impulso en el origen, $R_{XX}(\tau )={\sigma }^2\delta [\tau ]$. La transformada de un impulso es una constante, así que su PSD es plana y vale ${\sigma }^2$ en toda frecuencia. Con varianza unitaria ese nivel es $1$, y es el único número que hay que tener en la cabeza para leer todo lo que sigue: la respuesta correcta es una recta horizontal en uno. Cada vez que una estimación se aparta de esa recta, se está mirando error puro.

De la definición no se puede calcular nada: pide $R_{XX}$ exacta, que es una esperanza sobre el ensamble, y lo que hay es una realización de largo finito. El atajo natural es saltear la autocorrelación y transformar la señal directo:

$$Ŝ(f)=\frac{1}{N}{|\sum _{n=0}^{N-1}x[n]e^{-j2\pi fn}|}^2$$

Eso es el periodograma, y es lo que hace wgn_psd.m en dos líneas. Es insesgado en el sentido que importa: su valor esperado es la PSD verdadera. Un estimador insesgado suena a estimador bueno, y ahí está la trampa de esta nota — insesgado dice dónde está el centro de la distribución, y no dice una palabra sobre cuán ancha es.

02 · Por qué no converge

Mirá una sola frecuencia. La suma de arriba es una suma de $N$ variables aleatorias independientes, cada una multiplicada por un número complejo de módulo uno. Por el teorema central del límite, la parte real y la parte imaginaria de esa suma son gaussianas de media cero, con la misma varianza, y prácticamente independientes entre sí. Entonces $Ŝ(f)$ es la suma de dos cuadrados de gaussianas: una exponencial. Y una exponencial tiene una propiedad que arruina el plan entero — su desvío es igual a su media. En símbolos, $E[Ŝ]=S(f)$ pero también $\mathrm{Var}[Ŝ]=S^2(f)$, y el cociente entre desvío y media vale uno.

En esa cuenta $N$ no aparece por ningún lado, y no aparece por una razón concreta: cada muestra nueva entra en la suma, pero también entra en el denominador y agrega una frecuencia más que estimar. Duplicar las muestras no le da el doble de datos a cada bin; le da el doble de bins, cada uno tan malo como antes. La resolución mejora —el espaciado entre frecuencias es $\Delta f=f_s/N$ y se achica— pero la incertidumbre de cada valor no se mueve. Es un estimador insesgado y no consistente, que es exactamente la combinación que la intuición no tiene prevista.

Las tres medias rondan el 1 teórico, como corresponde a un estimador insesgado. Las tres columnas de la derecha rondan el 1 también, y ésa es la mala noticia: el error típico de cada valor es del tamaño del valor. Un periodograma de cien mil muestras no es una PSD con poco ruido, es una PSD con cien mil puntos y el mismo ruido que antes. En pantalla se ve como un pasto muy fino y exactamente igual de alto.

03 · Lo que sí funciona: promediar tramos

Si el problema es que cada bin tiene una sola realización de una exponencial, la salida es tenerlas varias y promediarlas. Bartlett lo hace de la manera más directa posible: partir la serie en $K$ tramos de largo $L$, calcular el periodograma de cada uno y promediarlos punto a punto.

$${Ŝ}_B(f)=\frac{1}{K}\sum _{i=0}^{K-1}{Ŝ}_i(f)$$

Los tramos no se solapan, así que los $K$ periodogramas son independientes. Promediar $K$ cosas independientes divide la varianza por $K$ y deja la media donde estaba, de modo que el cociente entre desvío y media —que valía uno— pasa a valer $1/\sqrt{K}$. Esta vez sí converge, y converge a la velocidad de siempre. Welch es la misma idea con dos refinamientos que el apunte agrega después —ventanear cada tramo y solaparlos— y es lo que arma psd_with_fft_average.m; el mecanismo que baja la varianza es el promedio, en los dos casos.

La coincidencia con la ley es apretada hasta la tercera cifra, y no está ahí por casualidad: el 0,496 y el 0,500 salen de dos lugares distintos —uno de contar el desvío de cuatrocientas corridas, el otro de una raíz cuadrada— y coinciden porque el mecanismo es el que la teoría dice que es. Lo que queda es el canje: con un presupuesto fijo de $N=L·K$ muestras, la resolución $\Delta f=f_s/L$ y la incertidumbre $1/\sqrt{K}$ tiran para lados opuestos. Elegir un estimador de PSD es elegir dónde pararse en ese canje, y no existe la opción de no elegir.