Search code examples
matlabbandpass-filterbutterworth

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

This is my expected output.

enter image description here

But, I am getting the following output.

enter image description here

What is wrong with my source code?

Source Code

main.m

clear_all();

I = gray_imread('woman.png');
d0 = 1;
d1 = 100;
n = 20;

J = butter_bandpass(I, d0, d1,n);

J = histeq(J);

K = {I, J};

imshowpair(I, J, 'montage');

butter_bandpass.m

function output_image = butter_bandpass(I, dL, dH,n)
    I_double = double(I);
    [Height, Width] = size(I_double);
    I_double = uint8(I_double);
    I_fft = fft2(I_double,2*Height-1,2*Width-1);
    I_fft = fftshift(I_fft);

    LowPass = ones(2*Height-1,2*Width-1);
    HighPass = ones(2*Height-1,2*Width-1);
    Kernel = ones(2*Height-1,2*Width-1);

    for i = 1:2*Height-1
        for j =1:2*Width-1
            D = ((i-(Height+1))^2 + (j-(Width+1))^2)^.5;        
            LowPass(i,j)= 1/(1 + (D/dH)^(2*n));
            HighPass(i,j)= 1 - (1/(1 + (D/dL)^(2*n)));

            % Create Butterworth filter.
            Kernel(i,j) = LowPass(i,j).*HighPass(i,j);
        end
    end

    % apply filter to image
    output_image = I_fft + Kernel.*I_fft;

    % revert the image to spatial domain
    output_image = ifftshift(output_image);
    output_image = ifft2(output_image,2*Height-1,2*Width-1);
    output_image = real(output_image(1:Height,1:Width));
    output_image = uint8(output_image);
end
Solution

  • I have solved this problem.

    Output

    enter image description here

    Refer to this low-pass-filter, and this high-pass-filter source codes.

    Since, we know that BPF(band-pass-filter) = LPF * HPF, we can implement a bandpass filter as follows,

    Source Code

    main.m

    clear_all();
I = gray_imread('woman.png');
Dl = 70;
Dh = 70;
n = 1;
[J, K] = butterworth_bpf(I, Dh, Dl, n);
J = histeq(J);
imshowpair(I, J, 'montage');

    butter_bp_kernel.m

    function k = butter_bp_kernel(I, Dh, Dl, n)
    hp_kernel = butter_hp_kernel(I, Dh, n);
    lp_kernel = butter_lp_kernel(I, Dl, n);    
    k = hp_kernel.*lp_kernel; % here is the multiplication

    butterworth_bpf.m

    function [out1, out2] = butterworth_bpf(I, Dh, Dl, n)
    Kernel = butter_bp_kernel(I, Dh, 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