El filtro de media móvil

La media móvil reemplaza cada muestra por el promedio de las últimas $M$. Todo lo demás sale de ahí: su respuesta al impulso es un rectángulo, la DTFT de un rectángulo es una sinc periódica, y el largo $M$ de la ventana decide dónde caen los ceros de esa sinc y qué tan angosto es el lóbulo que deja pasar.

01 · La definición

La ecuación en diferencias, en su versión causal:

$$y[n]=\frac{1}{M}\sum _{k=0}^{M-1}x[n-k]$$

Es el promediador del apunte, $\frac{1}{M_1+M_2+1}{\sum }_{k=-M_1}^{M_2}x[n-k]$, con $M_1=0$ —no mira muestras futuras, así que es causal— y $M_2=M-1$. Es lineal, porque es una suma con pesos fijos, e invariante en el tiempo, porque la ventana tiene la misma forma en cada $n$. Por ser LTI queda descripto entero por su respuesta al impulso.

02 · La respuesta al impulso

Cómo se obtiene: se pone un impulso a la entrada, $x[n]=\delta [n]$, y lo que sale es $h[n]$ por definición:

$$h[n]=\frac{1}{M}\sum _{k=0}^{M-1}\delta [n-k]=\{\begin{array}{ll}\frac{1}{M}& 0\le n\le M-1\\ 0& \text{en otro caso}\end{array}$$

Cada término de la suma es un impulso corrido: $\delta [n-k]$ vale 1 sólo en $n=k$. Entre los $M$ términos cubren $n=0,\dots ,M-1$ una vez cada uno, y el $1/M$ les pone la altura. Queda un rectángulo de $M$ muestras de alto $1/M$. Tres cosas se leen directo de él:

• Es FIR: dura $M$ muestras y después es cero. No hay realimentación, así que es estable siempre.
• Suma 1: ${\sum }_{n}h[n]=M·\frac{1}{M}=1$. El promedio de una constante es la misma constante, así que la ganancia en continua es exactamente 1.
• Es simétrico alrededor de $n=(M-1)/2$. Eso es lo que después da fase lineal (sección 06).

En MATLAB se obtiene igual que en papel: filter(ones(1,8)/8, 1, [1 zeros(1,15)]) devuelve 0.1250 ocho veces y después ceros. Es el rectángulo, medido.

03 · La DTFT, paso a paso

La DTFT de $h[n]$ es la suma de siempre. Como $h$ sólo es distinta de cero entre $0$ y $M-1$, la suma infinita se corta ahí:

$$H(e^{j\omega })=\sum _{n=-\infty }^{\infty }h[n]e^{-j\omega n}=\frac{1}{M}\sum _{n=0}^{M-1}{\left(e^{-j\omega }\right)}^n$$

Paso 1: es una serie geométrica. Con razón $r=e^{-j\omega }$ y $\sum _{n=0}^{M-1}{r}^n=\frac{1-r^M}{1-r}$, que vale mientras $r\ne 1$, o sea para $\omega \ne 0$:

$$H(e^{j\omega })=\frac{1}{M}·\frac{1-e^{-j\omega M}}{1-e^{-j\omega }}$$

Paso 2: sacar el ángulo mitad. Es el truco que convierte una resta de exponenciales en un seno. Para cualquier $\theta$:

$$1-e^{-j\theta }=e^{-j\theta /2}(e^{j\theta /2}-e^{-j\theta /2})=e^{-j\theta /2}·2j\mathrm{sen}(\theta /2)$$

Paso 3: aplicarlo arriba y abajo. Con $\theta =\omega M$ en el numerador y $\theta =\omega$ en el denominador, los $2j$ se cancelan y las dos exponenciales se juntan en una sola, $e^{-j\omega M/2}/e^{-j\omega /2}=e^{-j\omega (M-1)/2}$:

$$H(e^{j\omega })=\underset{\text{fase lineal}}{\underbrace{e^{-j\omega (M-1)/2}}}·\underset{\text{amplitud real},D_M(\omega )}{\underbrace{\frac{\mathrm{sen}(\omega M/2)}{M\mathrm{sen}(\omega /2)}}}$$

En $\omega =0$ la fórmula da $0/0$, pero el límite es 1, que es la suma de $h[n]$: la ganancia en continua de la sección anterior. Contra freqz, sobre 4096 frecuencias del círculo entero y con $M=8$, la forma cerrada difiere en 3,4×10⁻¹⁵ como máximo: es la misma función.

La fórmula tiene dos factores con papeles distintos. El primero tiene módulo 1 siempre: no toca la amplitud de nada, sólo atrasa. El segundo, $D_M(\omega )$, es real y cambia de signo. Así que $|H|=|D_M|$, y todo lo que el filtro hace con las amplitudes está en ese cociente de senos.

04 · Una sinc, pero periódica

A $D_M(\omega )$ se le dice sinc y se parece a una, pero no es: es el núcleo de Dirichlet, o sinc periódica. La diferencia está en el denominador. Cerca de $\omega =0$, $\mathrm{sen}(\omega /2)\approx \omega /2$, y queda la sinc de verdad:

$$D_M(\omega )\approx \frac{\mathrm{sen}(\omega M/2)}{\omega M/2}\text{para }\omega \text{ chico}$$

Lejos del origen la aproximación se rompe, y se rompe por la razón de siempre: $h[n]$ es discreta, así que $H(e^{j\omega })$ es periódica de período $2\pi$. La sinc de verdad decae para siempre; la de Dirichlet vuelve a subir y repite su lóbulo principal en cada múltiplo de $2\pi$. Evaluada con freqz en $\omega$ y en $\omega +2\pi$, la diferencia es de 9×10⁻¹⁵. Medida la distancia entre las dos curvas:

Diferencia máxima entre $D_M(\omega )$ y $\mathrm{sen}(\omega M/2)/(\omega M/2)$, en MATLAB sobre 4000 y 40 000 puntos. En el lóbulo principal el error cae como $1/M^2$: con $M$ grande, la parte que importa es una sinc. En las colas cae sólo como $1/M$, y en $\omega =\pi$ las dos curvas se separan por una cuenta exacta: $D_M(\pi )=\mathrm{sen}(M\pi /2)/M$, que es $0$ para $M$ par y $\pm 1/M$ para $M$ impar.

05 · Cómo la cambia la ventana

Todo lo que importa de $|H|$ depende de un solo número, $M$, y casi todo escala como $1/M$.

Los ceros. El numerador se anula cuando $\omega M/2=k\pi$, o sea en

$${\omega }_k=\frac{2\pi k}{M},k=1,\dots ,M-1$$

El $k=0$ no cuenta porque ahí también se anula el denominador. Son $M-1$ ceros por período, equiespaciados, y son ceros exactos: freqz evaluado en ellos, con $M=8$, da $|H|$ de 3,3×10⁻¹⁶ o menos. Tienen una lectura directa: una sinusoide de frecuencia $2\pi k/M$ completa $k$ ciclos enteros dentro de la ventana, y el promedio de ciclos enteros es cero.

El lóbulo principal va del primer cero negativo al primero positivo, de $-2\pi /M$ a $2\pi /M$: mide $4\pi /M$. Duplicar la ventana lo parte a la mitad.

El corte a −3 dB no tiene forma cerrada: hay que despejar $|D_M(\omega )|=1/\sqrt{2}$ numéricamente (con fzero). Medido como fracción del primer cero, ${\omega }_{3\mathrm{dB}}/(2\pi /M)$ da 0,4554 con $M=4$, 0,4460 con 8, 0,4437 con 16 y 0,4432 con 30, y se clava en 0,4430 para $M$ grande. O sea, ${\omega }_{3\mathrm{dB}}\approx \mathrm{2,78}/M$. Con el $M=30$ del script de la cátedra da 0,0928 rad.

El primer lóbulo lateral es el pico entre el primer cero y el segundo. Su altura en dB respecto del lóbulo principal, y dónde cae:

Los −13 dB de la tabla de ventanas del apunte son el límite de esta fila, no un valor fijo: con ventanas cortas el lóbulo lateral está más alto. Y al revés de lo que pasa con el ancho, no baja agregando muestras, porque para $M$ grande la forma del lóbulo ya no cambia, sólo se angosta. Es el mismo núcleo que aparece al truncar una respuesta ideal con ventana rectangular: la ventana rectangular es este mismo $h[n]$ sin el $1/M$, y por eso su rizado tampoco baja con más taps.

Al subir $M$ pasan tres cosas juntas y ninguna se puede separar de las otras. Arriba el rectángulo se ensancha y baja, porque su área sigue siendo 1. Al medio el lóbulo principal se angosta y aparecen más ceros, siempre equiespaciados cada $2\pi /M$. Abajo la pendiente de la fase se empina, porque el retardo crece. Con la sinc no periódica superpuesta y $M$ chico se ve dónde se despegan: cerca de $\pm \pi$, y del todo pasado $\pi$, donde la de Dirichlet vuelve a subir para repetir su lóbulo en $\pm 2\pi$.

06 · La fase

De la forma cerrada, la fase es $-\omega (M-1)/2$, una recta, más un salto de $\pi$ en cada cero, donde $D_M$ cambia de signo. Esos saltos no atrasan nada: son un $-1$ que multiplica a una frecuencia que el filtro igual anula. Por eso, con $M$ grande, la fase del gráfico parece un diente de sierra casi plano: entre dos ceros la recta baja $\pi (M-1)/M$, casi $\pi$, y el salto lo repone. Es lo mismo que devuelve angle en MATLAB. Lo que cuenta es la pendiente de cada tramo, y el retardo de grupo es menos la pendiente:

$$\tau (\omega )=-\frac{d\angle H}{d\omega }=\frac{M-1}{2}$$

Es constante: todas las frecuencias salen atrasadas lo mismo, y la forma de la señal no se deforma. grpdelay con $M=8$ da 3,5 muestras en las cuatro frecuencias en que se lo pregunta. Cuánto cuesta ese retardo, y por qué no hay manera de suavizar sin pagarlo, está en la nota del promediador.

07 · El análisis en MATLAB

El análisis completo, con $M=8$. Lo que figura como comentario es la salida real de correrlo en R2026a. El bloque 8 es el barrido que hace el script Ejemplo_filtros.m de la clase.

%% 1. El filtro: M coeficientes iguales a 1/M
M = 8;
h = ones(1, M) / M;

%% 2. La respuesta al impulso, obtenida: se le mete un delta
d = [1, zeros(1, 15)];
h_medida = filter(h, 1, d);
figure; stem(0:15, h_medida, 'filled'); grid on
xlabel('n'); ylabel('h[n]')
% 0.1250 ... 0.1250 (ocho veces), después 0

%% 3. La DTFT evaluada por MATLAB
Nw = 4096;
[H, w] = freqz(h, 1, Nw, 'whole');      % w en [0, 2*pi)
w = w - 2*pi*(w >= pi);                 % pasar a [-pi, pi)
[w, orden] = sort(w);  H = H(orden);

%% 4. La misma DTFT, por la forma cerrada
Hc = exp(-1j*w*(M-1)/2) .* sin(w*M/2) ./ (M*sin(w/2));
Hc(w == 0) = 1;                         % el límite 0/0 en w = 0
fprintf('max |freqz - cerrada| = %.1e\n', max(abs(H - Hc)))
% max |freqz - cerrada| = 3.4e-15

%% 5. Módulo, en veces y en dB, y fase
figure
subplot(311); plot(w/pi, abs(H)); grid on; ylabel('|H|')
subplot(312); plot(w/pi, 20*log10(abs(H))); grid on; ylabel('|H| [dB]'); ylim([-50 3])
subplot(313); plot(w/pi, angle(H)); grid on; ylabel('\angle H [rad]'); xlabel('\omega / \pi')

%% 6. Los ceros caen en 2*pi*k/M
wk = 2*pi*(1:M-1)/M;
fprintf('max |H| en los ceros = %.1e\n', max(abs(freqz(h, 1, wk))))
% max |H| en los ceros = 3.3e-16

%% 7. Retardo de grupo: constante e igual a (M-1)/2
disp(grpdelay(h, 1, [0.1 0.5 1.5 2.5]).')
% 3.5000    3.5000    3.5000    3.5000

%% 8. El barrido de M, como en el script de la cátedra
figure; hold on
for M = [4 8 16 32]
    [H, w] = freqz(ones(1, M)/M, 1, 8192);
    plot(w/pi, 20*log10(abs(H)), 'DisplayName', sprintf('M = %d', M))
end
yline(-3, '--k', 'HandleVisibility', 'off')
grid on; ylim([-50 3]); legend; xlabel('\omega / \pi'); ylabel('|H| [dB]')

Tres detalles que muerden si no se los conoce:

• Un escalar en freqz no es una frecuencia. freqz(h, 1, 0.3) no evalúa en $\omega =\mathrm{0,3}$: lee el escalar como la cantidad de puntos, así que con 0,3 da error y con 3 devuelve tres puntos cualesquiera sin avisar. grpdelay(h, 1, 0.3) es peor, porque devuelve un vector vacío sin error. Con $M=2$ hay un solo cero, así que el bloque 6 le pasaría un escalar. Hay que pasarles al menos dos frecuencias.
• Por omisión freqz devuelve sólo $[0,\pi )$. Para $h$ real alcanza, porque $|H|$ es par. Para ver el lóbulo principal entero, centrado, va 'whole' y el corrimiento del bloque 3.
• La forma cerrada da NaN en $\omega =0$, que es $0/0$ y no el límite. Hay que poner el 1 a mano. Y en dB, los ceros exactos dan -Inf: sin el ylim el gráfico se aplasta.

Para tener a mano

1. $h[n]=1/M$ para $0\le n\le M-1$: un rectángulo de área 1. Se obtiene metiendo un $\delta [n]$.
2. La DTFT sale de la serie geométrica más el ángulo mitad:$$H(e^{j\omega })=e^{-j\omega (M-1)/2}·\mathrm{sen}(\omega M/2)/(M\mathrm{sen}(\omega /2))$$
3. Es una sinc periódica (núcleo de Dirichlet): sinc cerca del origen, periódica de $2\pi$ como toda DTFT.
4. Ceros en $2\pi k/M$. Lóbulo principal de $4\pi /M$. Corte a −3 dB en $\approx \mathrm{2,78}/M$. Todo escala como $1/M$.
5. El primer lóbulo lateral tiende a −13,26 dB y no baja con $M$. Con ventanas cortas está más alto: −11,30 dB con $M=4$.
6. La fase es lineal, con saltos de $\pi$ en los ceros. El retardo de grupo es constante, $(M-1)/2$.