El filtro FIR

Un FIR calcula cada muestra de salida con las últimas $N$ muestras de entrada y nada más: no mira salidas anteriores. De esa sola restricción sale todo lo demás. Su respuesta al impulso dura $N$ muestras, así que no puede volverse inestable. Sus coeficientes son esa misma respuesta al impulso, así que se los puede elegir simétricos, y un FIR simétrico tiene fase exactamente lineal. El precio es el largo: una transición angosta pide muchos coeficientes, y cada uno es una multiplicación por muestra.

01 · La definición

La ecuación en diferencias de un FIR de $N$ coeficientes, como está en la filmina:

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

Es la ecuación en diferencias general del apunte sin la parte recursiva: ningún término $y[n-k]$ del lado derecho. Por eso en MATLAB se filtra con filter(h, 1, x): el denominador es 1.

Los coeficientes son la respuesta al impulso. Con $x[n]=\delta [n]$, cada término de la suma vale $h[k]\delta [n-k]$, que es distinto de cero sólo en $n=k$. Así que la salida en el instante $n$ es el coeficiente número $n$. No hay nada que calcular para pasar de la ecuación a $h[n]$: se lee.

Dura $N$ muestras. Para $n\ge N$ ningún término de la suma alcanza al impulso, y $h[n]=0$. De ahí el nombre, respuesta al impulso finita. A cada coeficiente se le dice tap, derivación, por el hardware que lo implementa: una línea de retardos con una derivación por coeficiente, que es el FIR transversal. «Un FIR de 41 taps» es un FIR con $N=41$.

El ejemplo más simple es la media móvil: $N$ coeficientes iguales a $1/N$. El del script de la clase, Ejemplo_filtros.m, tiene siete: [0.1 0.2 0.3 0.4 0.3 0.2 0.1], dividido por su suma para que la ganancia en continua sea 1.

02 · H(z) es un polinomio en z⁻¹

La transformada Z de $h[n]$ tiene tantos términos como coeficientes:

$$H(z)=\sum _{k=0}^{N-1}h[k]z^{-k}=h[0]+h[1]z^{-1}+\cdots +h[N-1]z^{-(N-1)}$$

No hay denominador: es un polinomio de grado $N-1$ en $z^{-1}$. Multiplicando y dividiendo por $z^{N-1}$ se ve qué raíces tiene:

$$H(z)=\frac{h[0]z^{N-1}+h[1]z^{N-2}+\cdots +h[N-1]}{z^{N-1}}$$

Arriba, $N-1$ ceros, que pueden caer en cualquier lugar del plano. Abajo, $N-1$ polos, todos en $z=0$. Todo lo que el filtro hace con las frecuencias lo deciden los ceros; los polos del origen sólo aportan retardo. Cómo leer esos ceros en el plano, y qué le hacen a $|H|$ según lo cerca que estén del círculo unitario, es de la nota de polos y ceros. La respuesta en frecuencia es este mismo polinomio evaluado en $z=e^{j\omega }$.

03 · Siempre estable

Un sistema LTI es estable BIBO si su respuesta al impulso es absolutamente sumable. Para un FIR la suma tiene $N$ términos finitos, así que es finita siempre, sean cuales sean los coeficientes. La cota sale directo de la ecuación: si la entrada cumple $|x[n]|\le B$,

$$|y[n]|\le \sum _{k=0}^{N-1}|h[k]||x[n-k]|\le B\sum _{k=0}^{N-1}|h[k]|<\infty$$

Visto en el plano Z es lo mismo: la condición de estabilidad pide todos los polos adentro del círculo unitario, y los del FIR están en el origen. No hay coeficiente que los mueva. Lo que se pierde por no realimentar es la memoria larga barata: un IIR de un polo, $y[n]=x[n]+\alpha y[n-1]$, recuerda para siempre con una multiplicación, y a cambio puede explotar. Eso es de la nota del filtro IIR.

04 · Fase lineal por simetría, paso a paso

La filmina clasifica los FIR por la simetría de sus coeficientes. Simetría positiva es $h[n]=h[N-n-1]$: el primero igual al último, el segundo igual al anteúltimo. Simetría negativa es $h[n]=-h[N-n-1]$. En los dos casos el centro de simetría está en

$$\alpha =\frac{N-1}{2}$$

que cae sobre una muestra si $N$ es impar y entre dos si es par. Se parte de la DTFT, que para un FIR es una suma finita:

$$H(e^{j\omega })=\sum _{n=0}^{N-1}h[n]e^{-j\omega n}$$

Paso 1: sacar la fase del centro. Se saca de factor común $e^{-j\omega \alpha }$, y cada exponente queda medido desde el centro:

$$H(e^{j\omega })=e^{-j\omega \alpha }\sum _{n=0}^{N-1}h[n]e^{j\omega (\alpha -n)}$$

Paso 2: juntar cada coeficiente con su espejo. El coeficiente $n$ está a distancia $\alpha -n$ del centro, y su espejo $N-1-n$ está a la misma distancia del otro lado, $-(\alpha -n)$. Con simetría positiva los dos tienen el mismo valor, y el par suma dos exponenciales conjugadas:

$$h[n](e^{j\omega (\alpha -n)}+e^{-j\omega (\alpha -n)})=2h[n]\cos (\omega (\alpha -n))$$

Con $N$ impar sobra el coeficiente del centro, $h[\alpha ]$, que está a distancia cero y entra tal cual.

Paso 3: lo que queda es real. La suma de los pares es una suma de cosenos con coeficientes reales:

$$H(e^{j\omega })=\underset{\text{retardo puro}}{\underbrace{e^{-j\omega \alpha }}}·\underset{\text{real}}{\underbrace{A(\omega )}}$$

Toda la parte imaginaria de $H$ vive en el primer factor, que tiene módulo 1 y fase $-\alpha \omega$. $A(\omega )$ es real pero puede cambiar de signo, y donde lo hace la fase salta $\pi$: son los saltos que se ven en la banda suprimida de cualquier FIR simétrico, y no atrasan nada, porque caen en frecuencias que el filtro anula. Entre salto y salto la pendiente es siempre la misma, así que el retardo de grupo es constante:

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

Con simetría negativa el paso 2 cambia un signo: el par da $2jh[n]\mathrm{sen}(\omega (\alpha -n))$, y el coeficiente del centro, que tendría que ser igual a menos sí mismo, es cero. Lo que queda es $j$ por una suma real de senos, y esa $j$ suma $\pi /2$ a la fase:

$$\theta (\omega )=\frac{\pi }{2}-\alpha \omega$$

La pendiente es la misma, así que el retardo de grupo también es $\alpha$. Es fase lineal generalizada: la envolvente llega entera y a tiempo, y la portadora sale girada 90°. Y la conclusión no es una aproximación: es álgebra. En el ensayo del elíptico está la cuenta con los siete coeficientes de la filmina, y lo que le pasa a un filtro que no tiene esta simetría.

05 · Los cuatro tipos, y qué prohíbe cada uno

Signo de la simetría por paridad de $N$ dan los cuatro tipos de la filmina. Lo que los distingue está en las distancias al centro, $\alpha -n$: son enteras si $N$ es impar y semienteras si es par. Evaluando la suma de cosenos o de senos en $\omega =0$ y en $\omega =\pi$ aparecen ceros que ningún coeficiente puede evitar:

• Tipo 2, en $\pi$. Cada término es $\cos (\pi (k+\frac{1}{2}))=0$. Todos a la vez, así que $H(e^{j\pi })=0$.
• Tipos 3 y 4, en $0$. Cada término es un $\mathrm{sen}(0)=0$. Dicho sin cuentas: los coeficientes se cancelan de a pares, y la ganancia en continua, que es su suma, da cero.
• Tipo 3, además en $\pi$. Con distancias enteras, $\mathrm{sen}(\pi k)=0$.
• Tipo 1, en ningún lado. $\cos (0)=1$ y $\cos (\pi k)=\pm 1$: nada se anula por fuerza.

Con los ejemplos de la tabla las cuentas se hacen a mano, porque $H(e^{j0})$ es la suma de los coeficientes y $H(e^{j\pi })$ la suma con signos alternados. Dan 9 y 1 para el tipo 1, 12 y 0 para el tipo 2, 0 y 0 para el tipo 3, 0 y 4 para el tipo 4. Lo que la tabla dice es más fuerte que los ejemplos: esos ceros están para cualquier elección de coeficientes con esa simetría. Un pasaaltos de fase lineal tiene que ser tipo 1 o tipo 4. Los tipos 3 y 4, con su giro de $\pi /2$ y su cero en continua, son los naturales para diferenciadores y transformadores de Hilbert. Y el retardo $\alpha$ es entero en los tipos 1 y 3 y semientero en los 2 y 4: la salida de un FIR de $N$ par cae entre dos muestras de la entrada, lo que complica alinearla con otra señal.

Ceros forzados sobre el círculo unitario en $z=1$ o en $z=-1$: es la misma afirmación leída en el plano Z, donde la simetría obliga además a que los ceros vengan en pares recíprocos. Eso se cuenta en polos y ceros.

06 · Lo que cuesta: retardo y multiplicaciones

El retardo. Un FIR de fase lineal atrasa todas las frecuencias $(N-1)/2$ muestras, o $(N-1)/(2f_s)$ segundos. Un filtro de 41 taps atrasa 20 muestras. No deforma, porque es el mismo para todas, y se compensa descontándolo. Pero no se puede evitar: un filtro causal que promedia hacia atrás devuelve el promedio tarde, y por eso el suavizado se paga en retardo.

Las multiplicaciones. Cada muestra de salida es una suma de $N$ productos: $N$ multiplicaciones y $N-1$ sumas por muestra. Con coeficientes simétricos se pueden sumar primero las dos muestras que van a multiplicar al mismo coeficiente, y las multiplicaciones bajan a la mitad; esa estructura, y la línea de retardos que la sostiene, son del FIR transversal. Lo que no baja es el crecimiento: para una misma atenuación, el largo crece como la inversa del ancho de transición. Por eso un IIR, que con un polo cerca del círculo hace una transición angosta en pocas multiplicaciones, puede ser mucho más barato, y cuánto más depende de la transición: el «diez veces» no es constante.

07 · El diseño por ventanas, como receta

Diseñar un FIR es elegir sus $N$ coeficientes. El método que la filmina desarrolla, y el que hace fir1 por dentro, son tres pasos.

Paso 1: la respuesta ideal. Se escribe la $H_D(\omega )$ que se querría y se anti-transforma. Para un pasabajos de corte ${\omega }_c$ sale un sinc:

$$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 infinito hacia los dos lados y no causal: ningún filtro real lo implementa.

Paso 2: truncar y correr. Se toman las $N$ muestras centrales y se corren $(N-1)/2$ para que el filtro sea causal. El sinc es par, así que lo que queda es simétrico, y el corrimiento es justamente el retardo $\alpha$ de la sección 04.

Paso 3: multiplicar por una ventana. Truncar de golpe deja un rizado junto al corte que no baja con más coeficientes. Para suavizarlo se multiplica por una ventana $w[n]$ que cae a cero de a poco:

$$h[n]=h_D[n-\frac{N-1}{2}]·w[n],n=0,\dots ,N-1$$

Las ventanas de la cátedra —rectangular, Hamming, Blackman— son simétricas, así que el producto también lo es: la fase lineal sale gratis, sin pedirla. Con $N$ impar el resultado es tipo 1; con $N$ par, tipo 2, con su cero en $\pi$. La ventana se elige por la atenuación que pide la máscara de tolerancias, y $N$ por el ancho de transición: más coeficientes angostan la transición, pero no bajan el rizado, que depende sólo de la ventana. Por qué pasa eso, y cuánto atenúa cada ventana, es del fenómeno de Gibbs.

En MATLAB los tres pasos son una línea: fir1(N-1, wc/pi, ventana). Pide el orden, $N-1$, no la cantidad de coeficientes; el corte va normalizado a $\pi$; y si no se le pasa ventana usa Hamming. Hace una sola cosa más que la receta: divide por la suma de los coeficientes, para que la ganancia en continua sea exactamente 1. Con 41 coeficientes y corte en $\mathrm{0,4}\pi$, la receta sin normalizar suma 0,998433, y una vez normalizada coincide con fir1 hasta 1,1×10⁻¹⁶.

Mientras los coeficientes sean espejo unos de otros, el retardo de grupo es una recta horizontal en $(N-1)/2$, se muevan como se muevan: el módulo cambia, el retardo no. La simetría la tienen los valores, no el botón: con «Ninguna» y los coeficientes todavía espejados, nada cambia hasta que se toca uno. Tocarlo basta. En fir1(8, 0.4), llevar $h[1]$ de −0,0136 a 0,09 hace que el retardo vaya de 0,49 a 4,77 muestras según la frecuencia, medido donde $|H|$ supera el 1 % de su máximo. Y con $N$ par y simetría positiva, $|H(\pi )|$ se queda en cero por más que se arrastre: el tipo 2 no llega a $\pi$.

08 · El análisis en MATLAB

Un pasabajos de 41 coeficientes con corte en $\mathrm{0,4}\pi$, diseñado, medido y usado. Lo que figura como comentario es la salida real de correrlo en R2026a.

%% 1. El diseño: 41 coeficientes, corte en 0,4π, Hamming (la ventana por omisión)
N = 41;
h = fir1(N-1, 0.4);                  % pide el ORDEN, N-1, y el corte normalizado a π
fprintf('taps = %d   suma = %.4f   asimetria = %g\n', length(h), sum(h), max(abs(h - fliplr(h))))
% taps = 41   suma = 1.0000   asimetria = 0

%% 2. Es la receta: sinc ideal por ventana, normalizado a ganancia 1 en continua
n  = -(N-1)/2 : (N-1)/2;
hD = sin(0.4*pi*n) ./ (pi*n);
hD(n == 0) = 0.4;                    % L'Hôpital: wc/pi
hw = hD .* hamming(N).';             % hamming devuelve una columna
fprintf('max |fir1 - receta| = %.1e\n', max(abs(h - hw/sum(hw))))
% max |fir1 - receta| = 1.1e-16

%% 3. La respuesta al impulso son los coeficientes
hi = impz(h, 1);
fprintf('largo de impz = %d   max |impz - h| = %g\n', length(hi), max(abs(hi.' - h)))
figure; stem(0:N-1, h, 'filled'); grid on; xlabel('n'); ylabel('h[n]')
% largo de impz = 41   max |impz - h| = 0

%% 4. El módulo
w = linspace(0, pi, 200001);
H = freqz(h, 1, w);
HdB = 20*log10(abs(H));
figure; plot(w/pi, HdB); grid on; ylim([-90 5]); xlabel('\omega / \pi'); ylabel('|H| [dB]')
Hc = freqz(h, 1, [0.4*pi, pi]);
fprintf('|H| en 0.4pi = %.2f dB   en pi = %.2f dB\n', 20*log10(abs(Hc)))
fprintf('desvio maximo en [0, 0.3pi] = %.4f\n', max(abs(abs(H(w <= 0.3*pi)) - 1)))
fprintf('peor valor en [0.5pi, pi] = %.2f dB\n', max(HdB(w >= 0.5*pi)))
% |H| en 0.4pi = -6.01 dB   en pi = -61.16 dB
% desvio maximo en [0, 0.3pi] = 0.0033
% peor valor en [0.5pi, pi] = -53.60 dB

%% 5. El retardo de grupo: constante, (N-1)/2
disp(grpdelay(h, 1, [0.1 0.2 0.3 0.35]*pi).')
%     20    20    20    20

%% 6. filter: un coseno en la banda de paso y otro en la suprimida
nn = 0:199;
x  = cos(0.1*pi*nn) + cos(0.7*pi*nn);
y  = filter(h, 1, x);
figure; plot(nn, x, ':', nn, y); grid on; xlabel('n')
e  = y(41:end) - cos(0.1*pi*(nn(41:end) - 20));
fprintf('max |y[n] - cos(0.1pi(n-20))| para n >= 40: %.4f\n', max(abs(e)))
% max |y[n] - cos(0.1pi(n-20))| para n >= 40: 0.0011

%% 7. Los cuatro tipos: qué fuerza cada simetría en 0 y en π
T = {[1 2 3 2 1], [1 2 3 3 2 1], [1 2 0 -2 -1], [1 2 3 -3 -2 -1]};
for t = 1:4
    Ht = freqz(T{t}, 1, [0 pi]);
    g  = grpdelay(T{t}, 1, [0.2 0.6]*pi);
    fprintf('tipo %d: |H(0)| = %.3g   |H(pi)| = %.3g   retardo = %.1f\n', ...
            t, abs(Ht(1)), abs(Ht(2)), g(1))
end
% tipo 1: |H(0)| = 9   |H(pi)| = 1   retardo = 2.0
% tipo 2: |H(0)| = 12   |H(pi)| = 2.45e-16   retardo = 2.5
% tipo 3: |H(0)| = 0   |H(pi)| = 0   retardo = 2.0
% tipo 4: |H(0)| = 0   |H(pi)| = 4   retardo = 2.5

%% 8. fir1 no te deja hacer un pasaaltos de tipo 2
hp = fir1(41, 0.4, 'high');
fprintf('fir1(41, 0.4, ''high'') devuelve %d coeficientes\n', length(hp))
% Warning: Odd order symmetric FIR filters must have a gain of zero at the
% Nyquist frequency. The order is being increased by one.
% fir1(41, 0.4, 'high') devuelve 43 coeficientes

Lo que dice cada bloque:

• Bloques 1 y 3. 41 coeficientes, exactamente simétricos, y impz devuelve los mismos 41 números: la respuesta al impulso de un FIR no se calcula, se lee.
• Bloque 4. El corte de fir1 es el punto de −6 dB, mitad de amplitud, no el de −3 dB: en $\mathrm{0,4}\pi$ da −6,01 dB. Hasta $\mathrm{0,3}\pi$ el módulo se aparta de 1 en 0,0033 como mucho, y desde $\mathrm{0,5}\pi$ atenúa al menos 53,60 dB. La transición entre los dos niveles de rizado, medida entre donde $|H|$ baja de $1-\delta$ y donde baja de $\delta$, con $\delta =\mathrm{0,00209}$ el peor lóbulo de la banda suprimida, va de $\mathrm{0,319}\pi$ a $\mathrm{0,483}\pi$: mide $\mathrm{0,164}\pi$.
• Bloque 5. 20 muestras en las cuatro frecuencias, que es $(41-1)/2$.
• Bloque 6. El coseno de $\mathrm{0,7}\pi$ desaparece y el de $\mathrm{0,1}\pi$ sale entero y corrido 20 muestras. Lo que sobra, 0,0011, es casi todo el coseno rechazado que se filtra: $|H(\mathrm{0,7}\pi )|$ vale 0,00113. Las primeras 40 muestras se descartan porque ahí la línea de retardos todavía se está llenando.
• Bloque 7. La tabla de la sección 05, medida. El 2,45×10⁻¹⁶ del tipo 2 es el cero exacto con el redondeo de punto flotante encima.

Y los detalles que muerden:

• fir1 pide el orden. fir1(40, ...) da 41 coeficientes. Y si el orden es impar y se pide un pasaaltos, lo sube solo, con un aviso (bloque 8): con $N$ par el filtro sería tipo 2, y el tipo 2 no puede dejar pasar $\pi$.
• Un escalar en freqz o grpdelay no es una frecuencia. freqz(h, 1, 0.4*pi) lee el tercer argumento como la cantidad de puntos. Por eso los bloques 4, 5 y 7 le pasan siempre un vector de al menos dos frecuencias.
• hamming(N) es una columna y fir1 devuelve una fila: sin la trasposición del bloque 2, el producto elemento a elemento arma una matriz de 41 por 41 sin avisar.
• filter no descuenta el retardo. La salida está corrida $(N-1)/2$ muestras respecto de la entrada, y para compararlas hay que correr una de las dos, como hace el bloque 6.

Para tener a mano

1. $y[n]={\sum }_{k=0}^{N-1}h[k]x[n-k]$, sin realimentación. Los coeficientes son $h[n]$, y $h[n]$ dura $N$ muestras.
2. $H(z)$ es un polinomio en $z^{-1}$: $N-1$ ceros y $N-1$ polos en el origen.
3. Siempre estable: ${\sum }_k|h[k]|$ es una suma finita, para cualquier coeficiente.
4. Simetría $h[n]=\pm h[N-n-1]$ implica fase lineal: se saca $e^{-j\omega (N-1)/2}$ de factor común y lo que queda es real (o $j$ por un real). Retardo de grupo constante, $(N-1)/2$.
5. Tipo 1: nada forzado. Tipo 2: $H(\pi )=0$. Tipo 3: $H(0)=H(\pi )=0$. Tipo 4: $H(0)=0$. Un pasaaltos de fase lineal es tipo 1 o 4.
6. Cuesta $N$ multiplicaciones por muestra y $(N-1)/2$ muestras de retardo.
7. Diseño por ventanas: $h[n]=h_D[n-(N-1)/2]·w[n]$. En MATLAB, fir1(N-1, wc/pi, ventana), normalizado a ganancia 1 en continua, con el corte en −6 dB.