Python: 图像的带通滤波器

Question

我有一张数据图像，其中存在一个正弦形式的成像伪影，我希望去除。由于它是单频正弦波，傅里叶变换并且进行带通滤波或“陷波滤波”似乎很自然（在这里我会在正负omega处使用高斯滤波器）。

在尝试解决这个问题时，我注意到两件事：

1）仅通过执行FFT和反变换，我已经减少了正弦波分量，如下所示。似乎只是通过来回转换就对数据进行了一些高通滤波？

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)

这是img_back的图像：

也许这里的滤波已经对我足够好了，但我对背景抑制没有很好的理解，所以不太自信。

2）为了更确保在非期望频率上的抑制，我制作了一个布尔型“带通”掩模并将其应用于数据，但傅里叶变换会忽略该掩模。

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)

由于np.fft忽略掩码，因此结果与之前相同。解决方法很简单:

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

不幸的是(但仔细想想也不足为奇)，结果如下所示:

左图显示了应用带通滤波器后的FFT振幅。它是中心（DC）分量周围的黑暗环。相位未显示。

显然，“砖墙”滤波器不是正确的解决方案。这种滤波器产生环的现象在这里得到了很好的解释：当你将砖墙滤波器应用于1D数据集时会发生什么。

所以现在我陷入了困境。也许最好使用其中一个内置的scipy方法，但它们似乎是针对1d数据的，就像这个butterworth滤波器的实现。也许正确的做法涉及使用fftconvolve()，就像这里用于模糊图像。我对fftconvolve的问题是:它是否需要两个'图像'(图像和过滤器)都在实空间中？我认为是的，但在示例中他们使用高斯，所以不确定(fft(gaussian)=gaussian)。如果是这样，那么尝试制作一个实空间带通滤波器似乎是错误的。也许正确的策略是使用convolve2d()与傅里叶空间图像和自制滤波器。如果是这样，你知道如何制作一个好的2d滤波器吗？

- Claire Thomas

---

尝试这样做：filtdat2 = filtdat.filled(0)，然后对filtdat2执行ifft。 - HYRY

是的，我可以做到这一点，并且它确实过滤掉了不需要的正弦波，但代价是数据的大量损失。这个滤波器被称为“砖墙”滤波器。当我创建的环中出现硬边缘时，进行ifft会产生类似于空气盘的东西。这里解释得非常好：为什么不要应用砖墙滤波器。 - Claire Thomas

Answer 1

因此，一个问题在于你的背景正弦波周期与你试图保留的信号分量并没有太大区别。也就是说，信号峰值的间距与背景的周期大致相同。这将使过滤变得困难。

我的第一个问题是，这个背景是否真的从实验到实验都是恒定的，还是取决于样品和实验设置？如果它是恒定的，那么背景帧减法比过滤更有效。

大多数标准的scipy.signal滤波函数（bessel、chebychev等）都是针对1-D数据设计的，但你可以很容易地将它们扩展到2-D的各向同性滤波。频率空间中的每个滤波器都是f的有理函数。两种表示形式是[a,b]，它们是分子和分母多项式的系数，或者[z,p,k]，它是多项式的因式分解表示，即：H(f) = k(f-z0)*(f-z1)/(f-p0)*(f-p1) 你只需从其中一个滤波器设计算法中获取多项式，将其作为sqrt（x^2+y^2）的函数求值，并将其应用于频域数据。

你能贴出原始图像数据的链接吗？

- Evan