如何在 Matlab 中实现低通巴特沃斯滤波器？

Question

来自this answer，我知道如何创建高通巴特沃斯滤波器。

从this video，我知道，lowpasskernel = 1 - highpasskernel。

所以，我创建了以下低通巴特沃斯滤波器，

function [out, kernel] = butterworth_lp(I, Dl, n)
    height = size(I, 1);
    width  = size(I, 2);

    [u, v] = meshgrid(-floor(width/2):floor(width/2)-1,-floor(height/2):floor(height/2)-1);

    % lp_kernel = 1 - hp_kernel
    kernel = 1 - butter_hp_kernel(u, v, Dl, n);

    % fft the image
    I_fft_shifted = fftshift(fft2(double(I)));

    % apply lowpass filter
    I_fft_shift_filtered = I_fft_shifted .* kernel;

    % inverse FFT, get real components
    out = real(ifft2(ifftshift(I_fft_shift_filtered)));

    % normalize and cast
    out = (out - min(out(:))) / (max(out(:)) - min(out(:)));
    out = uint8(255*out);

function k = butter_hp_kernel(u, v, Dh, n)
    uv =  u.^2+v.^2;
    D = sqrt(uv);
    frac = Dh./D;
    p = 2*n;
    denom =  frac.^p;
    k = 1./denom;

输出

这不是低通滤波器输出。

那么，我的代码有什么问题？

Answer 1

好的。我已经按照following formula (Page #8/48)解决了这个问题，

输出

源代码

butter_lp_kernel.m

function f = butter_lp_f(u, v, Dl, n)
    uv = u.^2+v.^2;
    Duv = sqrt(uv);
    frac = Duv./Dl;
    denom = frac.^(2*n);
    f = 1./(1.+denom);

function k = butter_lp_kernel(I, Dl, n) 
    Height = size(I,1); 
    Width = size(I,2); 

    [u, v] = meshgrid( ...
                    -floor(Width/2) :floor(Width-1)/2, ...
                    -floor(Height/2): floor(Height-1)/2 ...
                 ); 

    k = butter_lp_f(u, v, Dl, n);

ifftshow.m

function out = ifftshow(f)
    f1 = abs(f);
    fm = max(f1(:));
    out = f1/fm;
end

butterworth_lpf.m

function [out1, out2] = butterworth_lpf(I, Dl, n)

    Kernel = butter_lp_kernel(I, Dl, n);

    I_ffted_shifted = fftshift(fft2(I));

    I_ffted_shifted_filtered = I_ffted_shifted.*Kernel;

    out1 = ifftshow(ifft2(I_ffted_shifted_filtered));

    out2 = ifft2(ifftshift(I_ffted_shifted_filtered));
end

Answer 2

您没有正确复制高通滤波器中的公式：

 uv =  u.^2+v.^2;
 D = uv.^0.5;
 frac = Dh./D;
 p = 2*n;
 denom =  frac.^p;
 k = 1./(1+denom); --> you forgot the "1+"

新输出：