% Frequency response of an FIR filter, example 3.4
h = [1 1 1 1 1; 1 2 2 2 1; 1 2 3 2 1; 1 2 2 2 1; 1 1 1 1 1]/35;
[H,u,v] = freqz2(h,64,64,1);
[U V] = meshgrid(u,v);
%perspective plot of Frequency response
mesh(U,V,abs(H))
colormap([0 0 0]); 
set(gca,'ydir','reverse'); 
xlabel('horizontal frequency uX (c/sample)');
ylabel('vertical frequency vX (c/sample)')
print -deps 'FIRfilt_examp2_resp_persp.eps' %save as eps file when needed
%contour plot of Frequency response
vc=0:.1:1.;             % contours will be from 0 to 1 in steps of 0.1
figure, [C,hc]=contour(U,V,abs(H),vc); % generate the contour plot, including values to label contours
axis square
clabel(C,hc)             %label the contours
xlabel('Normalized horizontal frequency uX (c/sample)');
ylabel('Normalized vertical frequency vX (c/sample)');
set(gca,'ydir','reverse'); 
set(gca,'XAxisLocation','top');      %Xaxis labels on top
colormap([0 0 0]);     % use black only
%generate a zoneplate
dx=1/256;                     %sample spacing
x=-0.5+dx:dx:0.5; y=x;          %x and y values between -0.5 and +0.5 spaced by dx
[X,Y]=meshgrid(x,y);         % generate a 2D grid of xy values
r2=1.*dx;
Z=(1+cos(pi*(X.^2+Y.^2)./r2))./2;   % generate the zone plate
figure, imshow(Z)
imwrite(Z,'zoneplate.jpg','quality',100)
imwrite(Z,'zoneplate.tif')
B = conv2(Z,h,'same');
figure,imshow(B)
imwrite(B,'zoneplate_filtered.jpg','quality',100)
%filter Barbara image
Z2 = im2double(imread('barb.tif'));
figure, imshow(Z2)
B2 = conv2(Z2,h,'same');
figure, imshow(B2)