www.pudn.com > RadarSignals.rar > cross_ambfn2.m
% 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