clear all; close all;
addpath('../lib');

source = 'hw8prob5data';
a = imbytefile(source,512); close all;

subplot(121);
imagesc(a);
colormap gray;
axis equal;
xlim([1 512]); ylim([1 512]);
title('Original');

A = fft2(a);

p  = 0.0012; % Pixel pitch (rad)
Fm = 1/2*1/p; % Max frequency represented (inverse rad)
[X Y] = meshgrid(linspace(-Fm,Fm,512));
W = 62.5;
C = 5*chat(sqrt(X.^2+Y.^2)/W);
C = C + 1*chat(sqrt((X-3*W).^2+Y.^2)/W);
C = C + 1*chat(sqrt((X+3*W).^2+Y.^2)/W);
C = C + 1*chat(sqrt((X-5*W).^2+Y.^2)/W);
C = C + 1*chat(sqrt((X+5*W).^2+Y.^2)/W);
C = C + 1*chat(sqrt(X.^2+(Y-4*W).^2)/W);
C = C + 1*chat(sqrt(X.^2+(Y+4*W).^2)/W);
% figure; surf(X,Y,C); shading interp;
% xlabel('u');
% ylabel('v');
% zlabel('Spectral coverage (a.u.)');
% title('Problem 4');
% return;

C = fftshift(C);
B = C.*A;
b = abs(ifft2(B));

subplot(122);
imagesc(b);
colormap gray;
axis equal;
xlim([1 512]); ylim([1 512]);
title('P4 aperture (T-shaped)');