How can I implement a high-pass Butterworth filter in Matlab?

2.1k views Asked by At

According to Page#14 of this link, the equation for a high-pass Butterworth filter is,

enter image description here

And, according to Page#17, the output should be something like the following,

enter image description here

Now, I have looked at this answer in SO, and has written the following Matlab code using the formula given in the linked pdf document.

The output looks different than that of the one given above.

enter image description here

What is the possible problem in my source code?


Source Code

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);
1

There are 1 answers

0
user366312 On BEST ANSWER

I have solved this problem.

enter image description here

The key to solving this problem was the function ifftshow().


Source Code

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