python: bandpass filter of an image
Question: I have a data image with an imaging artifact that comes out as a sinusoidal
background, which I want to remove. Since it is a single frequency sine wave,
it seems natural to Fourier transform and either bandpass filter or "notch
filter" (where I think I'd use a gaussian filter at +-omega).

[![My data.  The red spots are what I want, the background sine wave in kx+ky
is
unwanted.](http://i.stack.imgur.com/AacAP.png)](http://i.stack.imgur.com/AacAP.png)

In trying to do this, I notice two things:

1) simply by performing the fft and back, I have reduced the sine wave
component, shown below. There seems to be some high-pass filtering of the data
just by going there and back?

    
    
    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)
    

This is an image of img_back:

[![The inverse fourier transform, no filter
applied.](http://i.stack.imgur.com/DKCEi.png)](http://i.stack.imgur.com/DKCEi.png)

Maybe the filtering here is good enough for me, but I'm not that confident in
it since I don't have a good understanding of the background suppression.

2) To be more sure of the suppression at the unwanted frequencies, I made a
boolean 'bandpass' mask and applied it to the data, but the fourier transform
ignores the mask.

    
    
    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)
    

Here the result is the same as before, because np.fft ignores masks. The fix
to that was simple:

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

Unfortunately, (but upon reflection also unsurprisingly) the result is shown
here:

[![The result of a brickwall bandpass
filter.](http://i.stack.imgur.com/9Klkd.png)](http://i.stack.imgur.com/9Klkd.png)

The left plot is of the amplitude of the FFT, with the bandpass filter
applied. It is the dark ring around the central (DC) component. The phase is
not shown.

Clearly, the 'brickwall' filter is not the right solution. The phenomenon of
making rings from this filter is well explained here: [What happens when you
apply a brick-wall filter to a 1D
dataset.](http://dsp.stackexchange.com/questions/724/low-pass-filter-and-fft-
for-beginners-with-python/725#725)

So now I'm stuck. Perhaps it would be better to use one of the built in scipy
methods, but they seem to be for 1d data, as [in this implementation of a
butterworth filter](http://wiki.scipy.org/Cookbook/ButterworthBandpass).
Possibly the right thing to do involves using fftconvolve() as is done [here
to blur an
image.](http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html#scipy.signal.fftconvolve)
My question about fftconvolve is this: Does it require both 'images' (the
image and the filter) to be in real space? I think yes, but in the example
they use a gaussian, so it's ambiguous (fft(gaussian)=gaussian). If so, then
it seems wrong to try to make a real space bandpass filter. Maybe the right
strategy uses convolve2d() with the fourier space image and a homemade filter.
If so, do you know how to make a good 2d filter?


Answer: So, one problem here is that your background sinusoid has a period not
terribly different from the signal components you are trying to preserve.
i.e., the spacing of the signal peaks is about the same as the period of the
background. This is going to make filtering difficult.

My first question is whether this background is truly constant from experiment
to experiment, or does it depend on the sample and experimental setup? If it
is constant, then background frame subtraction would work better than
filtering.

Most of the standard scipy.signal filter functions (bessel, chebychev, etc.)
are, as you say, designed for 1-D data. But you can easily extend them to
isotropic filtering in 2-D. Each filter in frequency space is a rational
function of f. The two representations are [a,b] which are the coefficiets of
the numerator and denominator polynomial, or [z,p,k] which is the factored
representation of the polynomial i.e.,: `H(f) = k(f-z0)*(f-z1)/(f-p0)*(f-p1)`
You can just take the polynomial from one of the filter design algorithms,
evaluate it as a function of sqrt(x^2+y^2) and apply it to your frequency
domain data.

Can you post a link to the original image data?

