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

clear,clc,close all,home

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


FL=input('Lower normalised cut-off frequency = ');
FU=input('Upper normalised cut-off frequency = ');

% use FL=0 and FU=0.25 for LP filter
% use FL=0.1 FU=0.25 for BP filter

N=input('Number of samples in single-sided impulse response = ');
% use N=7 to reproduce figures in handout
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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=['Rotated rectangular window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
	string2=['Rotated rectangular window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
	w2=zeros(2*U+1,2*U+1);
	for n1=1:2*U+1
		for n2=1:2*U+1
			u=sqrt((n1-U)^2+(n2-U)^2);
			if u<=U
				w2(n1,n2)=1.0;
			end
		end
	end

elseif option==2
	string1=['Rotated Hamming window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
	string2=['Rotated Hamming window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
	w2=zeros(2*U+1,2*U+1);
	for n1=1:2*U+1
		for n2=1:2*U+1
			u=sqrt((n1-U)^2+(n2-U)^2);
			if u<=U
				w2(n1,n2)=0.54+0.46*cos(pi*u/U);
			end
		end
	end
elseif option==3
%	beta=input('Sidelobe attenuation Kaiser window (dB) = ');
%	beta=0.1102*(beta-8.7);
%       beta=3 for figures in handout
	beta=input('beta for Kaiser window = ');
	w=kaiser(2*U+1,beta);
	string1=['Rotated Kaiser window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
	string2=['Rotated Kaiser window, N= '...
	,num2str(2*U+1),'*',num2str(2*U+1)];
	w2=zeros(2*U+1,2*U+1);
	for n1=1:2*U+1
		for n2=1:2*U+1
			u=sqrt((n1-U)^2+(n2-U)^2);
			if u<=U
				w2(n1,n2)=besseli(0,beta*sqrt(1-(u/U)^2))/besseli(0,beta);
			end
		end
	end

end

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

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
load lenna256.mat
y=filter2(h,x);
grayimage(y)


