% wavefft2lband
% usage: wavefft2lband(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 wavefft2lband(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',[n,l],[0,k],'-*r',[0,l],[m,k],'-*r', [n,l],[m,k],'-*r' );
   axis ij
   axis([0 n 0 m]);
   axis('square')
   title('frequency')
   hold off
   subplot(2,2,2)
   hold on
   mk = min(k,m-k);
   ml = min(l,n-l)
   plot(u,sqrt(mk^2+ml^2),'*r');
   axis xy   
   axis([0, K, 0, 1.1*K]);
   title('frequency');
   subplot(2,2,3)
   imshow(real(wave),[-1,1]);
   axis('square')
   title('real part');
   subplot(2,2,4)
   imshow(imag(wave),[-1,1]);
   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;

