¿Es numéricamente más estable implementar el filtrado como multiplicación o convolución?

12

Estoy escribiendo un programa para filtrar una señal de 20,000 muestras con un filtro Butterworth de quinto orden fuera de línea. Tengo acceso a una implementación de FFT. Parece que hay dos alternativas para implementar el filtrado:

  • convolucionando la señal con la respuesta al impulso en el dominio del tiempo, o
  • multiplicando la señal con la respuesta al impulso en el dominio de frecuencia y transformando inversamente el resultado

Estos métodos serían idénticos en el caso teórico de FT. Sin embargo, hacerlo en la vida real con el DFT, supongo que las cosas son diferentes. ¿Es uno de los métodos numéricamente más estable? ¿Hay otros problemas que debería tener en cuenta? El número de cálculos no es importante.

filters Andreas
fuente
El método FFT será mucho más rápido de calcular para la mayoría de las longitudes de señal. Solo las longitudes cortas son más rápidas con la convolución en el dominio del tiempo.
endolito el

Respuestas:

5

Con convolución, no se encontrará con ningún problema de estabilidad, ya que no existe un filtrado recursivo, por lo que no acumulará ningún error. En otras palabras, el sistema es todo ceros, no polos. He escuchado anecdóticamente, pero no verificado por mí mismo, que la convolución basada en FFT tiene un error menor que la convolución en el dominio del tiempo, simplemente porque tiene operaciones aritméticas O (n log n) en lugar de O (n ^ 2).

Por lo general, hasta donde yo sé, los filtros de Butterworth se implementan como filtros recursivos (IIR), por lo que ese es un tema diferente. Los filtros IIR tienen polos y ceros, por lo que puede haber problemas de estabilidad en la práctica. Además, para los filtros IIR, los métodos basados ​​en FFT no son una opción, pero por el lado positivo, los filtros IIR tienden a ser de muy bajo orden.

En cuanto a los problemas de estabilidad con los filtros IIR, tienden a tener problemas en los pedidos más altos: simplemente arrojaré un número y diré que aproximadamente el sexto orden lo está presionando. En cambio, generalmente se implementan como biquads en cascada (secciones de filtro de segundo orden). Para su filtro de quinto orden, escríbalo como una función de transferencia de dominio z (será una función racional de quinto grado), y luego factoríelo en sus 5 polos y 5 ceros. Recoge los conjugados complejos y tendrás dos biquads y un filtro de primer orden. En general, los problemas de estabilidad tienden a surgir a medida que los polos se acercan al círculo unitario.

También puede haber problemas con el ruido y los ciclos límite en los filtros IIR, por lo que hay diferentes topologías de filtro (es decir, forma directa I, forma directa II) que tienen diferentes propiedades numéricas, pero no pensaría demasiado en este punto, solo use doble- precisión y seguramente será lo suficientemente bueno.

schnarf
fuente
¿Qué quieres decir con que solo funciona para los filtros FIR? Supuse que los filtros IIR tendrían que ser muestreados de alguna manera de todos modos. ¿Los filtros IIR generalmente se implementan en el dominio del tiempo para evitar esto?
Andreas
1
Que yo sepa, los filtros IIR siempre se implementan en el dominio del tiempo. Un filtro IIR (aquí, por ejemplo, un filtro de segundo orden o "biquad") se define mediante una ecuación de diferencia como y(n) = b0 * x(n) + b1 * x(n-1) + b2 * x(n-2) - a1 * y(n-1) - a2 * y(n-2). Tenga en cuenta que esta es una combinación de las muestras de entrada anteriores (los valores de x) y las muestras de salida anteriores (los valores de y). Un filtro FIR solo depende de las entradas pasadas, por lo que admite una implementación eficiente del dominio de frecuencia. Un filtro IIR no lo hace, pero de todos modos es muy eficiente porque los filtros IIR tienden a ser mucho más bajos.
schnarf
1
La razón por la que los filtros IIR tienden a ser mucho más bajos es que los polos (la retroalimentación de las muestras de salida anteriores) permiten que el filtro se vuelva mucho más empinado con muy pocos coeficientes en comparación con un filtro FIR. Cuando digo un orden mucho más bajo, un filtro IIR típico podría ser de segundo orden (5 coeficientes), mientras que un filtro FIR típico para la misma tarea podría tener miles de coeficientes.
schnarf
4

En casi todos los casos, su mejor opción no es convolución ni FFT, sino simplemente aplicar el filtro IIR directamente (usando, por ejemplo, la función sosfilt ()). Esto será mucho más eficiente en términos de consumo de CPU y memoria.

Si hace una diferencia numérica depende del filtro específico. El único caso en el que puede surgir alguna diferencia es si los polos están muy, muy cerca del círculo unitario. Incluso hay algunos trucos que pueden ayudar. NO USE la representación de la función de transferencia y filter (), pero use polos y ceros con sosfilt (). Aquí hay un ejemplo de la diferencia.

n = 2^16;  % filter length
fs = 44100; % sample rate
x = zeros(n,1); x(1) = 1;
f0 = 15; % cutoff frequency in Hz
% design with poles and zeroes
[z,p,k] = butter(5,f0*2/fs);
clf
plot(sosfilt(zp2sos(z,p,k),x));
% design with transfer function
[b,a] = butter(5,f0*2/fs);
hold on
plot(filter(b,a,x),'k');

filter () falla en un corte de aproximadamente 15Hz a 44.1kHz. Para sosfilt (), el límite puede estar muy por debajo de 1/100 de Hz a 44,1 kHz sin ningún problema.

SI tiene problemas de estabilidad, la FFT tampoco ayuda mucho. Como su filtro es un filtro IIR, la respuesta al impulso es infinita y primero debería truncarse. A estas frecuencias muy bajas, la respuesta al impulso se prolonga tanto que la FFT tampoco resulta práctica.

Por ejemplo, si desea un límite de 1/100 Hz a 44,1 kHz y desea un rango dinámico en la respuesta al impulso de 100 dB, ¡necesita aproximadamente 25 millones de muestras! Eso es casi 10 minutos a 44.1 kHz y muchas, muchas veces más que su señal original

Hilmar
fuente
Esto realmente no responde la pregunta sobre problemas numéricos, pero no estaba al tanto de los problemas con filter- ¡gracias! Mi corte de paso alto es de 0,5 Hz a 250 Hz. ¿Cuál es la razón de los problemas con filter? Estoy escribiendo la implementación yo mismo.
Andreas
2

¿Por qué supones que las cosas serán diferentes? Los conceptos teóricos deberían traducirse en aplicaciones prácticas, con la única diferencia de los problemas de coma flotante, de los que no podemos escapar. Puede verificar fácilmente con un ejemplo simple en MATLAB:

x=randn(5,1);
y=randn(5,1);
X=fft(x,length(x)+length(y)-1);
Y=fft(y,length(x)+length(y)-1);

z1=conv(x,y);z2=ifft(X.*Y);
z1-z2

ans =

   1.0e-15 *

   -0.4441
   -0.6661
         0
   -0.2220
    0.8882
   -0.2220
         0
   -0.4441
    0.8882

Como puede ver, los errores son del orden de precisión de la máquina. No debería haber ninguna razón para no usar el método FFT. Como mencionó Endolith, es mucho más común usar el enfoque FFT para filtrar / convolucionar / correlacionar, etc. y mucho más rápido, excepto para muestras muy pequeñas (como en este ejemplo). No es que el procesamiento en el dominio del tiempo nunca se haga ... todo se reduce a la aplicación, necesidades y restricciones.

Lorem Ipsum
fuente
1
Creo que la pregunta original fue profundizar en los problemas de punto flotante inherentes al filtrado basado en FFT versus la implementación directa del filtro en el dominio del tiempo. Esto puede ser una preocupación real para el procesamiento de señal de punto fijo, si tiene filtros realmente largos, o si tiene una mala implementación de FFT, por ejemplo. Definitivamente no verá ningún efecto para secuencias de longitud 5 en coma flotante de doble precisión.
Jason R
@JasonR Los errores siguen siendo de precisión de máquina si extiende la longitud de las secuencias a 1e6 en el ejemplo anterior. Los errores que menciona surgen principalmente debido a un diseño de filtro deficiente o una mala implementación de FFT. Si están bien, no veo por qué la convolución en el dominio del tiempo debería dar una respuesta diferente de la del dominio de la frecuencia.
Lorem Ipsum
1
No debería dar una respuesta diferente según el dominio en el que realice los cálculos. Mi único punto es que las operaciones matemáticas reales difieren mucho entre los dos enfoques, por lo que dependiendo de las implementaciones de cada uno que tenga disponible, es posible que usted Podía ver diferencias tangibles. Para una precisión doble, suponiendo que tenga buenas implementaciones, no esperaría ninguna diferencia, excepto en casos extremos.
Jason R