aimtexrun23.m

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

%AIMTEXRUN23 Monte Carlo simulations of number-of-targets estimation in the AIMT application example using one CPI only.
%
%Monte Carlo simulations and statistical analysis in the AIMT application example. One CPI. Three targets. Originally the same file as "run23" in the directory "lob/dgaantal".
%

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


%[noSrcEst, meanNoSrc, stdNoSrc, detProbEq, detProbGE, realNoSrc, noTrials, snr, elapTime] = function aimtexrun23


%tic
%t = clock;
startClock = clock;

more off
% ----------------------------------------------------------------------- %
% Parameters.
% ----------------------------------------------------------------------- %

noTrials = 100
%noTrials = 5
snr = [-50:5:-0].'
%snr = [-40,-10, 20].'
beampos         = d2r([-6:3:6]);
  % Do conventional beamforming in these directions.
doaPos	 = d2r([-7:0.1:7]);
  % Calculate spectrum in these directions for 'amir' etc.
rangePos = 13;


% ----------------------------------------------------------------------- %
% Definition of the receiver antenna.
% ----------------------------------------------------------------------- %

errAnt = defaimtexant2(1e-4*rand(1,51));
idealAnt = defaimtexant2(zeros(1,51));
%errAnt = idealAnt;
noChannels     = errAnt.noElem;        % Number of digital antenna channels.

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

prf = 8e3;
pri = 1/prf;
noPulses = 128;
noCPI = 1;
%fRF =    (2997:1:3004).'*1e6;
%fRF =    (3000:1:3001).'*1e6;
fRF =    (3000).'*1e6;
c0 = speedoflight;
lambda    = c0 ./ fRF;
sampleTime = 1e-6;   % SubPulseLength [s]. This is not used and not existing
  % in FA's code.
%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 sources.
% ----------------------------------------------------------------------- %

%target = deftarget(trgType, trgPower, trgRange, trgDOA, trgVel);
range4 = 1; % [m] This means that the target before pulse compression is
  % located in range bin 1 and after pulse compression with Barker 13 in
  % range bin 13.
fprintf('Target range %d [bins].\n',ran2ranbin([range4].', sampleTime, speedoflight, noRangeBins)+13-1)
doa1 = [d2r(0); 0];
doa2 = [d2r(1); 0];
doa3 = [d2r(2); 0];
trg1 = deftarget('swerling0', 1, range4, doa1, 21);
trg2 = deftarget('swerling0', 1, range4, doa2, 22);
trg3 = deftarget('swerling0', 1, range4, doa3, 23);

%srcDef = defsources(srcList, trgCorrMx)
targetSources = defsources({trg1 trg2 trg3});
realNoSrc = 3;

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

sig1 = simradarsig(errAnt, targetSources, waveform, pri);
sigsize(sig1)

figure

noMethods = 6;
noSNR = length(snr);
noSrcEst = zeros(noTrials,noMethods,noSNR);

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

for snrLoop = 1:noSNR
  %fprintf('snrLoop = %d (%d)\n',snrLoop, noSNR);
  for trialLoop =1:noTrials
    % ----------------------------------------------------------------------- %
    % Simulation of signals.
    % ----------------------------------------------------------------------- %

    noisePower = 10^(-snr(snrLoop)/10);
    noiseSrc1 = defnoise('Gaussian', noisePower);
    noiseSources = defsources({noiseSrc1});
    sig2 = simradarsig(errAnt, noiseSources, waveform, pri);
    sig  = addsig(sig1,sig2);
    sig3 = sig;

    % ----------------------------------------------------------------------- %
    % Introduce antenna element position errors.
    % ----------------------------------------------------------------------- %

    sig.antenna = idealAnt;

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

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

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

    % Digital beamforming
    sig = spafilt(sig,[],[],[],beampos,[],[],[],[],1);

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

    % Doppler filter bank
    % sig = dfb2(sig);

    % ----------------------------------------------------------------------- %
    % Estimation of number of sources.
    % ----------------------------------------------------------------------- %
    R = ecorrm(sig,[],rangePos);
    %diagLoad = R.corrMat(1,1);
    %R = ecorrm(R,'diagload',[],[],diagLoad);
    %spect1 = sdoaspc('music',R,doaPos);
    %splot2(spect1), drawnow

    %guessNoSrc = 4;
    guessNoSrc = 3;
    if (trialLoop == 1)

      %pulse = 1;
      %channel = 1:25;
      %noisePower = mean(abs(sig3.signals(pulse,15:32,channel)).^2);
      %tgt1Power = mean(abs(sig3.signals(pulse,1:13,channel)).^2);
      %SNR = squeeze(10*log10(tgt1Power./noisePower)).'
        % Assumes that the SNR is high.

      noChan = size(beampos,2);
      spect1 = sdoaspc('music',R,doaPos,guessNoSrc);
      figure(1)
      splot2(spect1)
      spec{snrLoop} = spect1;
      figure(6)
      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('Done: snrL=%d(%d), trialL=%d(%d), TotalL=%d(%d), E.time=%.3f(%.3f) hour.\r', snrLoop, noSNR, trialLoop, noTrials, totalLoopCount, ...
      noTotalLoop, elapTime/3600, estTotTime/3600);

  end%for trialLoopf
  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

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

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

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

save run23 noSrcEst meanNoSrc stdNoSrc detProbEq detProbGE realNoSrc ...
  noTrials snr elapTime

printfm run23a1 1 8
printfm run23a6 6 8

aimtexrun23plot

%quit	% Exists MATLAB and releases memory.


figure,sigplot2(sig3,1,':',1, {'scaleType', 'dB'})
pulse = 1;
channel = 7:8;
noisePower = mean(abs(sig3.signals(pulse,15:32,channel)).^2);
tgt1Power = mean(abs(sig3.signals(pulse,1:13,channel)).^2);
SNR = squeeze(10*log10(tgt1Power./noisePower)).'