% EE 438 Digital Signal Processing with Applications - Spring 1997
% 2-D signals and systems
% 2-D grating pattern and spectrum
% grating.m
%---------------------------------------------------------------
% Initialize
clear all;
close all;
%---------------------------------------------------------------
% Parameters
P = 8;		% period of bar grating 
duty = 0.5;	% duty cycle of bar grating
W = 32; 	% width of window limiting bar grating
M = 256;  	% size of image
%---------------------------------------------------------------
%Initialize colormap
mcomp = 0 : 1.0 ./ 255 : 1.0; % output colors lie in [0,1]
map = [mcomp' mcomp' mcomp'];
index = (0:1:M-1);
x = ones(M,1)*index;
y = x(:,end:-1:1)';
degrees = 0;	% angle of bar grating
theta = pi*degrees/180;
g = abs(rem(cos(theta)*x + sin(theta)*y, P)) < duty*P;
figure(1);
subplot(121)
image(g*255);
colormap(map);
axis('off');
axis('image');
title('g(x,y)');

g_shift = fftshift(g);
G = fft2(g_shift);
G = fftshift(G);
Gmag = abs(G);
Gmax = max(max(Gmag));
Gmag_disp = conv2(Gmag, ones(2,2), 'same');
%figure(2);
subplot(122)
image(Gmag_disp*255/Gmax);
colormap(map);
axis('off');
axis('image');
title('|G(u,v)|');

%Limit vertical to width W
degrees = 90;		% angle of window
alpha = pi*degrees/180;
g1 = g.*(abs(cos(alpha)*(x-M/2+1) + sin(alpha)*(y-M/2+1)) <= W/2);
figure(3);
subplot(121)
image(g1*255);
colormap(map);
axis('off');
axis('image');
title('g_1(x,y)');

g_shift = fftshift(g1);
G = fft2(g_shift);
G = fftshift(G);
Gmag = abs(G);
Gmax = max(max(Gmag));
%figure(4);
subplot(122)
image(Gmag*255/Gmax);
colormap(map);
axis('off');
axis('image');
title('|G_1(u,v)|');


%Limit horizontal to width W
degrees = 0;		% angle of window
alpha = pi*degrees/180;
g2 = g.*(abs(cos(alpha)*(x-M/2+1) + sin(alpha)*(y-M/2+1)) <= W/2);
figure(5);
subplot(121)
image(g2*255);
colormap(map);
axis('off');
axis('image');
title('g_2(x,y)');

g_shift = fftshift(g2);
G = fft2(g_shift);
G = fftshift(G);
Gmag = abs(G);
Gmax = max(max(Gmag));
%figure(6);
subplot(122)
image(Gmag*255/Gmax);
colormap(map);
axis('off');
axis('image');
title('|G_2(u,v)|');


% generate grating pattern 
g3 = g.*g';
figure(7);
subplot(121)
image(g3*255);
colormap(map);
axis('off');
axis('image');
title('g_3(x,y)');

g_shift = fftshift(g3);
G = fft2(g_shift);
G = fftshift(G);
Gmag = abs(G);
Gmag_disp = conv2(Gmag, ones(2,2), 'same');
Gmax = max(max(Gmag));
%figure(8);
subplot(122)
image(Gmag_disp*255/Gmax);
colormap(map);
axis('off');
axis('image');
title('|G_3(u,v)|');



% truncated grating pattern 
g4 = g2.*g2';
figure(9);
subplot(121)
image(g4*255);
colormap(map);
axis('off');
axis('image');
title('g_4(x,y)');

g_shift = fftshift(g4);
G = fft2(g_shift);
G = fftshift(G);
Gmag = abs(G);
Gmax = max(max(Gmag));
%figure(10);
subplot(122)
image(Gmag*255/Gmax);
colormap(map);
axis('off');
axis('image');
title('|G_4(u,v)|');