run8cpi.m

雷达信号处理、或阵列信号处理中能够用上的重要的matlab工具箱——阵列信号处理工具箱

function [noSrcEst, meanNoSrc, stdNoSrc, detProbEq, detProbGE, realNoSrc, spec, noTrials, snr, elapTime] =  run8cpi(simName, broadside, noTgt, noTrials, snr, startFigNr, dummy)

%RUN8CPI Simulations and statistical analysis of number-of-targets estimation in the AIMT application example with several CPIs.
%
%--------
%Synopsis:
%  [noSrcEst, meanNoSrc, stdNoSrc, detProbEq, detProbGE, realNoSrc, spec, 
%    noTrials, snr,elapTime] =  run8cpi(simName, broadside, noTgt, noTrials, 
%    snr, startFigNr, dummy)
%
%Description:
%  Monte Carlo simulations and statistical analysis of number-of-targets 
%  estimation in the AIMT application example [2] with several CPIs with
%  frequency switching. Different carrier frequencies are used in the
%  different CPIs.  Targets directions are 0, 1 & 2 or 34, 35 & 36 degrees.
%  Target simulation, beamforming and estimation of the number of targets
%  are using subarrays scanned to 0 or 30 degrees.
%
%  No figure windows can be open when this function is invoked.
%
%Output and Input:
%
%  simName (StringT): First part of file names used for saving simulation 
%    results and figures.
%  broadside (BoolT): If == 1, then the targets DOAs  (0, 1 & 2 degrees)
%    and the examination direction of the antenna (about -7 to +7 degrees)
%    is broadside, i.e. perpendicular to the antenna array.
%    If == 0, then the targets DOAs  (34, 35 & 36 degrees)
%    and the examination direction of the antenna (about 23 to 37 degrees)
%    is non-broadside.
%  noTgt (IntScalarT): Number of targets present. Possible values are 1,2 & 3.
%  noTrials (IntScalarT): Number of trials for each SNR to simulate.
%  snr (RealVectorT): A vector with the signal-to-noise (SNR) values to
%    use in the simulation [dB]. For each SNR values, noTrials trials will
%    be executed and statistical analysis (mean and standard deviation)
%    will be performed.
%  startFigNr (IntScalarT): The first figure number to use when opening
%    figure windows. 
%
%--------
%Notations:
%  Data type names are shown in parentheses and they start with a capital
%  letter and end with a capital T. Data type definitions can be found in [1]
%  or by "help dbtdata".
%  [D] = This parameter can be omitted and then a default value is used.
%  When the [D]-input parameter is not the last used in the call, it must be
%  given the value [], i.e. an empty matrix.
%  ... = There can be more parameters. They are explained under respective
%  metod or choice.
%
%Examples:
%  noTrial = 50;
%  snr = [-50:5:10].';
%  
%  % Broadside
%  run8cpi('run42',1,1,noTrial,snr,1)
%  run8cpi('run43',1,2,noTrial,snr,11)
%  run8cpi('run44',1,3,noTrial,snr,21)
%
%  % Non-broadside
%  run8cpi('run45',0,1,noTrial,snr,31)
%  run8cpi('run46',0,2,noTrial,snr,41)
%  run8cpi('run47',0,3,noTrial,snr,51)
%  
%  pause
%  run8cpiplot('run42', 2,'onlySave')
%  run8cpiplot('run43',12,'onlySave')
%  run8cpiplot('run44',22,'onlySave')
%  
%  run8cpiplot('run45',32,'onlySave')
%  run8cpiplot('run46',42,'onlySave')
%  run8cpiplot('run47',52,'onlySave')
%
%Software Quality:
%  (About what is done to ascertain software quality. What tests are done.)
%  This function is not tested yet.
%
%Known Bugs:
%  Should use the antenna model without positioning errors in the conventional 
%  processing.
%
%References:
%  [1]: Bj鰎klund S.: "DBT, A MATLAB Toolbox for Radar Signal Processing.
%    Reference Guide", FOA-D--9x-00xxx-408--SE, To be published.
%  [2]: Athley F.:"Till鋗pningsexempel f鰎 studie av modellbaserad
%    lobformning", Ericsson Microwave Systems 1996, No FY/XU-96:299.In Swedish.
%
%See Also:
%  run8cpiplot


%
%     *    DBT, A Matlab Toolbox for Radar Signal Processing    *
% (c) FOA 1994-2000. See the file dbtright.m for copyright notice.
%
%    Start                : 990115 Svante Bj鰎klund (svabj).
%    Latest change: $Date: 2000/10/16 15:39:24 $ $Author: svabj $.
%    $Revision: 1.10 $
% *****************************************************************************

%tic
startClock = clock;
more off

if (nargin < 1)
  error('DBT-Error: To few input parameters.')
end

% ****************** Add missing input parameters ******************

arginNo=2;
if (nargin < arginNo)
  broadside = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  noTgt = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  noTrials = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  snr = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  startFigNr = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  dummy = [];
end
arginNo = arginNo +1;

% ****************** Default values ******************

if isempty(broadside)
  broadside = 1;
end%if
if isempty(noTgt)
  noTgt = 2;
end%if
if isempty(noTrials)
  noTrials = 3
end%if
if isempty(snr)
  %snr = [-50:5:-0].'
  snr = [-50,-0].'
end%if
if isempty(startFigNr)
  startFigNr = 1;
end%if
if isempty(startFigNr)
  dummy = 0;
end%if

% ----------------------------------------------------------------------- %
% Parameters.
% ----------------------------------------------------------------------- %

%noTrials = 10
%noTrials = 5
%snr = [-50:5:-0].'
%snr = [-70 -10].'

fprintf('noTrials = %d\n',noTrials)
fprintf('snr = ['), fprintf('%d ',snr.'), fprintf(']\n')
fprintf('broadside = %d\n',broadside)

if (broadside)
  beampos         = d2r([-6:3:6]);
    % Do conventional beamforming in these direction.
  subArrPoint = d2r([0; 0]);
    % Pointing direction of subarrays.
  doaPos	= d2r([-7:0.1:7]);
    % Calculate spectrum in these direction for 'amir' etc.
else
  beampos         = d2r([24:3:36]);
    % Do conventional beamforming in these direction.
  subArrPoint = d2r([30; 0]);
    % Pointing direction of subarrays.
  doaPos	 = d2r([22:0.1:38]);
    % Calculate spectrum in these direction for 'amir' etc.
end%if (broadside)

rangePos = 13;
  % Target range [range bins] after pulse compression.
dopplerPos = 72;	
  % The targets are located in this doppler channel. This is aproximatelly
  % a doppler frequency of 470 Hz.
freqNo = 4;
  % Choose wavelength ("freqNo") number 4 (3.00 GHz) as a compromise for
  % the signal processing.


% ----------------------------------------------------------------------- %
% Definition of the waveform.
% ----------------------------------------------------------------------- %

prf = 8e3;
pri = 1/prf;
noPulses = 128;
pulAxis = linspace(-prf/2,prf/2,noPulses);
noCPI = 8;
fRF =    (2997:1:3004).'*1e6;
%fRF =    (3000).'*1e6;
c0 = speedoflight;
lambda    = c0 ./ fRF;
sampleTime = 1e-6;   % SubPulseLength [s]
%noRangeBins = round(pri / sampleTime);
noRangeBins = 32;
pModulation = getmod('barker','13');  % Barker-13 has sidelobes 1/13 of the main peak.
%waveform = defwave(lambda, noRangeBins, noPulses, pModulation, sampleTime, ...
%  noCPI);

% ----------------------------------------------------------------------- %
% Definition of the receiver antenna.
% ----------------------------------------------------------------------- %
posErr = 1e-4*rand(1,51);	% Positioning errors of the antenna elements.
%errAnt = defaimtexant2(posErr,subArrPoint, lambda(freqNo));
noChannels = 25;        % Number of digital antenna channels.

% ----------------------------------------------------------------------- %
% Definition of the sources.
% ----------------------------------------------------------------------- %

%target = deftarget(trgType, trgPower, trgRange, trgDOA, trgVel);
range1 = 1; % [m]
range2 = 31; % [m]
range3 = 61; % [m]
  % This means that the targets before pulse compression are
  % located in range bin 1 and after pulse compression with Barker 13 in
  % range bin 13.
  % The size of one range bin is (c*sampletime)/2 = 150 m.
if (broadside)
  doa1 = d2r([0; 0]);
  doa2 = d2r([1; 0]);
  doa3 = d2r([2; 0]);
else
  doa1 = d2r([34; 0]);
  doa2 = d2r([35; 0]);
  doa3 = d2r([36; 0]);
end%if (broadside)

trg1 = deftarget('swerling0', 1, range1, doa1, 21);
trg2 = deftarget('swerling0', 1, range2, doa2, 22);
trg3 = deftarget('swerling0', 1, range3, doa3, 23);

fprintf('noTgt = %d\n',noTgt)
%srcDef = defsources(srcList, trgCorrMx)
if (noTgt == 1)
  targetSources = defsources({trg1});
  ranges = [range1];
elseif (noTgt == 2)
  targetSources = defsources({trg1 trg2});
  ranges = [range1 range2];
elseif (noTgt == 3)
  targetSources = defsources({trg1 trg2 trg3});
  ranges = [range1 range2 range3];
end%if
fprintf('Target range after pulse compression: %d [bins].\n',ran2ranbin(ranges, sampleTime, speedoflight, noRangeBins)+13-1)
realNoSrc = noTgt;

% ----------------------------------------------------------------------- %
% Simulation of constant signals.
% ----------------------------------------------------------------------- %

waveformProc = defwave(lambda(freqNo), noRangeBins, noPulses, pModulation, ...
  sampleTime, noCPI);
  % Waveform definition used in the signal processing. Waveform for all CPI:s.
antennaProc = defaimtexant2([],subArrPoint, lambda(freqNo));
  % Antenna definition used in the signal processing.

sig1 = RxRadarSigT(zeros(noPulses, noRangeBins, ...
  noChannels,1,noCPI,1),{antennaProc, waveformProc});
  % In "sig1" the constant target signals are stored. Different CPI:s of the target signals are simulated using different wavelengths.
sigParamProc = RxRadarSigT(zeros(1,1,1,1,1,1),{antennaProc, waveformProc});


for cpiLoop = 1: noCPI
  errAntCpi = defaimtexant2(posErr, subArrPoint, lambda(cpiLoop));
  noCpiPerCpi = 1;
  waveformCpi = defwave(lambda(cpiLoop), noRangeBins, noPulses, ...
    pModulation, sampleTime, noCpiPerCpi);
      % AntDefT and WaveformT for only one CPI. In this way different
      % wavelengths can be used for different CPIs.
  %sig11 = simradarsig(errAntCpi, targetSources, waveformCpi, pri, [], ...
  %  {subArrPoint,[]});
  sig11 = simradarsig(errAntCpi, targetSources, waveformCpi, pri, []);
  sig1.signals(:,:,:,:,cpiLoop,:) = sig11.signals;
end%for cpiLoop
  % This loop can also be done directly with the following code:

noMethods = 6;
noSNR = length(snr);
noSrcEst = zeros(noTrials,noMethods,noSNR);
spec = cell(noSNR);
%figure

estTotTime = Inf;
noTotalLoop = noSNR*noTrials;
clockPreLoop = clock;
elapTimePreLoop =  etime(clockPreLoop,startClock);

for snrLoop = 1:noSNR
  for trialLoop =1:noTrials
    % ----------------------------------------------------------------------- %
    % Simulation of noise signals.
    % ----------------------------------------------------------------------- %

    %noiseSrc = defnoise(noiseType, noisePower, corrMx);
    noisePower = 10^(-snr(snrLoop)/10);
    noiseSrc1 = defnoise('Gaussian', noisePower);
    noiseSources = defsources({noiseSrc1});
    %sig2 = simradarsig(errAnt, noiseSources, waveformProc, pri,[] , ...
    %  {subArrPoint,[]});
    sig2 = simradarsig(antennaProc, noiseSources, waveformProc, pri,[]);
      % The antenna definition is not used when simulating noise.
      % "sig2" contains the simulated noise s韌nals.
    sig  = addsig(sig1,sig2);
      % Add noise to the target signals. AntDefT, WaveformT etc. are copied
      % from "sig1" but this has no effect because only the field "signals"
      % is used.

    snapshots= zeros(size(beampos,2),noCPI);
    for cpiLoop = 1: noCPI
      antennaCpi = defaimtexant2(posErr, subArrPoint, lambda(cpiLoop));
      noCpiPerCpi = 1;
      waveformCpi = defwave(lambda(cpiLoop), noRangeBins, noPulses, ...
        pModulation, sampleTime, noCpiPerCpi);
        % AntDefT and WaveformT for only one CPI. In this way different
        % wavelengths can be used for different CPIs.
      sig22 = RxRadarSigT(sig.signals(:,:,:,:,cpiLoop,:), ...
        {antennaCpi,waveformCpi});
        % Pick out radar signal for only one CPI.

      % --------------------------------------------------------------------- %
      % Simulation of channel errors.
      % --------------------------------------------------------------------- %

      sig22 = simcherr(sig22,0.1,d2r(3),0.1);

      % --------------------------------------------------------------------- %
      % Signal processing.
      % --------------------------------------------------------------------- %

      % Digital beamforming
      %[sig22] = spafilt(sig22,[],[],{subArrPoint,[]},beampos,[],[],[],[],1);
      [sig22, spaTrans] = spafilt(sig22,[],[],[],beampos,[],[],[],[],1);

      % Pulse compression.
      sig22 = pulscomp(sig22, pModulation);

      % Doppler filter bank
      sig22 = dfb2(sig22);
      sig3 = sig22;
      snapshots(:,cpiLoop) = getm3(sig22.signals,3,[],...
        dopplerPos,rangePos,':',1,1);
    end%for cpiLoop

    % ----------------------------------------------------------------------- %
    % Estimation of number of sources.
    % ----------------------------------------------------------------------- %

    %R = basecorrm(getm3(sig.signals,3,[],dopplerPos,rangePos,':',1,':'), sig22)
    %R = basecorrm(snapshots, sig22);
    sigParamProc2 = setspatrans(sigParamProc, spaTrans);
    R = basecorrm(snapshots, sigParamProc2);

    %guessNoSrc = max(fix(noChan/2), noChan-4);
    guessNoSrc = 3;

    if (trialLoop == 1)
      noChan = size(beampos,2);
      spect1 = sdoaspc('music',R,doaPos,guessNoSrc);
      figure(1+startFigNr-1)
      splot2(spect1)
      spec{snrLoop} = spect1;
      figure(6+startFigNr-1)
      eigplot1(R,'dB')
      drawnow
    end%if

    noSrc1 = spanosrc(R,'aic',[],[],guessNoSrc,doaPos);
    noSrc2 = spanosrc(R,'mdl',[],[],guessNoSrc,doaPos);
    %noSrc3 = spanosrc(R,'amir',[],[],guessNoSrc,doaPos);
    noSrc4 = spanosrc(R,'amirbin',[],[],guessNoSrc,doaPos);
    %noSrc5 = spanosrc(R,'MA1',[],[],guessNoSrc,doaPos);
    noSrc6 = spanosrc(R,'fam',[],[],guessNoSrc,doaPos);
    noSrcEst(trialLoop,1,snrLoop) = noSrc1;
    noSrcEst(trialLoop,2,snrLoop) = noSrc2;
    %noSrcEst(trialLoop,3,snrLoop) = noSrc3;
    noSrcEst(trialLoop,4,snrLoop) = noSrc4;
    %noSrcEst(trialLoop,5,snrLoop) = noSrc5;
    noSrcEst(trialLoop,6,snrLoop) = noSrc6;

    totalLoopCount = (snrLoop-1)*noTrials+trialLoop;
    elapTime = etime(clock,startClock);
    if ((elapTime > 300) | (totalLoopCount >= 5))
      estTotTime = ((elapTime-elapTimePreLoop) * ...
        (noTotalLoop/totalLoopCount)) + elapTimePreLoop;
      %estTotTime = ((elapTime) * (noTotalLoop/totalLoopCount));
    end%if
    fprintf('%s:Done snrL=%d(%d),trialL=%d(%d),TotalL=%d(%d),E.time=%.3f(%.3f) hour.\r', simName,snrLoop, noSNR, trialLoop, noTrials, totalLoopCount, ...
      noTotalLoop, elapTime/3600, estTotTime/3600);

  end%for trialLoop
  fprintf('\n')
  fprintf('SNR = %d.\',snr(snrLoop))
  fprintf('\nMean of estimated number of sources.\n')
  meanNoSrc(snrLoop,:) = mean(squeeze(noSrcEst(:,:,snrLoop)))
  fprintf('\nStandard deviation of estimated number of sources.\n')
  stdNoSrc(snrLoop,:) = std(squeeze(noSrcEst(:,:,snrLoop)))
end%for snrLoop

B = (noSrcEst == realNoSrc);
detProbEq = (squeeze(sum(B))/noTrials).';

B = (noSrcEst >= realNoSrc);
detProbGE = (squeeze(sum(B))/noTrials).';

elapTime = etime(clock,startClock);	% Time in seconds.
fprintf('Elapsed time %d s.\n',elapTime)

eval(['save ', simName, ' noSrcEst meanNoSrc stdNoSrc detProbEq detProbGE realNoSrc spec noTrials snr elapTime'])

printfm([simName,'a1'], 1+startFigNr-1, 8)
printfm([simName,'a6'], 6+startFigNr-1, 8)
%printfm run21a1 1 8

run8cpiplot(simName,startFigNr+1)
%run21plot