如何在 Matlab 中实现高通巴特沃斯滤波器？ [复制]答案

【问题标题】：How can I implement a high-pass Butterworth filter in Matlab? [duplicate]如何在 Matlab 中实现高通巴特沃斯滤波器？ [复制]
【发布时间】：2017-12-22 06:50:50
【问题描述】：

根据this link 的第 14 页，高通巴特沃斯滤波器的方程是，

并且，根据第 17 页，输出应该类似于以下内容，

现在，我查看了这个answer in SO，并使用链接的pdf文档中给出的公式编写了以下Matlab代码。

输出看起来与上面给出的不同。

我的源代码中可能存在什么问题？

源代码

main.m

clear_all();

I = gray_imread('cameraman.png');

n = 1;
Dh = 10;

[J, Kernel] = butterworth_hp(I, Dh, n);

imshowpair(I, J, 'montage');

butterworth_hp.m

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

I_fft_shifted = fftshift(fft2(double(I)));

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

I_fft_shift_filtered = I_fft_shifted .* kernel;

out = real(ifft2(ifftshift(I_fft_shift_filtered)));

out = (out - min(out(:))) / (max(out(:)) - min(out(:)));
out = uint8(255*out);


function k = butter_hp_kernel(u, v, D0, n)
uv =  u.^2+v.^2;
Duv = uv.^0.5;
frac = D0./Duv;
p = 2*n;
denom =  frac.^p;
A = 0.414;
k = 1./(1+A*denom);

【问题讨论】：

@beaker，这个答案不正确。我还提供了链接。
我在您链接的帖子中更正了 Butterworth HPF 方程，并解释了为什么需要 sqrt(2) - 1 的因子。

标签： matlab highpass-filter butterworth

【解决方案1】：

我已经解决了这个问题。

解决这个问题的关键是函数ifftshow()。

源代码

main.m

clear_all();
I = gray_imread('cameraman.png');
Dh = 10;
n = 1;
[J, K] = butterworth_hpf(I, Dh, n);

imshowpair(I, K, 'montage');
%draw_multiple_images({I, J, K});

ifftshow.m

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

butter_hp_kernel.m

function k = butter_hp_kernel(I, Dh, 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_hp_f(u, v, Dh, n);

function f = butter_hp_f(u, v, Dh, n)
    uv = u.^2+v.^2;
    Duv = sqrt(uv);
    frac = Dh./Duv;
    %denom = frac.^(2*n);
    A=0.414; denom = A.*(frac.^(2*n));    
    f = 1./(1.+denom);

butterworth_hpf.m

function [out1, out2] = butterworth_hpf(I, Dh, n)

    Kernel = butter_hp_kernel(I, Dh, 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

【讨论】：