¿Cuál es la mejor aproximación de primer orden IIR (filtro AR) a un filtro de media móvil (filtro FIR)?

24

Suponga el siguiente filtro IIR de primer orden:

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

¿Cómo puedo elegir el parámetro st? El IIR se aproxima lo mejor posible al FIR, que es la media aritmética de las últimas muestras:αk

z[n]=1kx[n]+1kx[n1]++1kx[nk+1]

Donde , lo que significa que la entrada para el IIR podría ser más larga que k y, sin embargo, me gustaría tener la mejor aproximación de la media de las últimas k entradas.n[k,)kk

Sé que el IIR tiene una respuesta de impulso infinita, por lo tanto, estoy buscando la mejor aproximación. Estaría contento con la solución analítica, ya sea para la función de costo L2 o L1 .

¿Cómo se pueden resolver estos problemas de optimización con solo IIR de primer orden?

Gracias.

Royi
fuente
¿Tiene que seguir y[n]=αx[n]+(1α)y[n1] precisión]?
Phonon
1
Esto se convertirá en una aproximación muy pobre. ¿No puede permitirse más que un IIR de primer orden?
Leftaroundabout
Es posible que desee editar su pregunta para que no use y[n] para significar dos cosas diferentes, por ejemplo, la segunda ecuación mostrada podría leer z[n]=1kx[n]++1kx[nk+1] , y es posible que desee decir cuál es exactamente su criterio de "lo mejor posible", por ejemplo, ¿quiere |y[n]z[n]| será lo más pequeño posible para todos n , o |y[n]z[n]|2 será lo más pequeño posible para todos n .
Dilip Sarwate
@Phonon, sí, debe ser un IIR de primer orden. El criterio es simple, el resultado debe estar lo más cerca posible de la media de las últimas entradas al sistema donde . Estaría feliz de ver el resultado para ambos casos. Aunque supongo que la solución analítica solo es viable para . y[n]kn[k,inf]|y[n]z[n]|2
Royi

Respuestas:

10

No hay una solución analítica para que sea ​​un escalar (creo). Aquí hay un script que te da para una dada . Si lo necesita en línea, puede crear un LUT. El script encuentra la solución que minimizaα KααK

0πdw|H1(jw)H2(jw)|2

donde es la respuesta de frecuencia FIR y es la respuesta de frecuencia IIR.H 2H1H2

No especificó ningún rango para K. ¡Pero solo quiero dejar en claro que el siguiente sistema es equivalente a su filtro medio y tiene la misma complejidad computacional y su IIR de primer orden!

H(z)=1K1zK1z1

function a = find_a(K)

w = 0.0001:0.001:pi;
as = [-1:0.001:-0.001  0.001:0.001:1];

E = zeros(size(as));
for idx=1:length(as)
    fJ = J(w,as(idx),K);
    E(idx) = sum(fJ);
end

[Emin, indx] = min(E)
a = as(indx)

function f = J(w,a,K)
    num = 2*(2-a)*(1-cos(w*K)) + 2*(cos(w*(K-1)) - cos(w)) - 2*(1-a)*(cos(w)-cos(w*(K+1)));
    den = (2-a)^2 + 1 + (1-a)^2 + 2*(1-a)*cos(2*w) - 2*(2-a)^2*cos(w);
    f = -(a/K)*num./den;
    f = f+(1/K^2)*(1-cos(w*K))./(1-cos(w))+a^2./(1+(1-a)^2-2*(1-a)*cos(w));
end

end
niaren
fuente
@Drazick Es relativamente sencillo. Las dos expresiones para IIR y FIR están conectadas a la integral. La clave para encontrar la expresión alternativa para el filtro FIR es reconocer la serie / progresión geométrica. Encontrará todos los detalles aquí: en.wikipedia.org/wiki/Geometric_progression#Geometric_series . En el script, la función J calcula la expresión bajo el signo integral.
niaren
@niaren Sé que esta es una publicación antigua, así que si puedes recordar: ¿cómo se deriva tu función 'f'? He codificado algo similar, pero utilizando las funciones de transferencia complejas para FIR (H1) y IIR (H2) y luego haciendo sum (abs (H1 - H2) ** 2). Comparé esto con su suma (fj), pero obtuve diferentes resultados. Pensé en preguntar antes de estudiar las matemáticas.
Dom
@Dom Eso ES hace mucho tiempo y realmente no puedo recordar. Supongo que acabo de pasar por el proceso de calcular . No recuerdo cómo verifiqué la expresión. No me importa ir a través de la matemáticas de nuevo ...[H1(jω)H2(jω)][H1(jω)H2(jω)]
niaren
@niaren Hola, intenté deducir tu expresión pero me quedé atascado al sumar las fracciones complejas. Cometí un error en mi código ... su función parece dar resultados que son proporcionales a la suma (abs (H1 - H2) ** 2), lo que indica que es correcta.
Dom
16

Hay una buena discusión sobre este problema en el procesamiento de señal integrado con la arquitectura de micro señal , aproximadamente entre las páginas 63 y 69 . En la página 63 , incluye una derivación del filtro de promedio móvil recursivo exacto (que niaren dio en su respuesta ),

H(z)=1N1zN1z1.

Por conveniencia con respecto a la siguiente discusión, corresponde a la siguiente ecuación de diferencia:

yn=yn1+1N(xnxnN).

La aproximación que coloca el filtro en la forma que especificó requiere suponer que , porque (y cito de la página 68 ) " es el promedio de muestras ". Esa aproximación nos permite simplificar la ecuación de diferencia anterior de la siguiente manera: y n - 1 x nxnNyn1yn1xn

yn=yn1+1N(xnyn1)yn=yn11Nyn1+1Nxnyn=(11N)yn1+1Nxn.

Configurando , llegamos a su forma original, , que muestra que el coeficiente que desea ( con respecto a esta aproximación) es exactamente (donde es el número de muestras). yn=αxn+(1-α)yn-11α=1Nyn=αxn+(1α)yn1 N1NN

¿Es esta aproximación la "mejor" en algún aspecto? Ciertamente es elegante. Así es como se compara la respuesta de magnitud [a 44,1 kHz] para N = 3, y cuando N aumenta a 10 (aproximación en azul):

N = 3 N = [3,10]


Como sugiere la respuesta de Peter , aproximar un filtro FIR con un filtro recursivo puede ser problemático bajo una norma de mínimos cuadrados. Se puede encontrar una extensa discusión sobre cómo resolver este problema en general en la tesis de JOS, Técnicas para el diseño de filtros digitales e identificación del sistema con aplicación al violín . Él aboga por el uso de la Norma Hankel, pero en los casos en que la respuesta de fase no importa, también cubre el Método de Kopec, que podría funcionar bien en este caso (y utiliza una norma ). Aquí se puede encontrar una visión general de las técnicas de la tesis . Pueden producir otras aproximaciones interesantes.L2

Datageist
fuente
Esta es una forma "elegante" de decir algo sobre la memoria del filtro IIR de primer orden. Su memoria es equivalente a . Examinaré los otros métodos que mencionaste. Gracias. 1α
Royi
¿Podría explicar intuitivamente por qué bajo la norma LS ( ) no hay solución? L2
Royi
No estoy seguro de si hay una solución LS en este caso o no todavía, solo sé que tiende a tener problemas de convergencia para el problema general de "aproximación FIR basada en IIR". Actualizaré con más información cuando tenga la oportunidad.
Datageist
Bueno, si la función de costo que Peter sugirió (la primera) es correcta, hay una solución. Al menos según mis cálculos.
Royi
Excelente. Tengo curiosidad por ver cómo el enfoque "heurístico" se compara con algo más canónico. 1/N
Datageist
16

OK, tratemos de obtener lo mejor: tan que el coeficiente de es .

y[n]=αx[n]+(1α)y[n1]=αx[n]+(1α)αx[n1]+(1α)2y[n2]=αx[n]+(1α)αx[n1]+(1α)2αx[n2]+(1α)3y[n3]
x[nm]α(1α)m

La mejor aproximación al cuadrado medio minimizará:

J(α)=m=0k1(α(1α)m1k)2+m=kα2(1α)2m=m=0k1(α2(1α)2m2kα(1α)m+1k2)+α2(1α)2km=0(1α)2m=α21(1α)2k1(1α)2+2αk1(1α)k1(1α)+α2(1α)2k1(1α)2+1k=α21(1α)2+2k(1(1α)k)+1k=α22αα2+2k(1(1α)k)+1k=α2α+2k(1(1α)k)+1k
porque los coeficientes FIR son cero para .m>k1

El siguiente paso es tomar derivados y equiparar a cero.


Al observar una gráfica de la derivada para y de 0 a 1, parece que el problema (como lo configuré) está mal planteado, porque la mejor respuesta es .JK=1000αα=0

ingrese la descripción de la imagen aquí


Creo que hay un error aquí. La forma en que debería ser según mis cálculos es:

J(α)=m=0k1(α(1α)m1k)2+m=kα2(1α)2m=m=0k1(α2(1α)2m2kα(1α)m+1k2)+α2(1α)2km=0(1α)2m=α21(1α)2k1(1α)22αk1(1α)k1(1α)+1k+α2(1α)2k1(1α)2

Simplificándolo según los rendimientos de Mathematica :

J(α)=α2α+2(1α)k1k

El uso del siguiente código en MATLAB produce algo equivalente aunque diferente:

syms a k;

expr1 = (a ^ 2) * ((1 - ((1 - a) ^ (2 * k))) / (1 - ((1 - a) ^ 2)));
expr2 = ((2 * a) / k) * ((1 - ((1 - a) ^ (k))) / (1 - (1 - a)));
expr3 = (1 / k);
expr4 = ((a ^ 2) * ((1 - a) ^ (2 * k))) / (1 - ((1 - a) ^ (2)));

simpExpr = simplify(expr1 - expr2 + expr3 + expr4);

J(α)=2α2k2(1α)k+1k

De todos modos, esas funciones tienen un mínimo.


Supongamos que realmente solo nos importa la aproximación sobre el soporte (longitud) del filtro FIR. En ese caso, el problema de optimización es solo:

J2(α)=m=0k1(α(1α)m1k)2

El trazado de para varios valores de versus da como resultado la fecha en los gráficos y la tabla a continuación.J2(α)Kα

Para = 8. = 0.1533333 Para = 16. = 0.08 Para = 24. = 0.0533333 Para = 32. = 0.04 para = 40. = 0.0333333 para = 48. = 0.0266667 para = 56. = 0.0233333 para = 64. Kαmin
Kαmin
Kαmin
Kαmin
Kαmin
Kαmin
Kαmin
Kαmin = 0.02
Para = 72. = 0.0166667 Kαmin

ingrese la descripción de la imagen aquí

Las líneas discontinuas rojas son y las líneas verdes son , el valor de que minimiza (elegido de ).1/KαminαJ2(α)alpha=[0:.01:1]/3;

Peter K.
fuente
1
Solo iba a publicar exactamente lo mismo =)
Phonon
@Phonon: ¡Siéntete libre de continuar! Lo marqué como wiki de la comunidad para ese propósito.
Peter K.
La derivada wrt es una serie con un número infinito de términos (es decir, no un polinomio) que debe establecer igual a y luego resolver para , por lo que será necesario tener cuidado (o posiblemente una aproximación). α0α
Dilip Sarwate
¿Alguien puede verificar y / o corregir mi trabajo? :-)
Peter K.
@DilipSarwate, ¿Cuál sería la mejor aproximación? Gracias.
Royi
3

Basado en pruebas experimentales con un krango (2 a 100), el mejor ajuste (suma de error al cuadrado) da una relación de alfa = 1/k^0.865 ser el knúmero de muestras para el filtro MovAvg

chiliwili69
fuente
Utilicé una hoja de Excel interactiva para extraer la relación: drive.google.com/open?id=0B48gEYiKwYegY3NqQThMTTZ1OG8
chiliwili69
3

Me topé con esta vieja pregunta y me gustaría compartir mi solución. Como se menciona en otras respuestas, no existe una solución analítica, pero la función que se debe minimizar se comporta muy bien y el valor óptimo de se puede encontrar fácilmente con unas pocas iteraciones de Newton. También hay una fórmula para verificar la optimización del resultado.α

La respuesta al impulso del filtro de media móvil FIR de longitud viene dada porN

(1)hFIR[n]=1N(u[n]u[nN])

donde es la función de paso unitario. El filtro de primer orden IIRu[n]

(2)y[n]=αx[n]+(1α)y[n1]

tiene la respuesta impulsiva

(3)hIIR[n]=α(1α)nu[n]

El objetivo ahora es minimizar el error al cuadrado

(4)ϵ=n=0(hFIR[n]hIIR[n])2

Usando y , el error se puede escribir como(1)(3)

ϵ(α)=n=0N1(α(1α)n1N)2+n=Nα2(1α)2n=α2n=0(1α)2n2αNn=0N1(1α)n+n=0N11N2=α21(1α)22αN1(1α)N1(1α)+1N(5)=α2α2N(1(1α)N)+1N,0<α<2

Esta expresión es muy similar a la dada en esta respuesta , pero no es idéntica. La restricción en en asegura que la suma infinita converja, y es idéntica a la condición de estabilidad para el filtro IIR dada por .α(5)(2)

Establecer la derivada de a cero da como resultado(5)

(6)(1α)N1(2α)2=1

Tenga en cuenta que el óptimo debe estar en el intervalo porque los valores mayores de dan como resultado una respuesta de impulso alterna , que no puede aproximarse a la respuesta de impulso constante del filtro de promedio móvil FIR.α(0,1]α(3)

Tomando la raíz cuadrada de e introduciendo , obtenemos(6)β=1α

(7)β(N+1)/2+β(N1)/21=0

Esta ecuación no se puede resolver analíticamente para , pero se puede resolver para :NβN

(8)N=2log(1+β)log(β),β0

La ecuación se puede usar para verificar una solución numérica de ; debe devolver el valor especificado de .( 7 ) N(8)(7)N

La ecuación se puede resolver con unas pocas líneas de código (Matlab / Octave):(7)

N = 50; % de longitud de filtro deseada del filtro de promedio móvil FIR

if (N == 1)% sin iteración para caso trivial
    b = 0;
más
    % De iteración de Newton
    b = 1; % valor inicial
    Nit = 7;
    n = (N + 1) / 2;
    para k = 1: Nit,
        f = b ^ n + b ^ (n-1) -1;
        fp = n * b ^ (n-1) + (n-1) * b ^ (n-2);
        b = b - f / fp;
    fin

    % comprobar resultado
    N0 = -2 * log (1 + b) / log (b) + 1% debe ser igual a N
fin

a = 1 - b;

A continuación se muestra una tabla con los valores óptimos de para un rango de longitudes de filtro :NαN

   N alfa

   1 1.0000e + 00
   2 5.3443e-01
   3 3.8197e-01
   4 2.9839e-01
   5 2.4512e-01
   6 2.0809e-01
   7 1.8083e-01
   8 1.5990e-01
   9 1.4333e-01
  10 1.2987e-01
  20 6.7023e-02
  30 4.5175e-02
  40 3.4071e-02
  50 2.7349e-02
  60 2.2842e-02
  70 1.9611e-02
  80 1.7180e-02
  90 1.5286e-02
 100 1.3768e-02
 200 6.9076e-03 
 300 4.6103e-03
 400 3.4597e-03
 500 2.7688e-03
 600 2.3078e-03
 700 1.9785e-03
 800 1.7314e-03
 900 1.5391e-03
1000 1.3853e-03
Matt L.
fuente