Butterworth低通滤波器

Butterworth LPF

---

公式

利用一下方程得到一个和图像相同尺寸的滤波器$H$，进行滤波：

$$H(u,v)=\frac{1}{1+(\frac{D(u,v)}{D_0}{)}^{2n}}$$

其中$D_0$是一个常数，我们称为截止频率
D(u, v)定义如下:

$$D(u,v)=[(u-\frac{M}{2}{)}^2+(v-\frac{N}{2}{)}^2{]}^{\frac{1}{2}}$$

即在频率域中的点距离中心的距离。(通过对空域做中心变换后，我们在经过快速傅里叶变换之后得到的图像的中心将会是$(\frac{M}{2},\frac{N}{2})$，中心变换操作可以通过和$(-1{)}^{x+y}$做相关得到，可以通过傅里叶变换的平移性得到)

matlab步骤

（1）以$(-1{)}^{(x+y)}$ 乘以输入图像进行中心变换；

% uint8
f = imread('your_image_file');
% double
f = im2double(f);
[X, Y] = meshgrid(1:width, 1:height);
cent = f.*(-1).^(X + Y);

（2）直接以FFT2进行傅立叶变换；

F = fft2(cent);

（3）进行Butterworth滤波

% generate Butterworth filtering kernel
kernel = zeros(height, width);
for u = 1:height
    for v = 1:width
        D = sqrt((u - height / 2)^2 + (v - width / 2)^2);
        kernel(u, v) = 1 / (1 + (D / D0)^2);
    end
end
% fitering
F = F.*kernel;

（4）DFT反变换后取实部；

% get conjugate F
for u = 1:height
    for v = 1:width
        F(u, v) = (real(F(u, v)) - imag(F(u, v)) * 1i) / (height * width);
    end
end
% get conjugate f by doing 2-dimension FFT on conjugate F
cent = fft2(F);
% get real part
for u = 1:height
    for v = 1:width
        cent(u, v) = real(cent(u, v));
    end
end

（5）以$(-1{)}^{(x+y)}$ 乘以（3）中结果，反中心变换。

f = cent.* (-1).^(X + Y);

滤波效果

从频域的角度来看，观察我们用来构造滤波器的公式，当维度$n$保持不变时，随着截止频率$D_0$的增大，分母越小，分母的增加速度也越小，也就是滤波器从中心往周围扩散的下降趋势变小了，变得较为平坦，我们可以从下图图像的滤波器的变化可以看出，滤波器亮度高的部分逐渐展开变宽。这也意味着更多的高频成分被保留了，图像的细节也就越明显了。

从空域的角度来看，考虑高斯低通滤波器，再考虑高斯滤波器的傅里叶变换，高斯滤波器的傅里叶变换之后还是一个高斯滤波器，那么显然，高斯滤波器傅里叶反变换之后的结果也是一个高斯滤波器，也就是在空域中。Butterworth低通滤波器结构与效果和高斯低通滤波器相似，那么，当我们对Butterworth反变换与图像在空域中进行卷积，也将会类似在空域中进行高斯滤波后的效果，即对图片进行平滑操作，图像会变得模糊。参考一下空域卷积和到频域的转换式子：

$$f(x,y)\ast g(x,y)⟷F(x,y)\cdot G(x,y)$$

即，空域的卷积相当于频域率的逐点相乘的结果，反过来也就是说，在频域中进行Butterworth滤波，相当于对Butterworth滤波器反变换到空域中和图像做卷积。