tdbt.m

阵列信号处理的工具箱

function tdbt(testnr)

% [ ]Ok *** TEST ***
% of: DBT (A Matlab Toolbox for Digital Beamforming).
% Sign:     Datum:
%
% *****************************************************************************
%   *  DBT, A Matlab Toolbox for Radar Signal Processing  *
% (c) FOA 1994-2000. See the file dbtright.m for copyright notice.
%
%  A test program for the Toolbox.
%
%  Start        : 970314 Svante Bj鰎klund (svabj).
%  Latest change: $Date: 2000/10/16 15:22:16 $ $Author: svabj $.
%  $Revision: 1.52 $
% *****************************************************************************

global thSign
  % Defines the reference direction for the theta angle in DOA:s.
  % Defined in "defant".

if (nargin < 1)
  disp('Test of A Matlab Toolbox for Digital Beamforming.')
  %disp('Don''t forget to compile mex-files before the test.')
  disp('Test 11:  Test of PTMF with the AIMT application example.')
  disp('Test 12:  Test of the function dspastemat for the aimtEx antenna.')
  disp('Test 13:  Test of recursive antenna definitions.')
  disp('Test 14:  Test of the function dspastemat for recursive antenna definitions.')
  disp('Test 16:  Comparison between ULA and LA.')
  disp('Test 17:  Test of the antenna types "array" and "beamform".')
  disp('Test 20:  Test of spectral DOA estimation after modifications.')
  disp('Test 21: ')
  disp('Test 25:  Test of conventional beamforming.')
  disp('Test 26:  Test of tapering of antenna pattern.')
  disp('Test 40:  Test of parametric DOA-estimation after modifications.')
  disp('Test 50:  Test of "expAnt" antenna with different calibration.')
  disp('Test 51:  Test of "expAnt" antenna with different calibration of steering vec.')
  disp('Test 60:  Test of different antenna types and steering matrix calculations.')
  disp('Test 101: Test of pulse code modulations.')
  disp('Test 102: Test of exchange getm, getm2, getm3 and setm with direct indexing.')
  disp('Test 103: Test of function "sigsplitdim".')
  disp('Test 120: Test of wide band signal simulation.')
  disp('Test 130: Test of range steering matrix calculation.')
  disp('Test 131: Test of space-range steering matrix calculation.')
  disp('Test : ')
  testnr=input('Choose testno: ');
  while isempty(testnr)
    testnr=input('Choose testno: ');
  end
end%if


if isstr(testnr)
  testnr = str2num(testnr)
end%if

if (testnr==1)
  % Nothing.


elseif (testnr==11)
% *****************************************************************************
% [ ]  : Test 11
% Sign:     Datum:
%
%  Test of PTMF with the AIMT application example.
%
% *****************************************************************************
% Comments:

% [ ]
% Comments:
%  Target angles d2r([-5 0 5]).
%  Start values d2r([-5 0 5]).
disp('Test 11a')
T  = 200;         % Number of snapshots.
K  = 25;       % Number of channels.

theta = d2r([-5 0 5])'; % Target angles. The number of targets is
            %  given by the number of target angles.
phi   = zeros(size(theta));   % Target angles.
SNR   = [10 10 10]';    % Signal to noise ratio in dB at each
            %  antenna element!?!
alpha = d2r([20 210 33])';    % Start phases of the target signals.
dalpha   = d2r([0 -10 150])'; % A constant phase shift between snapshots.
            %  Means targets movements at constant velocity.
dist=Inf*ones(size(theta));   % Distances to the sources.
tgtModel       = 'const';
noiseModel        = 'nonoise';
lambda = 0.1;
ant = defant('aimtEx');
sig = compsim4(ant, lambda, T, tgtModel, [theta, phi, SNR, alpha, dalpha, dist, eye(size(theta,1))], noiseModel, eye(K));
sig = msigmat2(sig, T);
[edoa, sHat, Q, doaAll] = doapar1('ptmf',sig, d2r([-5 0 5]), Inf, [], [], 5*1e-3);
pdoapar(edoa)
pdoait(doaAll)
perrfun(Q)
% Expected result :
%  In the text window: about [-5 0 5]' degrees.
%  The plot which shows the convergence of the doa parameters shall
%  show a smooth convergence to about [-5 0 5]' degrees.
%  The error function may be rough but shall decrease.
% Real result :



pause

% [ ]
disp('Test 11b')
% Comments:
%  Target angles d2r([-5 0 5]).
%  Start values d2r([-1 0 1]).
T  = 200;         % Number of snapshots.
K  = 25;       % Number of channels.

theta = d2r([-5 0 5])'; % Target angles. The number of targets is
            %  given by the number of target angles.
phi   = zeros(size(theta));   % Target angles.
SNR   = [10 10 10]';    % Signal to noise ratio in dB at each
            %  antenna element!?!
alpha = d2r([20 210 33])';    % Start phases of the target signals.
dalpha   = d2r([0 -10 150])'; % A constant phase shift between snapshots.
            %  Means targets movements at constant velocity.
dist=Inf*ones(size(theta));   % Distances to the sources.
tgtModel       = 'const';
noiseModel        = 'nonoise';
lambda = 0.1;
ant = defant('aimtEx');
sig = compsim4(ant, lambda, T, tgtModel, [theta, phi, SNR, alpha, dalpha, dist, eye(size(theta,1))], noiseModel, eye(K));
sig = msigmat2(sig, T);
[edoa, sHat, Q, doaAll] = doapar1('ptmf',sig, d2r([-1 0 1]), Inf, [], [], 5*1e-2);
pdoapar(edoa)
pdoait(doaAll)
perrfun(Q)
% Expected result :
%  In the text window: about [-5 0 5]' degrees.
%  The plot which shows the convergence of the doa parameters shall show a  convergence to about [-5 0 5]' degrees.
%  The error function probably does not decrease so much.
% Real result :




% [ ]
disp('Test 11c')
% Comments:
%  Target angles d2r([-1 0 1]).
%  Start values d2r([-.5 0 .5]).
T  = 200;         % Number of snapshots.
K  = 25;       % Number of channels.

theta = d2r([-1 0 1])'; % Target angles. The number of targets is
            %  given by the number of target angles.
phi   = zeros(size(theta));   % Target angles.
SNR   = [10 10 10]';    % Signal to noise ratio in dB at each
            %  antenna element!?!
alpha = d2r([20 210 33])';    % Start phases of the target signals.
dalpha   = d2r([0 -10 150])'; % A constant phase shift between snapshots.
            %  Means targets movements at constant velocity.
dist=Inf*ones(size(theta));   % Distances to the sources.
tgtModel       = 'const';
noiseModel        = 'nonoise';
lambda = 0.1;
ant = defant('aimtEx');
sig = compsim4(ant, lambda, T, tgtModel, [theta, phi, SNR, alpha, dalpha, dist, eye(size(theta,1))], noiseModel, eye(K));
sig = msigmat2(sig, T);
[edoa, sHat, Q, doaAll] = doapar1('ptmf',sig, d2r([-.5 0 .5]), Inf, [], [], 5*1e-4);
pdoapar(edoa)
pdoait(doaAll)
perrfun(Q)
% Expected result :
%  In the text window: about [-1 0 1]' degrees.
%  The plot which shows the convergence of the doa parameters shall show
%   a smooth convergence to about [-1 0 1]' degrees.
%  The error function shall decrease monotonically.
% Real result :




elseif (testnr==12)
% *****************************************************************************
% [ ]  : Test 12
% Sign:     Datum:
%
% Test of the function dspastemat for the aimtEx antenna
%
% *****************************************************************************
% Comments

% [ ]
% Comments
%  Test of the function dspastemat for the 'aimtEx' antenna.
disp('Test 12a')
n = 10;           % Index to SubArray to be calculated
m = 1;            % Index to target angle to be calculated

dElem = 0.05;
dSubArray = 0.1;
ant = defant('aimtEx');       % Application example antenna
theta = d2r([-20 -5 0 5 20]');      % Target angles. The no of targets is
phi = zeros(size(theta));     %  given by the no of target angles.
doas = ([theta';phi']);       % Only theta is considered as target angle
lambda = 0.1;           % Wavelength (3 GHz)
w = [0.5 1 0.5]';       % Tapering weights vector
%l1 = -j*2*pi/lambda*dElem;      % Temp variable
%l2 = -j*2*pi/lambda*dSubArray;     % Temp variable
l1 = thSign*j*2*pi/lambda*dElem;    % Temp variable
l2 = thSign*j*2*pi/lambda*dSubArray;      % Temp variable

delta_theta = 0.0000001;
th1 = theta(m) + delta_theta;
th2 = theta(m) - delta_theta;
part1 = exp((n-1)*l2*sin(th1))*(w(1)+w(2)*exp(l1*sin(th1))+w(3)*exp(2*l1*sin(th1)))*sqrt(cos(th1));
part2 = exp((n-1)*l2*sin(th2))*(w(1)+w(2)*exp(l1*sin(th2))+w(3)*exp(2*l1*sin(th2)))*sqrt(cos(th2));
Expected_Value = (part1 - part2) / (2*delta_theta);   % Numerical value
disp(['Expected value = ' num2str(Expected_Value)])

%for m=1:5
%  for n=1:25
%    th1 = theta(m) + delta_theta;
%    th2 = theta(m) - delta_theta;
%    part1 = exp((n-1)*l2*sin(th1))*(w(1)+w(2)*exp(l1*sin(th1))+w(3)*exp(2*l1*sin(th1)))*sqrt(cos(th1));
%    part2 = exp((n-1)*l2*sin(th2))*(w(1)+w(2)*exp(l1*sin(th2))+w(3)*exp(2*l1*sin(th2)))*sqrt(cos(th2));
%    AExp(n,m) = (part1 - part2) / (2*delta_theta);   % Numerical value
%  end
%end

A = dspastemat(ant,doas,lambda);
Real_Result = A(n,m);
disp(['Real value = ' num2str(Real_Result)])
%Real_Result = sum(sum(abs(A - AExp)))

% Expected result :
%   The values shall be about the same.
% Real result :



% [ ]
disp('Test 12b')
n = 22;           % Index to SubArray to be calculated
m = 4;            % Index to target angle to be calculated

delta_theta = 0.0000001;
th1 = theta(m) + delta_theta;
th2 = theta(m) - delta_theta;
part1 = exp((n-1)*l2*sin(th1))*(w(1)+w(2)*exp(l1*sin(th1))+w(3)*exp(2*l1*sin(th1)))*sqrt(cos(th1));
part2 = exp((n-1)*l2*sin(th2))*(w(1)+w(2)*exp(l1*sin(th2))+w(3)*exp(2*l1*sin(th2)))*sqrt(cos(th2));
Expected_Value = (part1 - part2) / (2*delta_theta);   % Numerical value
disp(['Expected value = ' num2str(Expected_Value)])

A = dspastemat(ant,doas,lambda);
Real_Result = A(n,m);
disp(['Real value = ' num2str(Real_Result)])

% Expected result :
%   The value shall be about the same
% Real result :




elseif (testnr==13)
% *****************************************************************************
% [ ]  : Test 13
% Sign:     Datum:
%
% Test of a recursive antenna definitions.
% ==Testing of generalized function spastemat.m (a recursive ULA, where the
% elements can be ULA:s or one-channel elements) to calculate the steering
% vector.
%
% *****************************************************************************
% Comments: Testing of generalized function to calculate the steering vector

% [ ]
disp('Test 13a')
clear all
global thSign
thSign = 1;
w      = [0.5 1 0.5]';      % Tapering weights vector
theta  = pi/4;              % Target directions
ant    = defant('aimtEx');  % Application example antenna
lambda   = 0.1;             % Wavelength of simulated radarsignal.
%Expected_Value = exp(-j*0.1*(0:24)'*2*pi/lambda*sin(theta)) * ...
%                 w' * exp(-j*0.05*(0:2)'*2*pi/lambda*sin(theta)) * ...
%                 sqrt(cos(theta));
Expected_Value = exp(thSign*j*0.1*(0:24)'*2*pi/lambda*sin(theta)) * ...
                 w' * exp(thSign*j*0.05*(0:2)'*2*pi/lambda*sin(theta)) * ...
                 sqrt(cos(theta));
Real_Result = spastemat(ant,theta,lambda);
disp(['Expected value < ' num2str(max(abs(Expected_Value))*1.5E-12)])
disp(['Real value = ' num2str(sum(abs(Expected_Value - Real_Result)))])
% Expected result : 0
%  The value should be about the same and diff should be small
% Real result :


% [ ]
disp('Test 13b')
clear all
global thSign
thSign = 1;
w1          = [0.5 1 0.5]';          % Tapering weights vector 1
w2          = [0.1 0.6 1 0.6 0.1]';  % Tapering weights vector 2
theta       = 2*pi/7;                % Target directions
lambda      = 0.0842;                % Wavelength of simulated radarsignal.
noElemMain  = 25;
noElemSub1  = 3;
noElemSub2  = 5;
distElemMain   = lambda;
distElemSub1   = lambda/3;
distElemSub2   = lambda/2;
elem     = defant('pattFuncElem', 'sqrt(cos(x(1,:)))');
subarray2 = defant('ULA', [noElemSub2, distElemSub2], elem, w2);
subarray1 = defant('ULA', [noElemSub1, distElemSub1], subarray2, w1);
ant      = defant('ULA', [noElemMain, distElemMain], subarray1, ones(noElemMain:1));
%Expected_Value = exp(-j*distElemMain*(0:noElemMain-1)'*2*pi/lambda*sin(theta)) * ...
%                 (w1' * exp(-j*distElemSub1*(0:noElemSub1-1)'*2*pi/lambda*sin(theta))) * ...
%                 (w2' * exp(-j*distElemSub2*(0:noElemSub2-1)'*2*pi/lambda*sin(theta))) * ...
%                 sqrt(cos(theta));
Expected_Value = exp(thSign*j*distElemMain*(0:noElemMain-1)'*2*pi/lambda*sin(theta)) * ...
                 (w1' * exp(thSign*j*distElemSub1*(0:noElemSub1-1)'*2*pi/lambda*sin(theta))) * ...
                 (w2' * exp(thSign*j*distElemSub2*(0:noElemSub2-1)'*2*pi/lambda*sin(theta))) * ...
                 sqrt(cos(theta));
Real_Result = spastevec(ant,theta,lambda);
disp(['Expected value < ' num2str(max(abs(Expected_Value))*8E-12)])
%disp('Expected value = 0')
disp(['Real value = ' num2str(sum(abs(Expected_Value - Real_Result)))])
% Expected result : 0
%   The value should be about the same and diff should be small
% Real result :



if(0) % Not needed. Use 13d instead!
% [ ]
disp('Test 13c')
% Comments
%  A quantitative test of a recursive ULA, where the
%  elements can be ULA:s or one-channel elements.
%  Only one of tests 13c and 13d need to be done.
  clear all
  global thSign
  thSign = 1;
  lambda = 0.2;
  elem    = defant('pattFuncElem', 'cos(x(1,:))');
  subarray = defant('ULA', [3, lambda/2], elem, [1 1 1].');
  totArray = defant('ULA', [24, 3*lambda/2], subarray, ones(24,1));
  figure
  s1 = pantpat2(totArray,[],[],[],[],[],lambda);
  ant1=totArray.element;
  figure
  s2 = pantpat2(ant1,[],[],[],[],[],lambda);
  ant2=ant1.element;
  figure
  s3 = pantpat2(ant2,[],[],[],[],[],lambda);

  %ls1=s1.specSmpl;
  %ls2=s2.specSmpl;
  %ls3=s3.specSmpl;
  %save tdbt13 s1 s2 s3
% Expected result :
%  With the aid of the following code we obtain graphs in which it
%  is possible to measure.
%  The graphs in the "Real result" section below must be in accordance
%  with the following graphs.
if (0)
  x = -90:2:90;
  y =cos(d2r(x))^2;
  plot(x,10*log10(y))

  lambda = 0.2;
  smallArray = defant('isotropULA', [3, lambda/2]);
  pantpat(smallArray,[],[],[],[],[],[],lambda);

  bigArray = defant('isotropULA', [24, 3*lambda/2]);
  pantpat2(bigArray,[],[],[],[],[],lambda);
    % This antenna pattern shall have grating lobes.
end%if (0)
% Real result :

end%if(0)




% [ ]
disp('Test 13d')
% Comments
%  A quantitative test of a recursive ULA, where the
%  elements can be ULA:s or one-channel elements.
%  This test uses results saved to file from test 13c
%  to decide success.
%  Only one of tests 13c and 13d need to be done.
% Expected result :
  clear all
  global thSign
  thSign = 1;
  load tdbt13
  ls1=s1;
  ls2=s2;
  ls3=s3;

  lambda = 0.2;
  elem    = defant('pattFuncElem', 'cos(x(1,:))');
  subarray = defant('ULA', [3, lambda/2], elem, [1 1 1].');
  totArray = defant('ULA', [24, 3*lambda/2], subarray, ones(24,1));
  figure, hold on
  s1 = pantpat2(totArray,[],[],[],[],[],lambda);
  ant1=totArray.element;
  s2 = pantpat2(ant1,[],[],[],[],[],lambda);
  ant2=ant1.element;
  s3 = pantpat2(ant2,[],[],[],[],[],lambda);
  difference = abs(sum(ls1.specSmpl-s1.specSmpl) + ...
    sum(ls2.specSmpl-s2.specSmpl) + sum(ls3.specSmpl-s3.specSmpl) );
  limit = eps;      % Ok on magic and curare.
  limit = 7.7e-11;  % Needed for hawk.
if ( difference  > limit)
  error('DBT-Error: Difference to large.')
else
  disp('Test Ok.')
end%if
% Real result :




elseif (testnr==14)
% *****************************************************************************
% [ ]  : Test 14
% Sign:     Datum:
%
% Test of the function dspastemat for recursive antenna definitions.
% ==Test of the function dspastemat for generalized ULAs
%
% *****************************************************************************
% Comments

% [ ]
% Comments
%  Test of the function dspastemat for a generalized ULA antenna.
disp('Test 14a')
format long;

n = 10;           % Index to channel to be calculated
m = 1;            % Index to target angle to be calculated

ant = defant('aimtEx');       % Application example antenna
theta = d2r([-20 -10 -5 -2 0 2 5 10 20]');      % Target angles. The no of targets is
phi = zeros(size(theta));     %  given by the no of target angles.
lambda   = 0.1;
doas = ([theta';phi']);       % Only theta is considered as target angle

AExp = dspastemat(ant,doas,lambda);

noElemMain  = 25;
noElemSub   = 3;
lambda   = 0.1;
distElemMain   = lambda;
distElemSub = lambda/2;
taper = [0.5 1 0.5]';
elem     = defant('pattFuncElem', 'sqrt(cos(x(1,:)))');