Frequency-Domain Bandreject Filter(Butterworth)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Title : Bandreject Filter(Butterworth) %
% Domain: Frequency %
% Author: S.Ganesh Babu %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all ;clc;
f = im2double(imread('test_pattern.tif'));
figure,imshow(f),title('Original Image');
PQ = paddedsize(size(f)); % Get the size of the Image
% Pad the image with Zeros
W=50; %------------------------- Width of the band
D0 = 230; %------------------------- Cutoff Frequency
%----------------------------------
u = 0:(PQ(1)-1); %|
v = 0:(PQ(2)-1); %|
idx = find(u > PQ(1)/2); %| Meshgrid frequency matrices
u(idx) = u(idx) - PQ(1); %| for use in distance computation
idy = find(v > PQ(2)/2); %|
v(idy) = v(idy) - PQ(2); %|
[V, U] = meshgrid(v, u); %|
%-------------------------------------------
D = sqrt(U.^2 + V.^2); % ------------------ Computing the Distance Matrix
%-------------------------------------------
n=1; %------------------------------------- Order of the filter
H = 1./(1 + (((D.*W)./((D.*D)-(D0*D0))).^(2*n))); % Bandreject Butterworth Filter
H1 = fftshift(H); % Fourier version of Filter
%----------------------------------------
figure,surf(H1(1:10:PQ(1), 1:10:PQ(2))) %|
colormap([0 1 0.5]); %|Display the Surface plot of Filter
axis off %|
figure,imshow(H1,[]),title('Filter'); %|
%--------------------------------------------
F = fft2(f, size(H, 1), size(H, 2)); % Fourier version of Image
g = real(ifft2(H.*F)); % Multiplying filter Co-efficient with Image
% then takes the Inverse Fourier of Real Part
g = g(1:size(f, 1), 1:size(f, 2)); % Make the image in Original Size
%--------------------------------------------
figure, imshow(g,[]),title('Band Rejected Image');

function PQ = paddedsize(AB, CD, PARAM)
if nargin == 1
PQ = 2*AB;
elseif nargin == 2 & ~ischar(CD)
PQ =AB + CD - 1;
PQ = 2 * ceil(PQ / 2);
elseif nargin == 2
m = max(AB);
p = 2^nextpow2(2*m);
PQ = [p, p];
elseif nargin == 3
m = max([AB CD]);
p = 2^nextpow2(2*m);
PQ = [p, p];
else
error('Wrong number of inputs.');
end