% Chapter 3, Example 3.4, Fig. 3.10-3.12
% Eric Dubois, updated 2018-11-14
% Frequency response of a 5 x 5 FIR filter and filter output
% In this example,
% $X=1$ and 
%
%
% Generate the unit sample response and plot its frequency response
%
clear all, close all;
%load the unit sample response into matrix h
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;
%compute the frequency response using freqz2, spacing between samples is X=1
[H,u,v] = freqz2(h,64,64,1);
[U V] = meshgrid(u,v);

figure

%perspective plot of the magnitude of the Frequency response
mesh(U,V,abs(H))
colormap([0 0 0]); 
set(gca,'ydir','reverse'); 
xlabel('\itu \rm(c/px)');
ylabel('\itv \rm(c/px)')
set(gca,'FontName','times')          %set the plot font to times
set(gca,'FontSize',8)                %set the plot size to 8pt
set(gcf,'units','centimeters','position',[0.5,0.5,15,15]);
set(gcf,'Color',[1 1 1]);

%contour plot of Frequency response
figure

vc=0:.1:1.;             % contours will be from 0 to 1 in steps of 0.1
% generate the contour plot, including values to label contours
[C,hc]=contour(U,V,abs(H),vc); 
axis square
clabel(C,hc,'FontSize',6)             %label the contours
xlabel('\itu \rm(c/px)');
ylabel('\itv \rm(c/px)');
set(gca,'ydir','reverse'); 
set(gca,'XAxisLocation','top');      %Xaxis labels on top
colormap([0 0 0]);     % use black only
set(gca,'FontName','times')          %set the plot font to times
set(gca,'FontSize',8)                %set the plot size to 8pt
set(gcf,'units','centimeters','position',[0.5,0.5,15,15]);
set(gcf,'Color',[1 1 1]);

% Generate a zoneplate and filter it 
%generate a zoneplate
dx=1/512;                     %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
r0=.05;
Z=(1+cos(pi*(X.^2+Y.^2)./(r0^2)))./2;   % generate the zone plate
figure, imshow(Z)

B = conv2(Z,h,'same');
figure,imshow(B)

% Filter the Barbara image with the filter
%filter Barbara image
Z2 = im2double(imread('barb512_gray.tif'));
figure, imshow(Z2)

B2 = conv2(Z2,h,'same');
figure, imshow(B2)