% projective prolet

disp('Takes about an hour on a Decstation 5000/200')
disp('could speed up by not interpolating for the zeros')

begtime = fix(clock)

NW = 12;

M = 100;    % size of final output image
N = 100;

if ~exist('V')
  [V d] = slepian(200,NW);
end%if

SS_projf = zeros(M,N);   % size of image to be 100 by 100
t = linspace(-1,1,600).';  % this will really be f
tproj = project(t+1.34,pi/4,1,2); % 1.34 chosen to keep within [-1,1] and get 
                                  % projlet close to center
% actually this is a proj in freq domain, just use ``t'' for simplicity
for nw     = 1:2*NW   % only bother with the big fish
% for diagnostics for nw     = 23:24   % only bother with the big fish
  s_part   = V(:,nw); % slepian is one column of V
  s        = [zeros(s_part);s_part;zeros(s_part)];
  S        = fftshift(fft(s)); % get into freq domain
  disp(sprintf('%g  <---------------------------- PROJECTING SLEPIAN %g',nw,nw))
  S_projt = interp2(t(:),S(:),tproj(:));  % interpolate new projected pts in F
%  s_projf = ifft(fftshift(S_projt));
  s_projf = ifft(S_projt);
  SS_temp  = tf(s_projf,M,N);
  SS_projf = SS_projf + SS_temp.*conj(SS_temp);  % must use SQR magnitude
end%for

endtime = fix(clock)
disp('took this much time to generate the one projected prolet:')
endtime-begtime

disp('now about to display the projected prolet')
dispmatb(SS_projf)

