% pt2mov.m
% from pt2d_1.m with target motion (9/02)
% HAS NEGATIVE DOPPLER & FINDS MAX VELOCITY BIN
% processing a pulse train for doppler and range information
% multiple targets (stationary or moving)
 clear
 clf
% radar parameters
 N=64;        % number of pulses in train
 tau=1e-6;    % pulse width
 fp=10e3;     % PRF
 Tp=1/fp;     % interpulse period  
 T=N*Tp;      % total observation time
 dt=tau;      % range gate time
 df=1/T;      % frequency increment 
 fc=9e9;      % carrier frequency
 c=3e8;
 wave=c/fc;
 dr=c*tau/2;  % range bin size (range resolution)
 disp(['range resolution (m): ',num2str(dr)])
 M=floor(Tp/tau-1);  % number of range bins/gates  
% target velocities
 v(1,1:2)=-[30 30];  
 v(2,1:2)=[10 -5];
 v(3,1:2)=[10 -45];
 [ntarg,dum]=size(v);  %  get number of targets from velocity vector size
 for i=1:ntarg
    disp(['target number ',num2str(i),' velocity vector, [vx,vy] (m/s): ',num2str(v(i,1:2))])
 end
 dmx=5000;
% assign target start locations at time 0 in meters
% assume farthest range is at t=0
 xo(1)=6000; yo(1)=6000;
 xo(2)=3000; yo(2)=0;
 xo(3)=-4500; yo(3)=7500;
 pmx=max(max(xo),max(yo));  Ro=sqrt(xo.^2+yo.^2); Rmx=max(Ro)/1000; 
 for i=1:ntarg
 disp(['target number ',num2str(i),' start location, [xo,yo] (m): ',num2str([xo(i),yo(i)])])
 end
 % relative scattering amplitudes of targets
 Am=[1 1 1];
 Ph=deg2rad([0 0 0]);
 figure(1)
 clf
% time loop (for target motion) *****************************************************
t1=0;    % start time
t2=350;   % stop time
dt=5;   % time step
Nt=floor((t2-t1)/dt)+1;   % number of time steps
for kt=1:Nt
% initialize I and Q to zero for all range and velocity bins
    I=zeros(N,M); Q=zeros(N,M);
    t(kt)=t1+(kt-1)*dt;
    % loop through all targets
    I=zeros(N,M); Q=zeros(N,M);
    for i=1:ntarg
    x(i,kt)=xo(i)+v(i,1)*t(kt);
    y(i,kt)=yo(i)+v(i,2)*t(kt);
    R(i,1:2)=[x(i,kt) y(i,kt)];   % range to target
    magR=sqrt(dot(R(i,1:2),R(i,1:2))); RNG(kt)=magR;
    Rhat=R(i,1:2)/(magR+1);
    mbin(i)=floor(magR/dr)+1; MBIN(i,kt)=mbin(i);
    vr(i)=dot(Rhat,v(i,1:2));
    fd(i)=2*vr(i)/wave;
    if abs(fd(i)) > 5000
        disp(['target ',num2str(i),' amb vel at time step ',num2str(kt)])
    end
 
% pulse loop
  for n=1:N
     tn=(n-1)*Tp+t(kt);
     I(n,mbin(i))=I(n,mbin(i))+Am(i)*cos(2*pi*fd(i)*tn+Ph(i));
     Q(n,mbin(i))=Q(n,mbin(i))+Am(i)*sin(2*pi*fd(i)*tn+Ph(i));
  end
 end   % end of target loop
% add noise
 No=randn(N,M)*.1;
 I=I+No; Q=Q+No;

%------------- I and Q signals generated --------------------------------
% begin processing -- look in each range bin for a signal
 thresh=.2;
 it=0;
 Z=zeros(M,N);
 for m=1:M
   sum=0; 
   rng(m)=m*dr;
   A1=I(:,m)+j*Q(:,m);
   B1=fft(A1); As1=fftshift(B1);
% don't discard negative frequencies
% Ap1=As1';
    Ap1=As1; 
    if m==1, vel=-([0:N-1]-(N/2-1))*df*wave/2; end
    Z(m,:)=Ap1(:).';
 end 
   dop=vel*2/wave; rgate=[1:M]; 
   figure(1)
   subplot(121)
   for i=1:ntarg
   xymag(i,1:kt)=sqrt(x(i,:).^2+y(i,:).^2)/1000;
   xyang(i,1:kt)=angle(x(i,:)+j*y(i,:));
   polar(0,Rmx),axis square, axis equal, hold on
   polar(xyang(i,1:kt),xymag(i,1:kt),'k-')
   if i==1, ss='o'; end
   if i==2, ss='s'; end
   if i==3, ss='*'; end
   if i==4, ss='^'; end
   if i==5, ss='x'; end
   polar(xyang(i,kt),xymag(i,kt),ss)
   text(Rmx-.5*Rmx,Rmx+.3*Rmx,'km')
   end
   hold off
   text(-Rmx/2,Rmx+.55*Rmx,'Target Paths')
   subplot(122)
%   mesh(vel,rng/1000,abs(Z)),grid on,view(-20,20),axis square
%   xlabel('velocity (m/s)'), ylabel('range (km)'), zlabel('Magnitude')
   mesh(dop/1000,rgate*dr/1000,abs(Z)), view(-50,30), 
   xlabel('Doppler (kHz)'), ylabel('Range (km)'), zlabel('Magnitude')
   axis equal, axis square
   title('Doppler/Range Space')
   axis([-dmx/1000,dmx/1000,0,M*dr/1000,0,N])
   pause(0.01)
   K(kt)=getframe(gcf);
end  % end of time loop ***********************************************************
figure(2),axis off, movie(K,5);