2017-07-17 150 views
0

据页#14 this link,方程为一个高通巴特沃斯滤波器是,如何在Matlab中实现高通巴特沃斯滤波器?

enter image description here

而且,根据页#17,输出应该是像下面这样,

enter image description here

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

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

enter image description here

什么是我的源代码可能出现的问题?


源代码

的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); 
+0

@beaker,这个答案是不正确。我也提供了链接。 – anonymous

+0

我已经修正了我的帖子中的巴特沃斯HPF方程,并且解释了为什么需要'sqrt(2) - 1'的因子。 – rayryeng

回答

0

我已经解决了这个问题。

enter image description here

的关键,解决这个问题是功能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