2015-08-11 23 views
16

Ho un'immagine di dati con un artefatto di immagini che viene fuori come sfondo sinusoidale, che voglio rimuovere. Poiché si tratta di un'onda sinusoidale a frequenza singola, sembra naturale la trasformata di Fourier e il filtro passa-banda o "notch filter" (dove penso che utilizzerei un filtro gaussiano a + -omega).python: filtro passa-banda di un'immagine

My data. The red spots are what I want, the background sine wave in kx+ky is unwanted.

Nel tentativo di fare questo, ho notato due cose:

1) semplicemente eseguendo la FFT e ritorno, ho ridotto la componente onda sinusoidale, illustrato di seguito. Sembra che ci sia un filtro passa-alto dei dati solo andando avanti e indietro?

import numpy as np 

f = np.fft.fft2(img)     #do the fourier transform 
fshift1 = np.fft.fftshift(f)   #shift the zero to the center 
f_ishift = np.fft.ifftshift(fshift1) #inverse shift 
img_back = np.fft.ifft2(f_ishift)  #inverse fourier transform 
img_back = np.abs(img_back) 

Questa è l'immagine di img_back:

The inverse fourier transform, no filter applied.

Forse il filtraggio qui è abbastanza buono per me, ma io non sono così sicuro in essa dal momento che non ho una buona comprensione della soppressione dello sfondo.

2) Per essere più sicuro della soppressione alle frequenze indesiderate, ho creato una maschera "passabanda" booleana e l'ho applicata ai dati, ma la trasformata di Fourier ignora la maschera.

a = shape(fshift1)[0] 
b = shape(fshift1)[1] 

ro = 8 
ri = 5 
y,x = np.ogrid[-a/2:a/2, -b/2:b/2] 
m1 = x*x + y*y >= ro*ro 
m2 = x*x + y*y <= ri*ri 
m3=np.dstack((m1,m2))  
maskcomb =[] 
for r in m3: 
    maskcomb.append([any(c) for c in r]) #probably not pythonic, sorry 
newma = np.invert(maskcomb) 
filtdat = ma.array(fshift1,mask=newma) 
imshow(abs(filtdat)) 
f_ishift = np.fft.ifftshift(filtdat) 
img_back2 = np.fft.ifft2(f_ishift) 
img_back2 = np.abs(img_back2) 

Qui il risultato è lo stesso di prima, perché np.fft ignora le maschere. La correzione di che era semplice:

filtdat2 = filtdat.filled(filtdat.mean())

Sfortunatamente, (ma riflettendoci anche sorprendente) il risultato è mostrato qui:

The result of a brickwall bandpass filter.

La trama di sinistra è l'ampiezza della FFT , con il filtro passa-banda applicato. È l'anello scuro attorno al componente centrale (DC). La fase non è mostrata.

Chiaramente, il filtro "brickwall" non è la soluzione giusta. Il fenomeno della creazione di anelli da questo filtro è ben spiegato qui: What happens when you apply a brick-wall filter to a 1D dataset.

Quindi ora sono bloccato. Forse sarebbe meglio usare uno dei metodi scipy incorporati, ma sembrano essere per dati 1d, come in this implementation of a butterworth filter. Forse la cosa giusta da fare è usare fftconvolve() come è fatto here to blur an image. La mia domanda su fftconvolve è questa: Richiede che entrambe le 'immagini' (l'immagine e il filtro) siano nello spazio reale? Penso di sì, ma nell'esempio usano un gaussiano, quindi è ambiguo (fft (gaussiano) = gaussiano). Se è così, allora sembra sbagliato cercare di creare un vero filtro passa-banda spaziale. Forse la strategia giusta usa convolve2d() con l'immagine spaziale quadrupla e un filtro fatto in casa. Se è così, sai come creare un buon filtro 2D?

+0

provare questo: 'filtdat2 = filtdat.filled (0)', quindi eseguire ifft su 'filtdat2'. – HYRY

+0

Sì, posso farlo e filtra l'onda sinusoidale indesiderata, ma a grandi spese per i dati. Questo filtro è chiamato filtro "brick-wall". I bordi duri dell'anello che ho creato, quando faccio ifft, si traducono in qualcosa che assomiglia a dischi aerati. È spiegato molto bene qui: [perché non applicare i filtri brick-wall] (http://dsp.stackexchange.com/questions/724/low-pass-filter-and-fft-for-beginners-with-python/ 725 # 725) –

risposta

1

Quindi, un problema qui è che il sinusoide di sfondo ha un periodo non molto diverso dai componenti del segnale che si sta tentando di preservare. cioè, la spaziatura dei picchi del segnale è all'incirca uguale al periodo dello sfondo. Questo renderà il filtraggio difficile.

La mia prima domanda è se questo sfondo è veramente costante dall'esperimento all'esperimento, o dipende dall'esempio e dall'impostazione sperimentale? Se è costante, la sottrazione del frame in background funzionerebbe meglio del filtro.

La maggior parte delle funzioni standard del filtro scipy.signal (bessel, chebychev, ecc.) Sono, come dici tu, progettate per i dati 1-D. Ma puoi facilmente estenderli al filtro isotropico in 2D. Ogni filtro nello spazio di frequenza è una funzione razionale di f. Le due rappresentazioni sono [a, b] che sono i coefficienti del polinomio del numeratore e denominatore, o [z, p, k] che è la rappresentazione fattoriale del polinomio cioè, H(f) = k(f-z0)*(f-z1)/(f-p0)*(f-p1) Puoi semplicemente prendere il polinomio da uno dei algoritmi di progettazione del filtro, valutarlo come una funzione di sqrt (x^2 + y^2) e applicarlo ai dati del dominio della frequenza.

È possibile pubblicare un collegamento ai dati dell'immagine originale?