El fenómeno de Gibbs

Cuando una función con un salto se arma con una cantidad finita de términos, la suma no llega al salto derecho: se pasa, vuelve y oscila a los dos lados. Eso es el fenómeno de Gibbs. En filtros aparece cada vez que se corta la respuesta al impulso de un filtro ideal, y deja un sobrepico que tiende a un 8,95 % del salto, sin importar cuántos coeficientes se usen.

01 · La definición

Sea una función con una discontinuidad de salto $J$, y una aproximación que la reconstruye con $N$ términos: una serie de Fourier truncada, o un filtro con la respuesta al impulso ideal cortada a $N$ muestras. El fenómeno de Gibbs es que, cerca del salto, la aproximación se pasa del valor de llegada y oscila alrededor de él. El sobrepico es cuánto se pasa, medido como fracción del salto:

$$\text{sobrepico}=\frac{\mathrm{máx}H-1}{J}$$

En el pasabajos el salto es el de la plantilla ideal en ${\omega }_c$, de 1 a 0, así que $J=1$ y el sobrepico es directamente cuánto se pasa $|H|$ de 1 en la banda de paso.

Las dos situaciones son la misma. La respuesta en frecuencia $H_D(\omega )$ de un filtro discreto es periódica de período $2\pi$, y la DTFT $H_D(\omega )={\sum }_n{h}_D[n]e^{-j\omega n}$ es su serie de Fourier: los coeficientes de la serie son las muestras $h_D[n]$. Quedarse con $N$ muestras es quedarse con $N$ términos de la serie de un rectángulo, que es el caso clásico en el que se descubrió el fenómeno.

02 · De dónde sale en un filtro

El punto de partida es el de la filmina del método de ventana: el pasabajos ideal vale 1 entre $-{\omega }_c$ y ${\omega }_c$ y 0 en el resto del período. Su respuesta al impulso sale de la DTFT inversa:

$$h_D[n]=\frac{1}{2\pi }{\int }_{-{\omega }_c}^{{\omega }_c}{e}^{j\omega n}d\omega =\frac{\mathrm{sen}({\omega }_{c}n)}{\pi n},h_D[0]=\frac{{\omega }_c}{\pi }$$

Es la $2f_c\mathrm{sen}(n{\omega }_c)/(n{\omega }_c)$ de la filmina, con ${\omega }_c=2\pi f_c$, y el valor en $n=0$ sale por L'Hôpital. Es infinita hacia los dos lados. Para usarla hay que cortarla, y el corte es lo que produce el rizado.

Paso 1: truncar es multiplicar por una ventana. Quedarse con las $N$ muestras centrales (con $N$ impar) es multiplicar $h_D$ por una ventana $w[n]$ que vale 1 para $|n|\le (N-1)/2$ y 0 afuera. Es la ventana rectangular: no usar ventana es usar la más brusca de todas.

$$h[n]=h_D[n]w[n]$$

Después se corre $(N-1)/2$ muestras para que sea causal. El corrimiento sólo agrega una fase lineal y no toca el módulo, así que todo lo que sigue se hace con la versión centrada, cuya $H(\omega )$ es real.

Paso 2: multiplicar en el tiempo es convolucionar en frecuencia. Con $W(\omega )$ la DTFT de la ventana, la convolución es circular, sobre un período:

$$H(\omega )=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{H}_D(\theta )W(\omega -\theta )d\theta$$

Es el panel (c) de la filmina: $|H(\omega )|$ sale de pasar $W(\omega )$ por encima de $H_D(\omega )$. Para la rectangular centrada, $W$ es una serie geométrica que se suma igual que la de la media móvil:

$$W_R(\omega )=\sum _{n=-(N-1)/2}^{(N-1)/2}{e}^{-j\omega n}=\frac{\mathrm{sen}(\omega N/2)}{\mathrm{sen}(\omega /2)}$$

Es el núcleo de Dirichlet, derivado paso a paso en la media móvil: la misma sinc periódica, sin el $1/N$ y sin la fase porque acá está centrada. Tiene un lóbulo principal de alto $N$ con ceros en $\pm 2\pi /N$, y lóbulos laterales que alternan de signo y cuyo primero está a −13,26 dB.

Paso 3: contra un rectángulo, la convolución es una integral entre dos bordes. $H_D$ vale 1 en $[-{\omega }_c,{\omega }_c]$ y 0 afuera, así que sólo recorta el intervalo de integración. Cambiando la variable a $\phi =\omega -\theta$:

$$H(\omega )=\frac{1}{2\pi }{\int }_{\omega -{\omega }_c}^{\omega +{\omega }_c}W(\phi )d\phi$$

$H(\omega )$ es el área de $W$ que cae adentro de una ventana de ancho $2{\omega }_c$. Contra freqz, con $N=21$ y ${\omega }_c=\mathrm{0,4}\pi$, esta integral hecha con integral da lo mismo en cinco frecuencias con una diferencia de 7×10⁻¹⁶ como máximo.

Paso 4: cerca del borde sólo importa el borde. La integral de $W_R$ sobre un período entero vale $2\pi$, porque de toda la suma sólo sobrevive el término $n=0$. Cuando $\omega$ está cerca de ${\omega }_c$, el límite de arriba, $\omega +{\omega }_c$, cae lejos del lóbulo principal, donde $W$ ya es una cola chica que alterna. Despreciándola, y usando que $W$ es par:

$$H(\omega )\approx \frac{1}{2\pi }{\int }_{-\pi }^{{\omega }_c-\omega }W(\phi )d\phi$$

Cerca del corte, $H$ es la integral acumulada de la transformada de la ventana, hasta la distancia que queda al borde. Lejos adentro de la banda de paso esa integral ya juntó casi toda el área y vale casi 1; justo en ${\omega }_c$ juntó la mitad; lejos afuera, casi nada. Cada lóbulo lateral de $W$ que entra o sale suma o resta un poco, y eso es el rizado.

03 · El valor límite: la integral seno

Paso 1: con N grande el núcleo es una sinc. El lóbulo principal y los primeros laterales viven en $|\phi |$ del orden de $2\pi /N$, donde $\mathrm{sen}(\phi /2)\approx \phi /2$:

$$W_R(\phi )\approx \frac{\mathrm{sen}(N\phi /2)}{\phi /2}$$

Paso 2: cambiar a la variable del lóbulo. Con $u=N\phi /2$ es $d\phi =2du/N$ y $\phi /2=u/N$. Las $N$ se cancelan todas:

$$\frac{1}{2\pi }·\frac{\mathrm{sen}u}{u/N}·\frac{2du}{N}=\frac{1}{\pi }\frac{\mathrm{sen}u}{u}du$$

Paso 3: la integral acumulada es la integral seno. El límite de abajo, $-\pi$, pasa a $-N\pi /2$, que para $N$ grande es $-\infty$. Con $\mathrm{Si}(x)={\int }_0^x\mathrm{sen}u/udu$ y $\mathrm{Si}(\infty )=\pi /2$:

$$H(\omega )\approx \frac{1}{\pi }{\int }_{-\infty }^x\frac{\mathrm{sen}u}{u}du=\frac{1}{2}+\frac{\mathrm{Si}(x)}{\pi },x=\frac{N({\omega }_c-\omega )}{2}$$

Paso 4: el máximo cae en el primer cero. La derivada respecto de $x$ es $\mathrm{sen}x/(\pi x)$, que se anula en $x=k\pi$. De todos esos extremos el mayor es el primero, $x=\pi$: hasta ahí se sumó todo el lóbulo principal y nada del primer lóbulo negativo. Vale

$$H_{\mathrm{máx}}=\frac{1}{2}+\frac{\mathrm{Si}(\pi )}{\pi }=\frac{1}{2}+\frac{\mathrm{1,8519}}{\pi }\approx \mathrm{1,0895}$$

Un sobrepico de 8,95 % del salto. $\mathrm{Si}(\pi )$ no tiene forma más simple, pero se calcula a cualquier precisión: integral(@(x) sin(x)./x, 0, pi) da 1,8519, y con él 1,0895. El 9 % que dice el apunte es este número redondeado.

De la misma cuenta salen dos cosas más. Dónde está el pico: $x=\pi$ es ${\omega }_c-\omega =2\pi /N$, una mitad de lóbulo principal antes del corte. Cuánto rechaza la banda suprimida: del otro lado del corte, el extremo más grande es $x=-\pi$, y como $\mathrm{Si}$ es impar vale $\frac{1}{2}-\frac{\mathrm{Si}(\pi )}{\pi }=-\mathrm{0,0895}$. El primer lóbulo de la banda suprimida tiene exactamente la altura del sobrepico:

$$20{\log }_{10}\mathrm{0,0895}\approx -\mathrm{20,96}\text{ dB}$$

Los −21 dB de la ventana rectangular no son un dato de tabla: son el mismo 8,95 % visto del otro lado del salto. Y no son los −13 dB del primer lóbulo lateral de $W_R$. El filtro no copia la altura de los lóbulos de la ventana: acumula su área, y el área de un lóbulo lateral es mucho más chica que su pico.

04 · Por qué no baja con N

En la forma de la sección anterior, $N$ sólo aparece adentro de $x=N({\omega }_c-\omega )/2$. Cambiar $N$ estira o comprime el eje de frecuencia, pero la curva $\frac{1}{2}+\mathrm{Si}(x)/\pi$ tiene siempre la misma altura. Más taps no bajan el pico: lo acercan al corte, como $2\pi /N$. La tesis entera, con su propio widget, está en el ensayo sobre el rizado. Acá van los números. Pasabajos de ${\omega }_c=\mathrm{0,4}\pi$ con ventana rectangular, medido en MATLAB sobre 262 145 frecuencias entre 0 y $\pi$:

La columna del medio es la distancia del pico al corte, medida en mitades de lóbulo principal: se clava en 1, como predice el paso 4. La primera y la última columna se van a los valores de la integral seno. Pero no bajan de a poco: oscilan. De $N=13$ a $N=25$, los dos casos de la filmina, el sobrepico sube de 6,83 % a 10,38 %. Barriendo todos los $N$ impares de un tramo:

Con $N$ finito el paso 4 no es exacto: la cola de $W$ en el otro borde, a $2{\omega }_c$ de distancia, suma o resta un poco según dónde caiga el corte respecto de los ceros de $h_D$. Con ${\omega }_c=\mathrm{0,4}\pi$ esos ceros se repiten cada 2,5 muestras, y el patrón se repite cada 10 valores de $N$. Un caso extremo lo muestra: $N=19$ y $N=21$ dan el mismo 9,67 %, porque las dos muestras que agrega el 21 son $h_D[\pm 10]=\mathrm{sen}(4\pi )/(10\pi )=0$: es el mismo filtro. La franja se cierra sobre 8,95 % al crecer $N$, pero no se cierra hacia cero.

05 · Lo que sí lo baja: la ventana

Si el rizado es la integral acumulada de $W$, para bajarlo hay que cambiar $W$. El método de la filmina es el mismo $h[n]=h_D[n]w[n]$, con una ventana que en vez de cortar de golpe baja hasta casi cero en los bordes. Las dos de la cátedra, escritas para $n=0,\dots ,N-1$ como las genera MATLAB:

$$w_{\text{Hamming}}[n]=\mathrm{0,54}-\mathrm{0,46}\cos \frac{2\pi n}{N-1}$$ $$w_{\text{Blackman}}[n]=\mathrm{0,42}-\mathrm{0,5}\cos \frac{2\pi n}{N-1}+\mathrm{0,08}\cos \frac{4\pi n}{N-1}$$

Un corte suave tiene una transformada con lóbulos laterales mucho más bajos, que es lo que baja el rizado. Pero el lóbulo principal se ensancha: con $N=1001$ su primer cero cae en 1,00, 2,00 y 3,00 veces $2\pi /N$ para la rectangular, Hamming y Blackman. La transición de $H$ dura lo que tarda ese lóbulo en cruzar el borde, así que se ensancha con él. Las filminas 28 y 29 lo muestran con $N=57$. Medido sobre el pasabajos, con ${\omega }_c=\mathrm{0,4}\pi$ y el mismo $N=41$ para las tres:

La transición se midió como la distancia entre el último punto de la banda de paso con $|H|\ge 1-\delta$ y el primero de la suprimida con $|H|\le \delta$, con $\delta$ el peor lóbulo de la banda suprimida. En unidades de $2\pi /N$ es casi constante: con $N=161$ da 0,93, 3,35 y 5,62. La banda suprimida, en cambio, casi no se mueve con $N$ para ninguna de las tres: −53,61 dB con 21 taps y −53,97 dB con 161 para Hamming, −74,36 y −75,29 dB para Blackman. El canje queda así: la ventana decide cuánto rechaza el filtro, y $N$ decide qué tan angosta es la transición. Por eso con Blackman, para la misma transición que la rectangular, hacen falta unas seis veces más taps. Cómo se usa eso para cumplir una máscara es tema del diseño del FIR.

La última columna es el peor lóbulo lateral de la ventana sola, que es el número de las tablas. El apunte da −13, −41 y −57 dB. El de Hamming depende de $N$: medido, va de −40,64 dB con $N=21$ a −42,67 dB con $N=1001$. El de Blackman se queda entre −58,1 y −58,3 dB. En los tres casos el filtro rechaza bastante más de lo que dice ese lóbulo, por la razón de la sección 03: el filtro acumula áreas, no alturas. Y ojo al medirlo: con Hamming y Blackman el primer lóbulo lateral no es el peor, así que hay que tomar el máximo de toda la banda.

Lo que conviene mirar es el borde derecho del rectángulo, en ${\omega }_c$, mientras $W$ se le acerca. Con la rectangular, el máximo de abajo llega cuando el primer cero de $W$ toca el borde, con la lectura del borde en 1: todo el lóbulo principal adentro, el primer lóbulo negativo afuera. Un poco más y el lóbulo negativo empieza a entrar y $H$ baja. Con otro $N$ todo pasa igual, más apretado contra ${\omega }_c$. Con Hamming o Blackman y el mismo $N$, los lóbulos laterales casi no se ven, la curva de abajo llega a 1 sin pasarse visiblemente, y el lóbulo principal tarda mucho más en cruzar el borde. Al llegar a $\omega =\pi$ las tres lecturas del recorrido dan los números de la tabla.

06 · El análisis en MATLAB

Las mediciones de la nota, con fir1 y las tres ventanas. Lo que figura como comentario es la salida real de correrlo en R2026a.

%% 1. El valor límite, por la integral seno
Si = integral(@(x) sin(x)./x, 0, pi);
fprintf('Si(pi) = %.4f   1/2 + Si(pi)/pi = %.4f\n', Si, 0.5 + Si/pi)
fprintf('lóbulo en banda suprimida: %.2f dB\n', 20*log10(Si/pi - 0.5))
% Si(pi) = 1.8519   1/2 + Si(pi)/pi = 1.0895
% lóbulo en banda suprimida: -20.96 dB

%% 2. Sobrepico, atenuación y transición, para tres ventanas y cuatro N
wc = 0.4;                                  % corte, en unidades de pi
w  = linspace(0, pi, 2^18 + 1);            % vector: nunca un escalar
ventanas = {@rectwin, @hamming, @blackman};
for iv = 1:3
    fprintf('%s\n', func2str(ventanas{iv}))
    for N = [21 41 81 161]
        h  = fir1(N-1, wc, ventanas{iv}(N), 'noscale');   % orden N-1: N taps
        H  = abs(freqz(h, 1, w));
        sobrepico = (max(H(w <= wc*pi)) - 1) * 100;
        Hs = H(w > wc*pi);  ws_ = w(w > wc*pi);
        k  = find(diff(Hs) > 0, 1);             % primer mínimo: ahí empieza el rizado
        d  = max(Hs(k:end));                    % el peor lóbulo, no el primero
        wp = w(find(H >= 1 - d & w <= wc*pi, 1, 'last'));
        ws = ws_(find(Hs <= d, 1));
        fprintf('  N = %3d   sobrepico %7.4f %%   As %7.2f dB   transición %.2f x 2pi/N\n', ...
                N, sobrepico, 20*log10(d), (ws - wp) * N / (2*pi))
    end
end
% rectwin
%   N =  21   sobrepico  9.6716 %   As  -21.30 dB   transición 0.98 x 2pi/N
%   N =  41   sobrepico  9.2573 %   As  -21.17 dB   transición 0.95 x 2pi/N
%   N =  81   sobrepico  9.0905 %   As  -21.08 dB   transición 0.94 x 2pi/N
%   N = 161   sobrepico  9.0167 %   As  -21.02 dB   transición 0.93 x 2pi/N
% hamming
%   N =  21   sobrepico  0.1643 %   As  -53.61 dB   transición 3.50 x 2pi/N
%   N =  41   sobrepico  0.1680 %   As  -53.61 dB   transición 3.41 x 2pi/N
%   N =  81   sobrepico  0.1834 %   As  -53.81 dB   transición 3.37 x 2pi/N
%   N = 161   sobrepico  0.1900 %   As  -53.97 dB   transición 3.35 x 2pi/N
% blackman
%   N =  21   sobrepico  0.0338 %   As  -74.36 dB   transición 6.28 x 2pi/N
%   N =  41   sobrepico  0.0181 %   As  -75.30 dB   transición 5.73 x 2pi/N
%   N =  81   sobrepico  0.0173 %   As  -75.30 dB   transición 5.66 x 2pi/N
%   N = 161   sobrepico  0.0172 %   As  -75.29 dB   transición 5.62 x 2pi/N

%% 3. El sobrepico con la rectangular no baja: oscila y se cierra sobre 8,95 %
for tramo = {11:2:31, 101:2:121, 1001:2:1021}
    s = zeros(size(tramo{1}));
    for i = 1:numel(tramo{1})
        N = tramo{1}(i);
        H = abs(freqz(fir1(N-1, wc, rectwin(N), 'noscale'), 1, w(w <= wc*pi)));
        s(i) = (max(H) - 1) * 100;
    end
    fprintf('N de %4d a %4d: sobrepico entre %.2f %% y %.2f %%\n', ...
            tramo{1}(1), tramo{1}(end), min(s), max(s))
end
% N de   11 a   31: sobrepico entre 6.83 % y 11.46 %
% N de  101 a  121: sobrepico entre 8.69 % y 9.27 %
% N de 1001 a 1021: sobrepico entre 8.92 % y 8.98 %

%% 4. La trampa de la escala: fir1 por omisión normaliza la ganancia en continua
h = fir1(20, wc, rectwin(21));             % sin 'noscale'
H = abs(freqz(h, 1, w));
fprintf('N = 21 con escala: sobrepico %.2f %%, suma de h = %.4f\n', ...
        (max(H(w <= wc*pi)) - 1) * 100, sum(fir1(20, wc, rectwin(21), 'noscale')))
% N = 21 con escala: sobrepico 14.63 %, suma de h = 0.9568

Cuatro detalles que muerden si no se los conoce:

• fir1 pide el orden, no los taps. fir1(N-1, ...) da $N$ coeficientes, y la ventana que se le pasa tiene que medir $N$: hamming(N), no hamming(N-1). Y el corte va normalizado a $\pi$: el 0.4 es $\mathrm{0,4}\pi$.
• Por omisión, fir1 escala el filtro. Divide $h$ por su suma para que la ganancia en continua dé exactamente 1. Con la rectangular eso cambia el sobrepico, porque en $\omega =0$ el rizado no vale 1: con $N=21$ la suma es 0,9568, y el sobrepico escalado da 14,63 % en vez de 9,67 %. Para ver el fenómeno como es, 'noscale'. Con Hamming y Blackman la suma ya es casi 1 (0,9980 y 1,0003 con $N=21$) y el sobrepico cambia dos décimas de punto o menos.
• La atenuación es el peor lóbulo, no el primero. Con Hamming y Blackman los lóbulos de la banda suprimida no bajan en orden, así que el bloque 2 toma el máximo de toda la banda después del primer mínimo.
• La grilla importa. El pico de Gibbs mide $2\pi /N$ de ancho. Con freqz(h, 1, 512) y $N=1001$ cae un solo punto en el lóbulo del pico, y el sobrepico sale 6,26 % en vez de 8,96 %; con $N=161$, 8,83 % en vez de 9,02 %. Por eso el vector w de $2^{18}+1$ puntos, que además es un vector: un escalar como tercer argumento freqz lo lee como cantidad de puntos, no como frecuencia.

Para tener a mano

1. Truncar $h_D[n]=\mathrm{sen}({\omega }_{c}n)/(\pi n)$ es multiplicar por una ventana, y eso es convolucionar $H_D$ con $W$. Con la rectangular, $W$ es el núcleo de Dirichlet.
2. Contra el rectángulo ideal, $H(\omega )$ es el área de $W$ adentro de $[\omega -{\omega }_c,\omega +{\omega }_c]$, dividida por $2\pi$. Cerca del corte, la integral acumulada de $W$.
3. El sobrepico tiende a $\mathrm{Si}(\pi )/\pi -\frac{1}{2}\approx \mathrm{0,0895}$: 8,95 % del salto, a $2\pi /N$ antes del corte.
4. El primer lóbulo de la banda suprimida tiene la misma altura: −20,96 dB. No es el −13 dB de la ventana: el filtro acumula áreas, no alturas.
5. Más taps no bajan el sobrepico: lo comprimen contra el corte. Con $N$ chico oscila, de 6,83 % a 11,46 % entre 11 y 31 taps.
6. Lo baja la ventana, a costa de transición. Con $N=41$: rectangular −21 dB y 0,95 · 2π/N; Hamming −54 dB y 3,4 · 2π/N; Blackman −75 dB y 5,7 · 2π/N.
7. En MATLAB: fir1(N-1, wc/pi, ventana(N), 'noscale') para ver el fenómeno sin la normalización.