Filtro spaziale Fourier

Question

Ho una serie temporale vettoriale reale x di lunghezza T e un filtro h di lunghezza t < < T. h è un filtro nello spazio di Fourier, reale e simmetrico. È approssimativamente 1/f.

Mi piacerebbe filtrare x con h per ottenere y.

Supponiamo che t == T e FFT di lunghezza T potrebbero essere inseriti nella memoria (nessuno dei quali è vero). Per ottenere le mie x filtrati in pitone, lo farei:

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

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

Dal momento che le condizioni non sono titolari, ho provato il seguente trucco:

1. Selezionare una dimensione imbottitura P < t/2, selezionare una dimensione di blocco B tale che B + 2P sia una buona dimensione FFT
2. Scala h tramite interpolazione spline per essere di dimensione B + 2P> t (h_scaled)
3. y = []; Ciclo: blocco
   • asporto di lunghezza B + 2P da x (chiamato x_b)
   • Eseguire y_b = IFFT (FFT (x_b) * h_scaled)
   • goccia imbottitura P da entrambi i lati di y_b e concatenare con y
   • Seleziona accanto x_b sovrapposizione con l'ultimo da P

Si tratta di una buona strategia? Come faccio a selezionare il padding P in un buon modo? Qual è il modo corretto per farlo? Non conosco molta elaborazione del segnale. Questa è una buona occasione per imparare.

Sto usando cuFFT per velocizzare le cose, quindi sarebbe bello se la maggior parte delle operazioni fossero FFT. Il vero problema è il 3D. Inoltre, non sono preoccupato per gli artefatti da un filtro acausale.

Grazie, Paul.

Answer 1

Sei sulla strada giusta. La tecnica è chiamata overlap-save processing. t è abbastanza corto da consentire la memorizzazione di FFT di quella lunghezza? Se è così, è possibile scegliere la dimensione del blocco B tale che B > 2*min(length(x),length(h)) e rende una trasformazione veloce. Quindi quando si elabora, si rilascia la prima metà di y_b anziché cadere da entrambe le estremità.

Per capire perché si perde la prima metà, ricordare che la moltiplicazione spettrale è la stessa della convoluzione circolare nel dominio del tempo. Il coinvolgimento dello h a zero-zero crea strani transienti glitch nella prima metà del risultato, ma entro la seconda metà tutti i transienti sono spariti perché il punto di avvolgimento circolare in x è allineato con la parte zero di h. C'è una buona spiegazione di questo, con le immagini, in "Theory and Application of Digital Signal Processing" by Lawrence Rabiner and Bernard Gold.

È importante che il filtro del dominio del tempo si assottigli a 0 almeno su un'estremità in modo da non generare artefatti che squillano. Lei dice che h è reale nel dominio della frequenza, il che implica che ha tutta la fase 0. Solitamente, tale segnale sarà continuo solo in modo circolare e, se usato come filtro, creerà distorsioni per tutta la banda di frequenza. Un modo semplice per creare un filtro ragionevole è progettarlo nel dominio della frequenza con fase 0, trasformazione inversa e rotazione.Ad esempio:

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 

(Sono abbastanza nuovo in Python, per favore fatemi sapere se c'è un modo migliore per codificarlo).

Se si confrontano le risposte in frequenza di ht, htrot, e htwin, viene visualizzato il seguente (asse x è frequenza normalizzata fino a pi):

ht, in alto, ha un sacco di ripple . Ciò è dovuto alla discontinuità al limite. htrot, nel mezzo, è meglio, ma ha ancora increspature. htwin è bello e scorrevole, a scapito dello spianamento a una frequenza leggermente più alta. Notare che è possibile estendere la lunghezza della sezione rettilinea utilizzando un valore maggiore per N.

Ho scritto sul problema della discontinuità e ho anche scritto un esempio Matlab/Octave in another SO question se si desidera visualizzare maggiori dettagli.