% program prodwindowcircBP.m
% Window design of circular 2-D band-pass filters

clear,clc,close all,home

% cd c:\jl\4F8_imageprocessing\imageprocessing4\imageproc4mat


% load Lenna for use later on
load lenna256
% x now contains grayscale Lenna

FL=input('Lower normalised cut-off frequency = ');
FU=input('Upper normalised cut-off frequency = ');
% use FL=0 and FU=0.25 for LP filters
% use FL=0.1 FU=0.25 for BP filters

N=input('Number of samples in single-sided impulse response = ');
% use N=7 for a 15x15 support
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Choose 1-D window function');
disp('1. Rectangular');
disp('2. Hamming');
disp('3. Kaiser');
option=input('option = ');

U=N;
if option==1
	w=boxcar(2*U+1);
	string1=['Impulse response: Rectangular window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
	string2=['Frequency response: Rectangular window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
elseif option==2
	w=hamming(2*U+1);
	string1=['Impulse response: Hamming window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
	string2=['Frequency response: Hamming window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
elseif option==3
%	beta=input('Sidelobe attenuation Kaiser window (dB) = ');
%	beta=0.1102*(beta-8.7);
	beta=input('beta for Kaiser window = ');
	w=kaiser(2*U+1,beta);
	string1=['Impulse response: Kaiser window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
	string2=['Frequency response: Kaiser window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];

end

% Form 2-D window
w2=w*w';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FL=2*pi*FL;
FU=2*pi*FU;

n=(-N:N)+1e-8;
h=((n.^2)'*ones(1,2*N+1)+ones(2*N+1,1)*(n.^2)).^(1/2);
h=(FU/(2*pi)*besselj(1,FU.*h)-FL/(2*pi)*besselj(1,FL.*h))./h;
h=h.*w2;
% Plot impulse response
mesh(h)
title(string1);
rotate3d on;
pause

% Plot frequency response
pad=64;
H=abs(fftshift(fft2(h,pad,pad)));
f=(-pad/2:pad/2-1)./pad;
mesh(f,f,H)
title(string2);
rotate3d on;
pause;

% now filter Lenna with this filter to see the effect:
% for best effects use FU =0.1 or so.

figure(2);
LPz=filter2(h,x);
subplot(1,2,1),grayimage(x);
title('Original 256x256');
zoom on
subplot(1,2,2),grayimage(LPz);
zoom on
rotate3d on;
pause;


