Filtrado de espacio de Fourier

Question

Tengo una serie de tiempo vector real x de longitud T y un filtro h de longitud t < < T. h es un filtro en el espacio de Fourier, real y simétrico. Es aproximadamente 1/f.

Me gustaría filtrar x con h para obtener y.

Supongamos que t == T y FFT de longitud T pueden caber en la memoria (ninguno de los cuales es verdadero). Para conseguir mis x filtrados en Python, lo haría:

import numpy as np 
from scipy.signal import fft, ifft 

y = np.real(np.ifft(np.fft(x) * h))) 

Dado que las condiciones no se cumplen, he intentado este truco:

1. seleccionar un tamaño de relleno P < t/2, seleccione un tamaño de bloque B tal que B + 2P es un buen tamaño de FFT
2. Escala h mediante la interpolación de spline para ser del tamaño B + 2P> t (h_scaled)
3. y = []; Loop: bloque
   • Take de longitud B + 2P de x (llamado x_b)
   • Realizar y_b = IFFT (FFT (x_b) * h_scaled)
   • gota relleno P desde cualquier lado de y_b y concatenar con Y
   • Seleccionar junto con la superposición de x_b pasado por P

¿Es esta una buena estrategia? ¿Cómo selecciono el relleno P de una buena manera? ¿Cuál es la forma apropiada de hacer esto? No sé mucho procesamiento de señal. Esta es una buena oportunidad para aprender.

Estoy usando cuFFT para acelerar las cosas, por lo que sería genial si la mayor parte de las operaciones son FFT. El problema real es 3D. Además, no me preocupan los artefactos de un filtro acausal.

Gracias, Paul.

Answer 1

Estás en el camino correcto. La técnica se llama overlap-save processing. ¿Es t lo suficientemente corto como para que las FFT de esa longitud quepan en la memoria? Si es así, puede elegir su tamaño de bloque B de manera que B > 2*min(length(x),length(h)) y realice una transformación rápida. Luego, cuando procesa, suelta la primera mitad de y_b, en lugar de caer desde ambos extremos.

Para ver por qué se cae la primera mitad, recuerde que la multiplicación espectral es la misma que la convolución circular en el dominio del tiempo. Convolverse con h con cero relleno crea transitorios glitchy raros en la primera mitad del resultado, pero en la segunda mitad todos los transitorios se han ido porque el punto de envoltura circular en x está alineado con la parte cero de h. Hay una buena explicación de esto, con imágenes, en "Theory and Application of Digital Signal Processing" by Lawrence Rabiner and Bernard Gold.

Es importante que su filtro de dominio de tiempo se reduzca a 0 como mínimo en un extremo para que no se produzcan artefactos de timbre. Menciona que h es real en el dominio de frecuencia, lo que implica que tiene toda la fase 0. Usualmente, tal señal será continua solo de forma circular, y cuando se usa como un filtro creará distorsión a través de la banda de frecuencia. Una manera fácil de crear un filtro razonable es diseñarlo en el dominio de frecuencia con fase 0, transformación inversa y rotación.Por ejemplo:

def OneOverF(N): 
    import numpy as np 
    N2 = N/2; #N has to be even! 
    x = np.hstack((1, np.arange(1, N2+1), np.arange(N2-1, 0, -1))) 
    hf = 1/(2*np.pi*x/N2) 
    ht = np.real(np.fft.ifft(hf)) # discard tiny imag part from numerical error 
    htrot = np.roll(ht, N2) 
    htwin = htrot * np.hamming(N) 
    return ht, htrot, htwin 

(Soy bastante nuevo en Python, por favor avíseme si hay una mejor manera de codificar esto).

Si se comparan las respuestas de frecuencia de ht, htrot y htwin, vea el siguiente (eje x es la frecuencia normalizada de hasta pi):

ht, en la parte superior, tiene un montón de ondulación . Esto se debe a la discontinuidad en el borde. htrot, en el medio, es mejor, pero todavía tiene ondas. htwin es agradable y suave, a expensas de acoplarse a una frecuencia ligeramente mayor. Tenga en cuenta que puede ampliar la longitud de la sección de línea recta utilizando un valor mayor para N.

Escribí sobre el problema de la discontinuidad y también escribí un ejemplo de Matlab/Octave en another SO question si desea ver más detalles.