% 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 datau_basic=input(' Signal elements (row complex vector, each element last tb sec) = ? ');m_basic=length(u_basic);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) = ? ');endF=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);   phas=uamp*0;   phas=angle(u_basic);   if fcode==1      phas=phas+2*pi*cumsum(f_basic);   end   uexp=exp(j*phas);   u=uamp.*uexp;else                               % i.e., several samples within a bit   dt=1/r;	                       % interval between samples   ud=diag(u_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;   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;   endendt=[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;% plot the signal parametersfigure(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 ')% 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 dtdtau=ceil(T*m)*dt/N;tau=round([0:1:N]*dtau/dt)*dt;% calculate K+1 equally spaced grid points of Doppler axis with df spacingf=[0:1:K]*df; % 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);% 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 sidesu_padded=[zeros(1,ceil(T*m)),u,zeros(1,ceil(T*m))];% define column indexing and row indexing vectorscidx=[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 uindex = cidx(ones(N+1,1),:) + ridx(:,ones(1,m+ceil(T*m)));% calculate matrixmat2 = 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);% 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 1a_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 Dopplera=[flipud(conj(a_pos(1:K+1,:))) fliplr(a_pos(K+2:2*K+2,:))];% define new delay and Doppler vectors delay=[-fliplr(tau) tau];freq=f(K+2:2*K+2)*ceil(max(t));% exclude the zero Delay that was taken twicedelay=[delay(1:N) delay((N+2):2*(N+1))];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 mapcm=zeros(64,3);  		cm(:,3)=ones(64,1); 	   figure(2), clf, hold offmesh(delay, [0 freq], [zeros(1,amt);a])hold onsurface(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% % Prof. Nadav Levanon% % Department of Electrical Engineering - Systems% The Iby and Aladar Fleischman Faculty of Engineering% Tel Aviv University% % This site contains two basic MATLAB programs for plotting ambiguity functions. % The more basic program ambfn1.m prompts for the signal (a complex row vector), % for frequency coding if it exists, and for various parameters that define the % desired ambiguity function plot. The program produces an additional figure, which % includes three subplots: Amplitude, phase and frequency of the signal. % % The more elaborate program ambfn7.m creates a graphic user interface (GUI). % The GUI allows choosing one of many preset signals, or defining a new signal % through its amplitude, phase and frequency vectors. % % The GUI's four push buttons, call one or more of the following programs: % % cal_and_plot_amb_fn7.m % calplotsig7.m % cal_and_plot_acf_and_spec7.m % cal_and_plot_pamb7.m % % These four additional m-files calculate and plot, respectively: % (1) the ambiguity function, % (2) the signal, % (3) its autocorrelation function, periodic autocorrelation and frequency spectrum and % (4) the periodic ambiguity function. % The four additional m-files must be stored in the same directory as ambfn7.m. % % HINTS ON USING AMBFN7.M:% % The signal is defined by the vectors in the three slots in the GUI: Amplitude, Phase/pi and frequency*tb.% When you choose a preset signal the slots are filled automatically.% % You can define a signal externally (in the MATLAB Command window) by entering three vectors respectively: u_amp, u_phase and f_basic. (They must all be row vectors of the same length.)% % The u_amp vector must be defined. One or two of the other vectors can be avoided, but then you need to toggle off the dot next to the corresponding slot in the GUI (by clicking on it).% % There are 5 parameters in the GUI:  r   F*Mtb   T    N    K % % Of those, the only parameter that affects the signal and not only the plots is "r". Since the signal is defined by a vector, with a well defined length (number of elements , referred to as M), it is often necessary to increase the number of samples (repeats) during each of these elements ("bits"), in order to meet the Nyquist criterion. This is the function of  r . % % Suggested values for r:% % In a Costas signal of M elements the signal bandwidth is approximately M/tb. Therefore, the sampling interval should be % % ts < tb/(2M)% % Hence % % tb/ts = r > 2M % % % In a phase-coded signal, the main spectral lobe ends at f = 1/tb. However, the spectral sidelobes extend much further at a rate of approximately 6 dB per octave. A typical spectral skirt crosses the -30 dB level at f = 10/tb.  Hence, choosing r = 2 is the minimum setting, but using r > 10 is recommended. % % Any time you change r  you need to click on the "Cal&Plot Sig." button to recalculate the signal. Only then you can click on any of the other buttons to get the different plots.% % When you change any of the other GUI parameters (F*Mtb   T   N    K ), you can usually click on the different plot buttons without recalculating the signal.% % The GUI parameters (F*Mtb  T  N   K ), are associated only with the plots.% % F*Mtb - defines the extent of the plotted Doppler axis in normalized Doppler (Doppler multiplied by the entire duration of the signal (Mtb)). % F*Mtb also defines the extent of the frequency axis in the spectrum plot (the lower subplot created when you click on the "ACF. & SPEC Plot.").  While the maximum value in the F*Mtb ruler is 60, you can type in a higher value .% % T - is the extent of the (positive) delay axis in units of the entire duration of the signal, so it is actually T/Mtb.% For example, if you choose the signal "pulse train, 6 pulses" and choose T=1, the ambiguity plot will include all the 5 recurrent lobes. If you choose T=1/6 you will get exactly one repetition period.% % N - is the number of grid points on the (positive) delay axis of the plot.% % K - is the number of grid points on the Doppler axis of the plot.% % % % Regarding the periodic ambiguity function (PAF): % % If you want to choose one of the preset signals, but wish to plot the PAF when the processor is matched to several periods of the signal (e.g., 4 periods), then you can do the following:% % In the GUI  choose  the desired signal and click on the "Cal&Plot Sig." button. % % Next, in the MATLAB command window enter three commands as follows:% % u_amp=[u_amp u_amp u_amp u_amp];% u_phase=[u_phase u_phase u_phase u_phase];% f_basic=[f_basic f_basic f_basic f_basic];% % Next click again on the "Cal&Plot Sig." button.% % In this way you have modified the signal to 4 repeats of the original signal. If you now plot the periodic ambiguity function, you will see the difference. (The volume is concentrated near Doppler values of n/period.) % % If you wish to display exactly 2 periods of the PAF, you can change T from 1 to T=0.25 .% 1．帶狀（對角）稀疏矩陣% 函數 spdiags% 格式 [B,d] = spdiags(A) %從矩陣A中提取所有非零對角元素，這些元素保存在矩陣B中，向量d表示非零元素的對角線位置。% B = spdiags(A,d) %從A中提取由d指定的對角線元素，并存放在B中。% A = spdiags(B,d,A) %用B中的列替換A中由d指定的對角線元素，輸出稀疏矩陣。% A = spdiags(B,d,m,n) %產生一個m×n稀疏矩陣A，其元素是B中的列元素放% 在由d指定的對角線位置上。% 例1-100% >>A = [11 0 13 0% 0 22 0 24% 0 0 33 0% 41 0 0 44% 0 52 0 0% 0 0 63 0% 0 0 0 74];% >>[B,d] = spdiags(A) % B =% 41 11 0% 52 22 0% 63 33 13% 74 44 24% d =% -3 %表示B的第1列元素在A中主對角線下方第3條對角線上% 0 %表示B的第2列在A的主對角線上% 2 %表示B的第3列在A的主對角線上方第2條對角線上% 例1-101% >> B=[1 2 3 4% 5 6 7 8% 9 10 11 12% 13 14 15 16]; % >> d=[-2 0 1 3];% >> A=spdiags(B,d,4,4);% >> full(A)% ans =% 2 7 0 16% 0 6 11 0% 1 0 10 15% 0 5 0 14