% Chapter 10, Fig. 10.15, 10.16
%Eric Dubois, updated 2018-12-20
%Window design, example of chapter 10
%Compute the ideal response and multiply by a circularly symmetric Hamming
%window
clear all, close all;
L=20; %the filter will be of size (2L+1)*(2L+1)
%Compute the analytical ideal impulse response
%W = 0.171; %use for L=6
%W = 0.1498; %use for L=10
W = 0.1366; %use for L=20
Npt = 128; % extent of the ideal impulse response 
xI = -Npt/2:Npt/2 -1; yI = xI;
[XI,YI] = meshgrid(xI+.0000001,yI+.0000001);  %offset grid to avoid divide by 0 at (0,0)
hI = (W./sqrt(XI.^2 + YI.^2)).*besselj(1,2*pi*W*sqrt(XI.^2+YI.^2));
figure, mesh(XI,YI,hI);
colormap([0 0 0]);     % use black only
%Create a 2D circularly symmetric Hamming window
win = zeros(2*L+1,2*L+1); %the zero point corresponds to (L+1,L+1)
x = -L:L; y = x;
[X,Y] = meshgrid(x,y);
win =  0.54 + 0.46*cos(pi*sqrt((X).^2 + (Y).^2)/L);
win((X).^2 + (Y).^2 > L^2) =  0;
figure, mesh(X,Y,win);
colormap([0 0 0]);     % use black only
%the zero point corresponds to (Npt/2+1, Npt/2+1)
h = hI(Npt/2+1-L:Npt/2+1+L,Npt/2+1-L:Npt/2+1+L).*win;
DCgain = sum(sum(h));
h = h/DCgain;
%Compute and plot the frequency response
[H,fx,fy] = freqz2(h,64,64,[1 1]);
Hmag = abs(H);
[Fx,Fy] = meshgrid(fx,fy);
%contour plot
figure;
v=0:.1:1.;             % contours will be from 0 to 1 in steps of 0.2
[C,h1]=contour(Fx,Fy,abs(H),v); % generate the contour plot, including values to label contours
axis square
v1 = 0:.2:1.;
clabel(C,h1,v1,'FontSize',6);           %label the contours  
xlabel('\it u \rm (c/px)'), ylabel('\it v \rm (c/px)');
set(gca,'ydir','reverse'); 
colormap([0 0 0]);     % use black only
hold on
% plot a circle of radius 1/8
rectangle('Position',[-0.125, -0.125, 0.25, 0.25],'Curvature',[1,1],'EdgeColor','r','LineWidth',1);
axis square
set(gca,'ydir','reverse'); 
set(gca,'fontname','times')          %set the plot font to Times
colormap([0 0 0]);     % use black only
hold off
set(gcf,'Color',[1 1 1]);
%perspective plot
figure;
mesh(Fx(1:2:64,1:2:64),Fy(1:2:64,1:2:64),abs(H(1:2:64,1:2:64)));            % generate the perspective plot
colormap([0 0 0]);     % use black only
xlabel('\it u \rm (c/px)'), ylabel('\it v \rm (c/px)');
set(gca,'ydir','reverse'); 
set(gca,'fontname','times')          %set the plot font to Times
set(gcf,'Color',[1 1 1]);
%Filter the barbara image
A = im2double(imread('barb512_gray.tif'));
B = conv2(A,h,'same');
figure;
imshow(B)