Bitácora · DCSS · Concepto

El filtro IIR

Unidad 1 · Filtrado digital · realimentación, el filtro de un polo, estabilidad, biquads y análisis en MATLAB

Un filtro IIR calcula cada salida con entradas y también con salidas que ya calculó. Esa realimentación le da una respuesta al impulso que no termina nunca, hecha con un puñado de coeficientes. Tiene tres costos: el filtro puede ser inestable, la fase no es lineal y los errores de redondeo vuelven a entrar al cálculo.

01 · La definición

La ecuación en diferencias general tiene una parte que mira entradas y otra que mira salidas pasadas:

y[n]=∑k=0Nbkx[n−k]−∑k=1Maky[n−k]

Transformando los dos lados, cada retardo de k muestras se vuelve un z−k, y la función de transferencia queda como un cociente de polinomios:

H(z)=B(z)A(z)=∑k=0Nbkz−k1+∑k=1Makz−k

Los bk multiplican entradas y los ak multiplican salidas pasadas, con a0=1. Es el orden en que los recibe filter(b, a, x) y el que usan los scripts de clase; la filmina escribe la misma ecuación con las dos letras intercambiadas. Las raíces de B(z) son los ceros del filtro, las de A(z) sus polos.

Si todos los ak son cero, no hay realimentación: es un FIR, y su h[n] son los propios bk. Con algún ak distinto de cero, la salida depende de salidas anteriores: es un IIR, Infinite Impulse Response.

02 · Por qué la respuesta al impulso no se apaga

Un impulso entra una sola vez, en n=0. Desde ahí la entrada es cero para siempre, pero la salida no: y[1] se calcula con y[0], y[2] con y[1], y así. Cada salida le pasa una parte de sí misma a la siguiente. Puede achicarse, pero no tiene de dónde sacar un cero exacto.

Dicho en z: si h[n] se anulara a partir de algún n, su transformada sería un polinomio en z−1. Un cociente B/A con un denominador que no se simplifica no es un polinomio. La memoria infinita sale de ahí, y cuesta sólo M multiplicaciones.

La condición es que el denominador no se simplifique: recursivo no es lo mismo que IIR. La media móvil se puede calcular recursivamente, y[n]=y[n−1]+(x[n]−x[n−M])/M, y sigue siendo FIR: su H(z)=(1−z−M)/(M(1−z−1)) es la serie geométrica de la media móvil sumada, y el polo en z=1 se cancela con uno de los ceros del numerador. Es la idea de la estructura de muestreo en frecuencia de la filmina.

03 · El IIR de un polo de la cátedra

El IIR más chico que existe tiene un solo coeficiente de realimentación. La guía del TP1 lo escribe así, con α<1:

y[n]=αy[n−1]+(1−α)x[n]

Cada salida es la anterior, pesada por α, más un poco de la entrada nueva. El factor (1−α) sólo escala: está para que la ganancia en continua sea 1. Toda la dinámica vive en α. En la notación de la sección anterior, b0=1−α y a1=−α.

Ejemplo_filtros.m arma el mismo filtro con b = alpha; a = [1, -(1-alpha)]. Ahí la letra nombra al otro coeficiente: alpha pesa la entrada y el polo queda en 1 - alpha. El alpha = 0.1 del script es el α=0,9 de la guía. Corridos en MATLAB, los dos dan la misma h[n] a 2,8×10⁻¹⁷, que es el redondeo de la máquina. En esta nota α es siempre el de la guía: la posición del polo.

04 · La respuesta al impulso, por dos caminos

Camino 1: meter un impulso y seguir la recurrencia. Con x[n]=δ[n] la salida es h[n], y el filtro arranca en reposo, h[−1]=0:

h[0]=α·0+(1−α)·1=1−αh[1]=α(1−α)+(1−α)·0=α(1−α)h[2]=α·α(1−α)+0=α2(1−α)

Desde n=1 el impulso ya no aporta nada y cada paso sólo multiplica por α. El patrón queda a la vista:

h[n]=(1−α)αnu[n]

Es la exponencial real, escalada. Sin el factor de ganancia, la recurrencia y[n]=x[n]+αy[n−1] da h[n]=αnu[n] pelada.

Camino 2: la serie geométrica de H(z). Paso 1, transformar la ecuación: Y(z)=αz−1Y(z)+(1−α)X(z). Paso 2, despejar el cociente:

H(z)=Y(z)X(z)=1−α1−αz−1

Paso 3, reconocer la serie geométrica. Con razón r=αz−1 y 1/(1−r)=∑n=0∞rn, que vale para |r|<1, o sea para |z|>|α|:

H(z)=∑n=0∞(1−α)αn⏟h[n]z−n

Paso 4, leer los coeficientes: por definición de transformada, el que acompaña a z−n es h[n]. Es el mismo resultado por otra cuenta: uno itera la ecuación, el otro desarrolla una fracción. Contra filter con un delta a la entrada, la forma cerrada difiere en 1,4×10⁻¹⁷ sobre 60 muestras.

05 · La respuesta en frecuencia

Evaluando en el círculo unitario, z=ejω, y tomando módulo al cuadrado con |1−αe−jω|2=(1−αcosω)2+α2sen2ω:

H(ejω)=1−α1−αe−jω,|H|2=(1−α)21−2αcosω+α2

Contra freqz sobre 4096 frecuencias, con α=0,9, difiere en 9,4×10⁻¹⁵. Tres valores se leen directo. En continua, H(ej0)=1: para eso estaba el (1−α). En ω=π, |H|=(1−α)/(1+α): con α>0 es un pasabajos, con α<0 deja pasar más lo alto que lo bajo. Y el corte a −3 dB sale de igualar |H|2 a 1/2, que deja 1−2αcosω+α2=2(1−α)2. Despejando:

cosω3dB=4α−1−α22α

Tiene solución mientras el lado derecho no baje de −1, o sea para α≥3−22≈0,1716. Por debajo, |H(ejπ)| ya es mayor que 1/2 y el filtro no llega a cortar en ningún lado. Con α=0,9 da 0,10546 rad, y freqz evaluado ahí devuelve −3,0103 dB.

06 · Estabilidad: el polo tiene que quedar adentro

Un sistema LTI es estable BIBO —entrada acotada, salida acotada— si y sólo si su respuesta al impulso es absolutamente sumable. Para el filtro de un polo, la suma es otra serie geométrica:

∑n=0∞|h[n]|=|1−α|∑n=0∞|α|n=|1−α|1−|α|si |α|<1

Con |α|≥1 la serie diverge. El filtro es estable si y sólo si |α|<1, y α es justo donde está el polo: A(z)=1−αz−1 se anula en z=α. Es la regla general enunciada para un polo: un IIR causal es estable si todos sus polos quedan adentro del círculo unitario. Por qué vale para cualquier cantidad de polos, y cómo se lee en el plano z, está en Polos y ceros.

Los tres casos que quedan afuera se ven distintos:

Hay una trampa. freqz no mira la estabilidad: evalúa B/A en el círculo y devuelve lo que salga. Con α=1,05 entrega un pasabajos prolijo, |H|=1 en continua y 0,0244 en π, sin avisar. Esa curva no es la respuesta en frecuencia del filtro, porque la DTFT de una h[n] que explota no converge. En el paso 3 de la sección 04 la serie pedía |z|>|α|, y con |α|>1 el círculo unitario queda afuera de esa región.

07 · Cómo lo cambia el polo

Con 0<α<1, la respuesta al impulso es una exponencial que cae. Escrita como αn=e−n/nτ, su constante de tiempo es nτ=−1/lnα muestras. Cerca de 1, lnα≈−(1−α) y queda la regla de bolsillo nτ≈1/(1−α). La memoria y el ancho de banda son la misma perilla: acercar el polo a 1 alarga la memoria y angosta el pasabajos.

αcorte −3 dB|H(π)|nτ = −1/ln α1/(1−α)αⁿ < 1 % desde
0,50,7227 rad−9,5 dB1,442n = 7
0,80,2241 rad−19,1 dB4,485n = 21
0,90,1055 rad−25,6 dB9,4910n = 44
0,950,05130 rad−31,8 dB19,5020n = 90
0,990,01005 rad−46,0 dB99,50100n = 459

Todas las columnas son forma cerrada. El corte, cruzado en MATLAB con fzero sobre el módulo, coincide en las cinco filas al quinto decimal; |H(π)| contra freqz, al milésimo de dB. La última columna es el primer n que cumple n>ln0,01/lnα. Llegar a 10⁻⁶ con α=0,9 lleva 132 muestras, y cero exacto no llega nunca.

Es el ejercicio 7 del TP1 visto al revés: con α=0,95 el corte cae en 0,0513 rad, y la media móvil que corta más cerca de ahí es la de M=54, que corta en 0,0516 rad (con 55, en 0,0506). Son dos multiplicaciones por muestra contra cincuenta y cuatro, calculando la media como FIR directo. Cuánto vale de verdad esa ventaja cuando la máscara se pone exigente está en «Diez veces más barato» no es una constante.

Simulación Un polo que se corre hasta el círculo — arriba, h[n]/(1−α) = αⁿ, lo que devuelve la recurrencia; abajo, |H(e^jω)| de la forma de la guía entre 0 y π. Arranca en α = 0,5, bien adentro: llevalo hasta 1 y pasate.
Σ|h[n]| — constante de tiempo — corte a −3 dB — retardo de grupo en ω = 0 — h[60]/h[0] — con 8 bits —

Mientras el polo se acerca a 1, arriba la exponencial se estira y abajo el pasabajos se angosta. En α=1 la exponencial se vuelve un escalón, y un paso más allá empieza a crecer. Abajo la curva sigue ahí, y hasta tiene la misma forma: es lo que devuelve freqz, que no sabe que el filtro explotó. Con 8 bits aparece otra cosa, más escondida: con la vista en dB, la respuesta cae y de golpe se clava en un piso, a unos pocos escalones del cero. Con α negativo se clava oscilando. Es el efecto de palabra finita de la sección 09, en su versión más chica.

08 · La fase no es lineal

La fase de H(ejω) es menos la fase del denominador, 1−αcosω+jαsenω:

θ(ω)=−arctanαsenω1−αcosω

No es una recta. Derivando, el retardo de grupo es

τ(ω)=−dθdω=αcosω−α21−2αcosω+α2

En continua vale α/(1−α) y en π, −α/(1+α). Con α=0,9, grpdelay da 9, −0,3232, −0,4475 y −0,4737 muestras en 0, π/4, π/2 y π, los mismos cuatro números que la fórmula. Lo bajo sale atrasado nueve muestras y lo alto sale un poco adelantado: una señal con las dos cosas llega deformada.

No es un defecto de este filtro en particular. Fase lineal exacta pide una h[n] simétrica, y una respuesta causal e infinita no tiene respecto de qué ser simétrica. Ningún IIR causal tiene fase lineal. En un elíptico de orden alto el retardo no se queda en nueve muestras, se dispara justo en el borde de la banda: eso está en La fase no lineal mide 60 muestras justo donde te importa.

09 · Palabra finita, y por qué se implementa en biquads

En hardware los coeficientes y las muestras tienen B bits, no los 53 de mantisa de un double. El IIR es el que más lo sufre, por las dos razones de siempre: el error que se comete en y[n−1] vuelve a entrar en y[n], y los polos se mueven cuando se redondean los coeficientes.

En el filtro de un polo, el polo es el coeficiente. Redondear α lo corre a lo sumo medio escalón, y eso importa cuando está cerca del borde. Con 8 bits, un escalón de 2⁻⁷, α=0,999 se redondea a 1 exacto: el polo queda sobre el círculo y el filtro deja de ser estable.

Redondear el estado frena la caída. Si y[n] también se guarda en 8 bits, cada paso le saca (1−|α|)|y|, y cuando eso es menos de medio escalón el redondeo lo devuelve adonde estaba. La respuesta se clava en k escalones, con k el mayor entero que cumple k(1−|αq|)≤1/2, donde αq es el coeficiente ya redondeado. Es la banda muerta. Simulado en MATLAB, con la recurrencia sin ganancia y un impulso de altura 1, se clava en 1, 2, 4 y 10 escalones para α = 0,5, 0,8, 0,9 y 0,95; con α=−0,9 queda saltando entre +4 y −4 escalones para siempre, que es un ciclo límite. En double pasa lo mismo, sólo que mucho más abajo: filter(0.1, [1 -0.9], x), con un impulso como x, entra en los subnormales en n=6702 y desde n=7028 se queda quieta en 5·2⁻¹⁰⁷⁴. La h[n] que en papel no se apaga tampoco se apaga en la máquina, pero por otra razón.

Con orden alto, la forma directa se rompe. Un denominador de grado 8 con los polos amontonados cerca del círculo tiene raíces hipersensibles a sus coeficientes. El elíptico ellip(8, 0.5, 80, 0.3) tiene en a coeficientes de hasta 21,98, que se cancelan entre sí para dejar polos de radio 0,978. Redondeando todos los coeficientes a múltiplos de 2⁻¹¹, medido sobre 20 001 frecuencias:

estructurapeor rechazobanda de pasopolo más corrido
diseño, en double−80,00 dB−0,50 a 0 dB—
forma directa, 2⁻¹¹−48,44 dB−0,81 a +1,42 dB0,0127
4 biquads, 2⁻¹¹−79,99 dB−0,56 a +0,01 dB0,00069

El rechazo se mide desde donde el diseño llega a −80 dB, en 0,3748π; la forma directa falla justo ahí, y desde 0,4π todavía rechaza 56,20 dB. Con la misma cuantización y el mismo filtro, sólo cambió cómo se escribe. En la cascada cada polo depende de los dos coeficientes de su sección, no de los ocho del polinomio entero, y se corre unas dieciocho veces menos. La cascada tampoco sale perfecta: la banda de paso se pasa 0,06 dB del rizado pedido. Por eso el cuarto paso del diseño, en la filmina, es analizar los efectos de palabra finita.

Cada sección es un biquad, un cociente de dos cuadráticas, y el filtro entero es su producto:

H(z)=g∏k=1Lb0k+b1kz−1+b2kz−21+a1kz−1+a2kz−2

Se agrupan de a dos porque un par de polos complejos conjugados da una sección de coeficientes reales. En MATLAB cada sección es una fila [b0 b1 b2 1 a1 a2] de la matriz sos, second-order sections.

Y la forma directa ni siquiera necesita cuantización para romperse. [b, a] = butter(16, 0.05), en double, devuelve un a cuyo polo más lejano tiene radio 1,0668, afuera del círculo, cuando el Butterworth tiene todos los polos adentro: pedido como ceros, polos y ganancia, el más lejano está en 0,9848. Con orden 15 ya da 1,0267. Depende del corte, porque con el corte bajo los polos se amontonan cerca de z=1: el Butterworth de orden 17 de la sección siguiente, con corte en 0,208π, sale bien en [b, a], a 8×10⁻⁹ de la cascada.

10 · Los diseños clásicos, como receta

Un IIR selectivo no se diseña a mano. Se parte de un filtro analógico clásico, que la teoría ya dejó resuelto, y se lo digitaliza con la transformación bilineal, que es lo que hacen por dentro butter, cheby1 y ellip. Las familias se distinguen por dónde aceptan rizado:

La receta tiene cuatro pasos: pedir el orden mínimo con buttord, cheb1ord o ellipord; diseñar pidiendo ceros, polos y ganancia; pasar a secciones con zp2sos; y verificar contra la máscara con freqz. Para una misma máscara —banda de paso hasta 0,2π con 1 dB, suprimida desde 0,3π con 60 dB—:

familiaordenbiquadspeor en la banda de pasopeor rechazo
butter179−0,89 dB−60,00 dB
cheby195−1,00 dB−67,93 dB
ellip63−1,00 dB−60,00 dB

Las tres cumplen la máscara, medido con freqz sobre las secciones. Un orden impar deja una sección de primer orden, con el tercer coeficiente en cero: por eso 17 da 9 biquads. El elíptico hace con 3 secciones lo que el Butterworth hace con 9, a cambio de la fase más torcida cerca del corte.

11 · El análisis en MATLAB

El análisis completo del filtro de un polo, más la receta y la trampa de la forma directa. Lo que figura como comentario es la salida real de correrlo en R2026a. El bloque 1 es el IIR de Ejemplo_filtros.m, tal cual está escrito ahí.

%% 1. El IIR de un polo, escrito como en Ejemplo_filtros.m
alpha = 0.1;
b = alpha;
a = [1, -(1-alpha)];                 % el polo queda en 1 - alpha = 0.9
p = 1 - alpha;

%% 2. La respuesta al impulso, obtenida: se le mete un delta
d = [1, zeros(1, 59)];
h = filter(b, a, d);
disp(h(1:5))
%     0.1000    0.0900    0.0810    0.0729    0.0656
n = 0:59;
fprintf('max |filter - cerrada| = %.1e\n', max(abs(h - alpha*p.^n)))
% max |filter - cerrada| = 1.4e-17

%% 3. impz elige solo cuántas muestras devolver
[hi, ni] = impz(b, a);
fprintf('impz: %d muestras, la última vale %.2e\n', numel(hi), hi(end))
% impz: 93 muestras, la última vale 6.17e-06

%% 4. La respuesta en frecuencia, contra la forma cerrada
[H, w] = freqz(b, a, 4096, 'whole');
Hc = alpha ./ (1 - p*exp(-1j*w));
fprintf('max |freqz - cerrada| = %.1e\n', max(abs(H - Hc)))
% max |freqz - cerrada| = 9.4e-15
w3 = acos((-1 + 4*p - p^2) / (2*p));
fprintf('w3 = %.5f rad, |H| ahí = %.4f dB\n', w3, 20*log10(abs(freqz(b, a, [w3 pi])) * [1; 0]))
% w3 = 0.10546 rad, |H| ahí = -3.0103 dB

%% 5. El retardo de grupo no es constante
wg = [0 pi/4 pi/2 pi];
disp(grpdelay(b, a, wg).')
%     9.0000   -0.3232   -0.4475   -0.4737
disp((p*cos(wg) - p^2) ./ (1 - 2*p*cos(wg) + p^2))
%     9.0000   -0.3232   -0.4475   -0.4737

%% 6. Un polo afuera del círculo: filter explota y freqz no avisa
p2 = 1.05;
h2 = filter(1 - p2, [1, -p2], d);
fprintf('h[59] = %.3f\n', h2(end))
% h[59] = -0.889
disp(abs(freqz(1 - p2, [1, -p2], [0 pi])))
%     1.0000    0.0244

%% 7. La receta: orden mínimo, ceros-polos-ganancia, secciones de orden 2
Wp = 0.2;  Ws = 0.3;  Rp = 1;  Rs = 60;
[N, Wn] = ellipord(Wp, Ws, Rp, Rs);
[z, pz, k] = ellip(N, Rp, Rs, Wn);
[sos, g] = zp2sos(z, pz, k);
[Hs, ws] = freqz(sos, 8192);
Hs = g * Hs;
fprintf('orden %d, %d secciones\n', N, size(sos, 1))
% orden 6, 3 secciones
fprintf('paso: %.3f dB   rechazo: %.2f dB\n', ...
    min(20*log10(abs(Hs(ws <= Wp*pi)))), max(20*log10(abs(Hs(ws >= Ws*pi)))))
% paso: -1.000 dB   rechazo: -60.00 dB
y = g * sosfilt(sos, randn(1000, 1));  % así se filtra con la cascada

%% 8. La forma directa se rompe sola, sin cuantizar nada
[bb, ab] = butter(16, 0.05);
[zb, pb, kb] = butter(16, 0.05);
fprintf('polo más lejano: [b,a] %.4f   zpk %.4f\n', max(abs(roots(ab))), max(abs(pb)))
% polo más lejano: [b,a] 1.0668   zpk 0.9848

Los detalles que muerden:

Para tener a mano

  1. IIR: hay realimentación, algún ak≠0, y H=B/A no se simplifica a un polinomio. Por eso h[n] no se apaga. Recursivo no implica IIR.
  2. El de un polo de la guía: y[n]=αy[n−1]+(1−α)x[n], con h[n]=(1−α)αnu[n] y H(z)=(1−α)/(1−αz−1). En Ejemplo_filtros.m, alpha es el peso de la entrada y el polo es 1 - alpha.
  3. Estable si y sólo si |α|<1: el polo adentro del círculo. freqz no lo chequea.
  4. Corte: cosω3dB=(4α−1−α2)/(2α). Constante de tiempo: −1/lnα≈1/(1−α) muestras.
  5. Fase no lineal: el retardo de grupo va de α/(1−α) en continua a −α/(1+α) en π. Ningún IIR causal tiene fase lineal.
  6. Palabra finita: el coeficiente redondeado mueve el polo, el estado redondeado deja una banda muerta o un ciclo límite, y en orden alto la forma directa pierde la máscara. Se implementa en biquads: [z, p, k] → [sos, g] = zp2sos(z, p, k) → sosfilt.
  7. Receta: buttord/cheb1ord/ellipord, diseño en ceros-polos-ganancia, zp2sos, verificar con freqz. Para una misma máscara el elíptico da el menor orden.
← Anterior Siguiente → Inicio
Relacionados El filtro FIR Polos y ceros La media móvil La fase del elíptico