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.
La ecuación en diferencias general tiene una parte que mira entradas y otra que mira salidas pasadas:
Transformando los dos lados, cada retardo de muestras se vuelve un , y la función de transferencia queda como un cociente de polinomios:
Los multiplican entradas y los multiplican salidas pasadas, con . 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 son los ceros del filtro, las de sus polos.
Si todos los son cero, no hay realimentación: es un FIR, y su son los propios . Con algún distinto de cero, la salida depende de salidas anteriores: es un IIR, Infinite Impulse Response.
Un impulso entra una sola vez, en . Desde ahí la entrada es cero para siempre, pero la salida no: se calcula con , con , 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 : si se anulara a partir de algún , su transformada sería un polinomio en . Un cociente con un denominador que no se simplifica no es un polinomio. La memoria infinita sale de ahí, y cuesta sólo 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 sigue siendo FIR: su es la serie geométrica de la media móvil sumada, y el polo en se cancela con uno de los ceros del numerador. Es la idea de la estructura de muestreo en frecuencia de la filmina.
El IIR más chico que existe tiene un solo coeficiente de realimentación. La guía del TP1 lo escribe así, con :
Cada salida es la anterior, pesada por , más un poco de la entrada nueva. El factor 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, y .
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 de la guía. Corridos en MATLAB, los dos dan la misma 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.
Camino 1: meter un impulso y seguir la recurrencia. Con la salida es , y el filtro arranca en reposo, :
Desde el impulso ya no aporta nada y cada paso sólo multiplica por . El patrón queda a la vista:
Es la exponencial real, escalada. Sin el factor de ganancia, la recurrencia da pelada.
Camino 2: la serie geométrica de H(z). Paso 1, transformar la ecuación: . Paso 2, despejar el cociente:
Paso 3, reconocer la serie geométrica. Con razón y , que vale para , o sea para :
Paso 4, leer los coeficientes: por definición de transformada, el que acompaña a es . 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.
Evaluando en el círculo unitario, , y tomando módulo al cuadrado con :
Contra freqz sobre 4096 frecuencias, con , difiere en 9,4×10⁻¹⁵. Tres valores se leen directo. En continua, : para eso estaba el . En , : con es un pasabajos, con deja pasar más lo alto que lo bajo. Y el corte a −3 dB sale de igualar a , que deja . Despejando:
Tiene solución mientras el lado derecho no baje de −1, o sea para . Por debajo, ya es mayor que y el filtro no llega a cortar en ningún lado. Con da 0,10546 rad, y freqz evaluado ahí devuelve −3,0103 dB.
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:
Con la serie diverge. El filtro es estable si y sólo si , y es justo donde está el polo: se anula en . 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 , está en Polos y ceros.
Los tres casos que quedan afuera se ven distintos:
filter da , que es , y sigue creciendo un 5 % por muestra. El signo menos viene del .filter devuelve ceros.Hay una trampa. freqz no mira la estabilidad: evalúa en el círculo y devuelve lo que salga. Con entrega un pasabajos prolijo, en continua y 0,0244 en , sin avisar. Esa curva no es la respuesta en frecuencia del filtro, porque la DTFT de una que explota no converge. En el paso 3 de la sección 04 la serie pedía , y con el círculo unitario queda afuera de esa región.
Con , la respuesta al impulso es una exponencial que cae. Escrita como , su constante de tiempo es muestras. Cerca de 1, y queda la regla de bolsillo . 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,5 | 0,7227 rad | −9,5 dB | 1,44 | 2 | n = 7 |
| 0,8 | 0,2241 rad | −19,1 dB | 4,48 | 5 | n = 21 |
| 0,9 | 0,1055 rad | −25,6 dB | 9,49 | 10 | n = 44 |
| 0,95 | 0,05130 rad | −31,8 dB | 19,50 | 20 | n = 90 |
| 0,99 | 0,01005 rad | −46,0 dB | 99,50 | 100 | n = 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; contra freqz, al milésimo de dB. La última columna es el primer que cumple . Llegar a 10⁻⁶ con lleva 132 muestras, y cero exacto no llega nunca.
Es el ejercicio 7 del TP1 visto al revés: con el corte cae en 0,0513 rad, y la media móvil que corta más cerca de ahí es la de , 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.
Mientras el polo se acerca a 1, arriba la exponencial se estira y abajo el pasabajos se angosta. En 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.
La fase de es menos la fase del denominador, :
No es una recta. Derivando, el retardo de grupo es
En continua vale y en , . Con , 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 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.
En hardware los coeficientes y las muestras tienen 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 vuelve a entrar en , 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⁻⁷, 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 también se guarda en 8 bits, cada paso le saca , y cuando eso es menos de medio escalón el redondeo lo devuelve adonde estaba. La respuesta se clava en escalones, con el mayor entero que cumple , donde 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 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 y desde se queda quieta en 5·2⁻¹⁰⁷⁴. La 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:
| estructura | peor rechazo | banda de paso | polo 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 dB | 0,0127 |
| 4 biquads, 2⁻¹¹ | −79,99 dB | −0,56 a +0,01 dB | 0,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:
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 : 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.
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:
butter): sin rizado en ningún lado, la banda de paso lo más plana posible. Paga con la transición más lenta.cheby1): rizado en la banda de paso, caída monótona en la suprimida.ellip): rizado en las dos bandas. Es el de menor orden para una misma máscara.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—:
| familia | orden | biquads | peor en la banda de paso | peor rechazo |
|---|---|---|---|---|
| butter | 17 | 9 | −0,89 dB | −60,00 dB |
| cheby1 | 9 | 5 | −1,00 dB | −67,93 dB |
| ellip | 6 | 3 | −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.
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:
freqz no avisa si el filtro es inestable. El bloque 6 devuelve una curva de pasabajos normal para un filtro que explota. La estabilidad se mira aparte: max(abs(roots(a))) < 1, o isstable(b, a).impz sin largo elige uno. Con el filtro del bloque 1 cortó en 93 muestras, cuando la última todavía valía 6,2×10⁻⁶. Para comparar contra una forma cerrada hay que fijarlo: impz(b, a, 60). Y filter con un delta da lo mismo sin sorpresas.[b, a] para un IIR de orden alto. Se pide [z, p, k] y se pasa a sos. Y zp2sos con dos salidas: con una sola, mete la ganancia del elíptico de la sección 09, 1,5×10⁻³, en el numerador de la primera sección, y esos coeficientes chiquitos se redondean mal. Con los mismos 2⁻¹¹ de antes, el rechazo cae de −79,99 a −78,89 dB.freqz o grpdelay no es una frecuencia. freqz(b, a, pi) lee el escalar como cantidad de puntos. Por eso los bloques 4 a 6 le pasan siempre al menos dos frecuencias.Ejemplo_filtros.m, alpha es el peso de la entrada y el polo es 1 - alpha.freqz no lo chequea.[z, p, k] → [sos, g] = zp2sos(z, p, k) → sosfilt.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.