2010-09-23 9 views
5

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.Filtro spaziale Fourier

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.

risposta

6

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): 1/f filters with length 64

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.

+0

Grazie per il riferimento di sovrapposizione. Ne avevo letto su Press et al., Numerical Recipes, per quanto riguarda il filtraggio del dominio del tempo e non ero sicuro di come mapparlo al dominio della frequenza. Non sono sicuro di lasciar perdere: 1) perché la seconda metà di y_b piuttosto che le estremità, 2) nel tuo altro post SO, si rilascia la prima metà. – Paul

+0

Il mio filtro h è derivato da una media su dati non elaborati, con h (f) ~ 1/f e fasi impostate su 0. Sto filtrando un segnale sintetico con questo filtro per conferirgli uno spettro più simile ai miei dati grezzi. Non sono sicuro di come progettare questo filtro nel dominio del tempo. Una cosa che hai sottolineato è che ifft (h) meglio essere zero ad una estremità per evitare artefatti che suonano. Lo controllerò perché è molto probabile che no. Esiste un analogo all'applicazione di una finestra di Hamming nel dominio del tempo a qualche metodo di finestra nel dominio della frequenza (il tuo primo esempio nel tuo altro post SO)? – Paul

+0

Yeesh - mi dispiace per aver rovinato la prima metà/seconda metà del numero. Ho aggiornato con quella correzione e alcuni pensieri su come generare un 'h' ben educato. – mtrw