% Chapter 12, Fig.12.3
% Eric Dubois, updated 2019-01-15
% Figure to illustrate the frequency response of an FIR Gaussian filter on a
% hexagonal lattice. 
% Required functions:  Lattice_Points_2d
clear all, close all;
%create an FIR gaussian filter with zero insertion to show several periods
%of the frequency response.
x = -3:.25:3;
y = (-3:.25:3)*sqrt(3)/2;
[X, Y] = meshgrid(x,y);
r0 = .65;
h = exp(-(X.^2 + Y.^2)/(2*r0^2));
h(X.^2+Y.^2>9)=0;
mask = zeros(25,25);
mask(1:8:25,3:4:25)=1; mask(5:8:25,1:4:25)=1;
h = h.*mask;
h = h /sum(sum(h));
[H,u,v] = freqz2(h,64, 64,[.25 sqrt(3)/8]);
[U,V] = meshgrid(u,v);
figure;
lev=0:.1:1.;             % contours will be from 0 to 1 in steps of 0.2
[C,h1]=contour(U,V,abs(H),lev); % generate the contour plot, including values to label contours
axis equal
lev1 = 0:.2:1.;
clabel(C,h1,lev1,'FontSize',6);           %label the contours  
xlabel('\it u \rm (c/X)'), ylabel('\it v \rm (c/X)');
set(gca,'ydir','reverse'); 
colormap([0 0 0]);     % use black only
hold on;
%plot points of the reciprocal lattice
VLs = [1 0; -1/sqrt(3) 2/sqrt(3)];
UL = [-2 -4/sqrt(3)]';
LR = [2 4/sqrt(3)]';
%get the points of the lattices in the desired region 
xoutLs = Lattice_Points_2d(VLs,UL,LR);
%plot the points
points = plot(xoutLs(:,1),xoutLs(:,2),'k.');
points.MarkerSize= 9;
%draw axes of symmetry
line([-2 2],[0 0],'Color','b','LineWidth',.55);
line([-4/3 4/3],[4/sqrt(3) -4/sqrt(3)],'Color','b');
line([4/3 -4/3],[4/sqrt(3) -4/sqrt(3)],'Color','b');
line([0 0],[-4/sqrt(3) 4/sqrt(3)],'Color','r','LineWidth',.55);
line([-2 2],[2/sqrt(3) -2/sqrt(3)],'Color','r');
line([-2 2],[-2/sqrt(3) 2/sqrt(3)],'Color','r');
%draw unit cell of reciprocal lattice
hex = polyshape([2/3 1/3 -1/3 -2/3 -1/3 1/3],[0 1/sqrt(3) 1/sqrt(3) 0 -1/sqrt(3) -1/sqrt(3)]);
plot(hex,'FaceColor','white','FaceAlpha',0);
set(gcf,'Color',[1 1 1]);