% EE 438 Digital Signal Processing with Applications
% 2-D signals and systems
% Sampling a 2-D cosine wave
% sine_sampling.m
%---------------------------------------------------------------
% Initialize
clear all;
close all;
%---------------------------------------------------------------
% Parameters
u_0 = 1/8	% spatial frequencies
v_0 = 1/32
C = 8;  	% sampling interval in x and y
M = 256;  	% size of image, must be even
%---------------------------------------------------------------
%Initialize colormap
mcomp = 0 : 1.0 ./ 255 : 1.0; % output colors lie in [0,1]
map = [mcomp' mcomp' mcomp'];

%gamma = 1.0;
%gamma_inv = 1/gamma;
index = (0:1:M-1);
x = ones(M,1)*index;
y = x(:,end:-1:1)';
%degrees = 10;
%theta = pi*degrees/180;
%P = 10;
%rho = 1/P;
%u_0 = rho*cos(theta)
%v_0 = rho*sin(theta)
g = 0.5 + 0.5*cos(2*pi*(u_0*x + v_0*y));
figure(1);
subplot(2,3,1)
%image(g.^gamma_inv*255);
image(g*255);
colormap(map);
axis('off');
axis('image');
title('original cosine wave')

g = fftshift(g);
G = fft2(g);
G = fftshift(G);
Gmag = abs(G);
Gmag = conv2(Gmag, ones(2,2), 'same');
Gmax = max(max(Gmag));
figure(1);
subplot(2,3,4)
image(Gmag*255/Gmax);
colormap(map);
axis('off');
axis('image');
title('original spectrum')
clear G Gmag;
%----------------------------------------
%Downsample, and then upsample by (C,D)
g_s = zeros(M,M);
g_s(1:C:M,1:C:M) = g(1:C:M,1:C:M);
g_s_disp = conv2(g_s, ones(2,2), 'same');
%g_s_disp = g_s;
figure(1);
subplot(2,3,2)
image(g_s_disp*255);
colormap(map);
axis('off');
axis('image');
title('sampled cosine')

g_s = fftshift(g_s);
G_s = fft2(g_s);
G_s = fftshift(G_s);
G_smag = abs(G_s);
G_smag = conv2(G_smag, ones(2,2), 'same');
G_smax = max(max(G_smag));


%-----------------------------------------
% Show baseband region
U = round(M/C);
u0 = round(M/2 - U/2);
v0 = round(M/2 - U/2);
u1 = u0+U-1;
v1 = v0;
u2 = u1;
v2 = v0+U-1;
u3 = u0;
v3 = v2;
for k = 0:U-2,	
	G_smag(u0+k+1,v0+1) = G_smax;
	G_smag(u1+1,v1+k+1) = G_smax;
	G_smag(u2-k+1,v2+1) = G_smax;
	G_smag(u3+1,v3-k+1) = G_smax;
end

figure(1);
subplot(2,3,5)
image(G_smag*255/G_smax);
colormap(map);
axis('off');
axis('image');
title('spectrum of sampled cosine')

clear g_s;
clear G_smag;
%-----------------------------------------
% Reconstruct
G_r = zeros(M,M);
k1=round(M/2);
k2=round(U/2);
G_r(k1-k2+1:k1+k2,k1-k2+1:k1+k2) = ...
C.*C.*G_s(k1-k2+1:k1+k2,k1-k2+1:k1+k2);
G_rmag = abs(G_r);
G_rmag = conv2(G_rmag, ones(2,2), 'same');
G_rmax = max(max(G_rmag));

figure(1);
subplot(2,3,6)
image(G_rmag*255/G_rmax);
colormap(map);
axis('off');
axis('image');
title('spectrum of reconstruction')

clear G_s;
clear G_rmag;

G_r = fftshift(G_r);
g_r = ifft2(G_r);
g_r = fftshift(g_r);

figure(1);
subplot(2,3,3)
%image(g_r.^gamma_inv*255);
image(abs(g_r)*255);		% g_r will have a tiny imaginary component
colormap(map);
axis('off');
axis('image');
title('reconstruction')