% wavefft2cband
% usage: wavefft2cband(m,n,q,steps);
% show circular bands of dft wave forms
%
% m,n: number of rows and columns of the image matrix
% q: frequency 
% steps: number of directions 

function wavefft2cband(m,n,q,steps);

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

%define frequency domain parameters 
hm = floor(m/2);
hn = floor(n/2);

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

 
start = 0;
for u = 0:1:steps
   k = round(q*cos(2*pi*u/steps));
   l = round(q*sin(2*pi*u/steps));
   wave = exp(2*pi*i*k*tm')*exp(2*pi*i*l*tn);
   subplot(2,2,1)
   hold on
   plot([0,l],[0,k],'.-r');
   axis ij
   axis([-hn hn -hm hm]);
   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, steps, 0, 1.1*q]);
   title('frequency');
   hold off 
   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 ~start 
      disp('press a key');
      pause
      start =1;
   end
   disp(['k=',num2str(k),', l=',num2str(l),', q=',num2str(sqrt(k^2+l^2))]) 
   pause(0.01)
   clear wave;
 end;