% cross_ambfn2.m - amodification of "ambfn1.m" for plotting cross ambiguity
%                   between two signals OF THE SAME LENGTH
%                   Designed to allow mismatch caused by FFT Doppler processing a pulse train.
%                   The two signals differ only in phases (if there is
%                   frequency mode it is the same in both signals.
% ambfn1.m - plots ambiguity function of a signal u_basic (row vector) 
%
% The m-file returns a plot of quadrants 1 and 2 of the ambiguity function of a signal 
% The ambiguity function is defined as:
%									  
% a(t,f) = abs ( sumi( u(k)*u'(i-t)*exp(j*2*pi*f*i) ) )
%					 
% The user is prompted for the signal data:
% u_basic is a row complex vector representing amplitude and phase
% f_basic is a corresponding frequency coding sequence
%
% The duration of each element is tb (total duration of the signal is tb*(m_basic-1))
%
% F is the maximal Dopler shift
% T is the maximal Delay
% K is the number of positive Doppler shifts (grid points)
% N is the number of delay shifts on each side (for a total of 2N+1 points)
% The code allows r samples within each bit
%
% Written by Eli Mozeson and Nadav Levanon, Dept. of EE-Systems, Tel Aviv University

% clear all

% prompt for signal data
u_basic=input(' Signal elements (row complex vector, each element last tb sec) = ? ');
m_basic=length(u_basic);
v_basic=input(' 2nd Signal elements (row complex vector, each element last tb sec) = ? ');

fcode=input(' Allow frequency coding (yes=1, no=0) = ? ');
if fcode==1
    f_basic=input(' Frequency coding in units of 1/tb (row vector of same length) = ? ');
end

F=input(' Maximal Doppler shift for ambiguity plot [in units of 1/Mtb] (e.g., 1)= ? ');
K=input(' Number of Doppler grid points for calculation (e.g., 100) = ? ');
df=F/K/m_basic;

T=input(' Maximal Delay for ambiguity plot [in units of Mtb] (e.g., 1)= ? ');
N=input(' Number of delay grid points on each side (e.g. 100) = ? ');

sr=input(' Over sampling ratio (>=1) (e.g. 10)= ? ');
r=ceil(sr*(N+1)/T/m_basic);

if r==1
   dt=1;
   m=m_basic;
   uamp=abs(u_basic);
                    vamp=abs(v_basic);
   phas=uamp*0;
   phas=angle(u_basic);
                   phasv=angle(v_basic);  
   if fcode==1
      phas=phas+2*pi*cumsum(f_basic);
                   phasv=phas+2*pi*cumsum(f_basic);

   end
   uexp=exp(j*phas);
   u=uamp.*uexp;
                  vexp=exp(j*phasv);
                  v=vamp.*vexp;

   
else                               % i.e., several samples within a bit
   dt=1/r;	                       % interval between samples
   ud=diag(u_basic);
                       vd=diag(v_basic);
   ao=ones(r,m_basic);
   m=m_basic*r;
   u_basic=reshape(ao*ud,1,m);    % u_basic with each element repeated r times
   uamp=abs(u_basic);
   phas=angle(u_basic);
   u=u_basic;
   
                        v_basic=reshape(ao*vd,1,m);    % v_basic with each element repeated r times
                        vamp=abs(v_basic);
                        phasv=angle(v_basic);
                        v=v_basic;
   if fcode==1
       ff=diag(f_basic);
       phas=2*pi*dt*cumsum(reshape(ao*ff,1,m))+phas;
       uexp=exp(j*phas);
       u=uamp.*uexp;
       
       phasv=2*pi*dt*cumsum(reshape(ao*ff,1,m))+phasv;
       vexp=exp(j*phasv);
       v=vamp.*vexp;
   end
end

t=[0:r*m_basic-1]/r;
tscale1=[0 0:r*m_basic-1 r*m_basic-1]/r;
dphas=[NaN diff(phas)]*r/2/pi;
                   dphasv=[NaN diff(phasv)]*r/2/pi;

% plot the signal parameters
figure(1), clf, hold off 
subplot(3,1,1)
plot(tscale1,[0 abs(uamp) 0],'linewidth',1.5)
ylabel(' Amplitude ')
axis([-inf inf 0 1.2*max(abs(uamp))])

subplot(3,1,2)
plot(t, phas,'linewidth',1.5)
axis([-inf inf -inf inf])
ylabel(' Phase [rad] ')

subplot(3,1,3)
plot(t,dphas*ceil(max(t)),'linewidth',1.5)
axis([-inf inf -inf inf])
xlabel(' \itt / t_b ')
ylabel(' \itf * Mt_b ')


               % plot the 2nd signal parameters
                 figure(2), clf, hold off 
                 subplot(3,1,1)
                 plot(tscale1,[0 abs(vamp) 0],'linewidth',1.5)
                 ylabel(' Amplitude ')
                 axis([-inf inf 0 1.2*max(abs(vamp))])

                 subplot(3,1,2)
                 plot(t, phasv,'linewidth',1.5)
                 axis([-inf inf -inf inf])
                 ylabel(' Phase [rad] ')

                 subplot(3,1,3)
                 plot(t,dphasv*ceil(max(t)),'linewidth',1.5)
                 axis([-inf inf -inf inf])
                 xlabel(' \itt / t_b ')
                 ylabel(' \itf * Mt_b ')


% calculate a delay vector with N+1 points that spans from zero delay to ceil(T*t(m))
% notice that the delay vector does not have to be equally spaced but must have all
% entries as integer multiples of dt

dtau=ceil(T*m)*dt/N;
% tau=round([0:1:N]*dtau/dt)*dt;
tau=round([0:1:2*N]*dtau/dt)*dt;

% calculate K+1 equally spaced grid points of Doppler axis with df spacing
f=[0:1:K]*df;
ff=f;

% duplicate Doppler axis to show also negative Dopplers (0 Doppler is calculated twice)
f=[-fliplr(f) f];

% calculate ambiguity function using sparse matrix manipulations (no loops)

% define a sparse matrix based on the signal samples u1 u2 u3 ... um
% with size m+ceil(T*m) by m (notice that u' is the conjugate transpose of u)
% where the top part is diagonal (u*) on the diagonal and the bottom part is a zero matrix
%
%			[u1*  0   0  0 ...  0  ] 
%			[ 0  u2*  0  0 ...  0  ]
%			[ 0   0  u3* 0 ...  0  ]	m rows
%			[ .				 .	  .  ]
%			[ .				 .	  .  ]
%			[ .   0   0	 . ...  um*]
%			[ 0					  0  ]		
%			[ .					  .  ]   N rows
%			[ 0   0   0  0 ...  0  ]
%
%  mat1=spdiags(u',0,m+ceil(T*m),m);  <====== replaced by the 2nd signal
mat1=spdiags(v',0,m+ceil(T*m),m);


% define a convolution sparse matrix based on the signal samples u1 u2 u3 ... um
% where each row is a time(index) shifted versions of u.
% each row is shifted tau/dt places from the first row 
% the minimal shift (first row) is zero
% the maximal shift (last row) is ceil(T*m) places
% the total number of rows is N+1
% number of columns is m+ceil(T*m)

% for example, when tau/dt=[0 2 3 5 6] and N=4
%
%			[u1 u2 u3 u4  ...               ... um  0  0  0  0  0  0]
%			[ 0  0 u1 u2 u3 u4  ...               ... um  0  0  0  0]
%			[ 0  0  0 u1 u2 u3 u4  ...               ... um  0  0  0]
% 			[ 0  0  0  0  0 u1 u2 u3 u4  ...               ... um  0]
%			[ 0  0  0  0  0  0 u1 u2 u3 u4  ...               ... um] 
 
% define a row vector with ceil(T*m)+m+ceil(T*m) places by padding u with zeros on both sides
% u_padded=[zeros(1,ceil(T*m)),u,zeros(1,ceil(T*m))];
u_padded=[zeros(1,ceil(T*m)),u,zeros(1,2*ceil(T*m))];

% define column indexing and row indexing vectors
cidx=[1:m+ceil(T*m)];
ridx=round(tau/dt)';

% define indexing matrix with Nused+1 rows and m+ceil(T*m) columns 
% where each element is the index of the correct place in the padded version of u
% index = cidx(ones(N+1,1),:) + ridx(:,ones(1,m+ceil(T*m)));
index = cidx(ones(2*N+1,1),:) + ridx(:,ones(1,m+ceil(T*m)));

[mmm,nnn]=size(index);
% calculate matrix
mat2 = sparse(u_padded(index)); 


% calculate the ambiguity matrix for positive delays given by 
%
%	[u1 u2 u3 u4  ...               ... um  0  0  0  0  0  0] [u1*  0   0  0 ...  0  ]
%	[ 0  0 u1 u2 u3 u4  ...               ... um  0  0  0  0] [ 0  u2*  0  0 ...  0  ]
%	[ 0  0  0 u1 u2 u3 u4  ...               ... um  0  0  0]*[ 0   0  u3* 0 ...  0  ]
% 	[ 0  0  0  0  0 u1 u2 u3 u4  ...               ... um  0] [ .			 .	  .  ]
%	[ 0  0  0  0  0  0 u1 u2 u3 u4  ...               ... um] [ .			 .	  .  ]
%                                                             [ .   0   0  . ...  um*]
%       												      [ 0		   	      0  ]		
%													          [ .		    	  .  ]  
%			                                                  [ 0   0   0  0 ...  0  ]
%
% where there are m columns and N+1 rows and each element gives an element 
% of multiplication between u and a time shifted version of u*. each row gives
% a different time shift of u* and each column gives a different entry in u.
%
uu_pos=mat2*mat1;

% clear mat2 mat1

% calculate exponent matrix for full calculation of ambiguity function. the exponent
% matrix is 2*(K+1) rows by m columns where each row represents a possible Doppler and
% each column stands for a differnt place in u.
% e=exp(-j*2*pi*f'*t);
e=exp(-j*2*pi*ff'*t);

% calculate ambiguity function for positive delays by calculating the integral for each
% possible delay and Doppler over all entries in u.
% a_pos has 2*(K+1) rows (Doppler) and N+1 columns (Delay)
a_pos=abs(e*uu_pos');

% normalize ambiguity function to have a maximal value of 1
a_pos=a_pos/max(max(a_pos));

% use the symmetry properties of the ambiguity function to transform the negative Doppler
% positive delay part to negative delay, positive Doppler
% a=[flipud(conj(a_pos(1:K+1,:))) fliplr(a_pos(K+2:2*K+2,:))];
a=a_pos;

% define new delay and Doppler vectors 
delay0=[-fliplr(tau) tau];
% freq=f(K+2:2*K+2)*ceil(max(t));
% freq=f*ceil(max(t));
freq=ff*ceil(max(t));

% exclude the zero Delay that was taken twice
% delay=[delay0(1:N) delay0((N+2):2*(N+1))];
delay=[delay0(N+1:2*N+1) delay0(2*N+3:(3*N+2))];

% a=a(:,[1:N (N+2):2*(N+1)]);

% plot the ambiguity function and autocorrelation cut
[amf amt]=size(a);

% create an all blue color map
cm=zeros(64,3);  		
cm(:,3)=ones(64,1); 	
   
figure(3), clf, hold off
mesh(delay, [0 freq], [zeros(1,amt);a])

hold on
surface(delay, [0 0], [zeros(1,amt);a(1,:)])

colormap(cm)
view(-40,50)
axis([-inf inf -inf inf 0 1])
xlabel(' {\it\tau}/{\itt_b}','Fontsize',12);
ylabel(' {\it\nu}*{\itMt_b}','Fontsize',12);
zlabel(' |{\it\chi}({\it\tau},{\it\nu})| ','Fontsize',12);
hold off