% wavefft2lband
% usage: wavefftlband(m,n,K);
% show linear bands of wave forms
%
% m,n: number of rows and columns of the image matrix
% K: sum of the vertical and horizontal frequencies 

function wave2dlband(m,n,K);

%define sample points
tm = (0:1:(m-1))/m;
tn = (0:1:(n-1))/n;

% set up graphics
figure(1)
clf;
set(gcf,'DefaultLineEraseMode','none');
set(gcf,'NextPlot','add');


%set up initial frequencies
if K <= m 
   k = K; l = 0;
else
   k = m; l = K-m;
end

u=0;
while (k >= 0) & (l<n)
em = zeros(m,1); em(k+1)=1; en = zeros(1,n); en(l+1)=1; wave = fft2(em*en);
%   wave = exp(2*pi*i*k*tm')*exp(2*pi*i*l*tn);
   subplot(2,2,1)   
   plot([0,l],[0,k],'-*r');
   axis ij
   axis([0 n 0 m]);
   axis('square')
   title('frequency vector')
   hold off
   subplot(2,2,2)
   hold on
   plot(u,sqrt(k^2+l^2),'*r');
   axis xy   
   axis([0, K, 0, 1.1*K]);
   title('computed (aliased) frequency');
   subplot(2,2,3)
   imagesc(real(wave));
   colormap(gray);
   axis('square')
   title('real part');
   subplot(2,2,4)
   imagesc(real(wave));
   colormap(gray);
   axis('square')
   title('imaginary part')
   if u < 1 
      disp('press a key');
      pause
   end
   u=u+1;
   disp(['k=',num2str(k),', l=',num2str(l)]) 
   pause(0.01)
   k=k-1; l=l+1;
end;