modcorrmkernel.m

阵列信号处理的工具箱

function [ROut,T] = modcorrmkernel(RIn, method, in5, in6, in7)

%MODCORRMKERNEL Modifies a correlation matrix.
%
%--------
%Synopsis:
%  [ROut,T] = modcorrmkernel(RIn, method, in5, in6, in7)
%
%Description:
%  This function performs the actual modification of a correlation matrix.
%
%  This function is not aimed at the end user but at the extension programmer.
%
%Output and Input:
%  ROut (CxMatrixT): Output correlation matrix.
%  T (CxMatrixT): If the method changes the size of the correlation matrix,
%    this output parameter is the transformation matrix which is stored in
%    the output correlation matrix (RxCorrMatT) and specifies the change
%    of size.
%
%  RIn (CxMatrixT): Input correlation matrix.
%  method [D](StringT): Type of correlation matrix estimation. See below for
%    sections about each method.
%    = 'ml': ML-estimate (default). See below.
%    = 'forback': Forward-backward averaging. See below.
%    = 'spatsmooth': Spatial smoothing. See below.
%    = 'diagload': Diagonal loading. See below.
%    = 'diagloadeig': Diagonal loading relative noise eigenvalues. See below.
%    = 'equeig': Equal Noise Eigenvalues. See below.
%    = 'normeig': Norming relative the noise eigenvalues. See below.
%  in5, in6, in7 (): Depends on selected method.
%    See the help text of function "ecorrm" for sections about each method.
%
%--------
%Forward-Backward Averaging
%Synopsis:
%  Rout = ecorrm(inData, 'forback')
%  Rout = ecorrm(inData, 'forback', rangeIx, cpiIx)
%
%Description:
%  This method requires an ULA. When used on the 'expAnt' antenna type
%  this probably means that the calibration and near field compensations
%  must be performed on the signals (with "sigcomp2") and not on the
%  steering matrix (with "setcal1"). This method can only handle two
%  coherent signals. For more information about this method, see [2] p. 75-76.
%
%Algoritm:
%  Rfb = 0.5*(R+J*conj(R)*J);
%  where J is a (L x L) exchange matrix, whose components are zero except for
%  ones in the anti-diagonal, R is the ML-estimate, and Rfb is the modified
%  estimate.
%
%--------
%Spatial Smoothing
%Synopsis:
%  Rout = ecorrm(inData, 'spatsmooth', [], [], noElemSub)
%  Rout = ecorrm(inData, 'spatsmooth', [], [], noElemSub, channelIx)
%  Rout = ecorrm(inData, 'spatsmooth', rangeIx, cpiIx, noElemSub, channelIx)
%
%Output and Input:
%  noElemSub (IntScalarT): Number of elements in the subarrays.
%  channelIx [D](IndexT): Specifies which channels of the original steering
%    vector to use. The number of channels to use must be equal to
%    "noElemSub". If all channels are equal, which is a requirement for
%    this method, it does not matter which channels are chosen.
%    Default value is (1:noElemSub).
%
%Description:
%  In this method the array antenna is split into a number of overlapping
%  subarrays whose correlation matrices are averaged. The number K of
%  subarrays used is K = N - noElemSub +1, where N is the total number
%  of channels (elements) in the array. For each additional subarray the
%  method can handle one extra coherent source signal. A drawback with this 
%  method is that the size of the antenna decreases, since the subarrays
%  are smaller than the original array.
%
%  This method requires a regular array antenna with a translational
%  invariant property, e.g. an ULA. When used on the 'expAnt' antenna type
%  this probably means that the calibration and near field compensations
%  must be performed on the signals (with "sigcomp2") and not on the
%  steering matrix (with "setcal1").
%
%  For more information about this method, see [2] p. 75-76.
%
%Algoritm:
%  Rss = zeros(noElemSub,noElemSub);
%  for n=1:noSubArray
%    Rss = Rss + R(n:n+noElemSub-1, n:n+noElemSub-1);
%  end
%  Rout = (1/noSubArray)*Rss;
%  where noSubArray is the number of subarrays.
%
%--------
%Diagonal Loading
%Synopsis:
%  Rout = ecorrm(inData, 'diagload', [], [], loadVal)
%  Rout = ecorrm(inData, 'diagload', rangeIx, cpiIx, loadVal)
%
%Output and Input:
%  loadVal (RealScalarT):
%
%Description:
%  An amount, which is determined by the input parameter "loadVal", of
%  the identity matrix is added to the correlation matrix. Diagonal loading
%  reduces the fluctuation of the small eigenvalues of the correlation
%  matrix [5,3,4]. A suitable size of "loadVal" is in the order of the noise
% variance [5].
%
%Algoritm:
%
%--------
%Diagonal Loading relative noise eigenvalues
%Synopsis:
%  Rout = ecorrm(inData, 'diagloadeig', [], [], relLoadVal)
%  Rout = ecorrm(inData, 'diagloadeig', [], [], relLoadVal, noSrc)
%  Rout = ecorrm(inData, 'diagloadeig', rangeIx, cpiIx, relLoadVal)
%  Rout = ecorrm(inData, 'diagloadeig', rangeIx, cpiIx, relLoadVal, noSrc)
%
%Output and Input:
%  relLoadVal (RealScalarT): Load relative the noise eigenvalues (noise power).
%  noSrc [D](IntScalarT): Number of signal sources. This means that there
%    are ((number-of-channels) - noSrc) noise eigenvalues. If this value
%    is too low, the estimated noise power will be somewhat too low. If
%    this value is too high, the estimated noise power will be too high
%    (maybe much too high). Default value is = 1, which means that only the
%    smallest eigenvalue is used.
%
%Description:
%  An amount, which is determined by the input parameter "relLoadVal"
%  multiplied with the noise power, of the identity matrix is added to the 
%  correlation matrix. Diagonal loading reduces the fluctuation of the small 
%  eigenvalues of the correlation matrix [5,3,4]. A suitable size of 
%  "loadVal" is in the order of the noise variance (power) [5]. The noise
%  power is calculated as the mean value of the smallest eigenvalues of the 
%  spatial correlation matrix
%
%Algoritm:
%
%--------
%Equal Noise Eigenvalues.
%Synopsis:
%  Rout = ecorrm(inData, 'equeig')
%  Rout = ecorrm(inData, 'equeig', [], [], [], noSrc)
%  Rout = ecorrm(inData, 'equeig', rangeIx, cpiIx)
%  Rout = ecorrm(inData, 'equeig', rangeIx, cpiIx, [], noSrc)
%
%Output and Input:
%  noSrc [D](IntScalarT): Number of signal sources. This means that there
%    are ((number-of-channels) - noSrc) noise eigenvalues. If this value
%    is too low, the estimated noise power will be somewhat too low. If
%    this value is too high, the estimated noise power will be too high
%    (maybe much too high). Default value is = 1, which means that only the
%    smallest eigenvalue is used.
%
%Description:
%  All noise eigenvalues, as determined by "noSrc", are set equal to the
%  noise power. The noise power is calculated as the mean value of the
%  smallest eigenvalues of the spatial correlation matrix.
%
%Algoritm:
%
%--------
%Norming relative the noise eigenvalues.
%Synopsis:
%  Rout = ecorrm(inData, 'normeig')
%  Rout = ecorrm(inData, 'normeig', [], [], [], noSrc)
%  Rout = ecorrm(inData, 'normeig', rangeIx, cpiIx)
%  Rout = ecorrm(inData, 'normeig', rangeIx, cpiIx, [], noSrc)
%
%Output and Input:
%  relLoadVal (RealScalarT):
%  noSrc [D](IntScalarT): Number of signal sources. This means that there
%    are ((number-of-channels) - noSrc) noise eigenvalues. If this value
%    is too low, the estimated noise power will be somewhat too low. If
%    this value is too high, the estimated noise power will be too high
%    (maybe much too high). Default value is = 1, which means that only the
%    smallest eigenvalue is used.
%
%Description:
%  All eigenvalues (both signal and noise) are set equal to the noise power.
%  The noise power is calculated as the mean value of the smallest eigenvalues 
%  of the spatial correlation matrix.
%
%Algoritm:
%
%--------
%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:
%
%Software Quality:
%  (About what is done to ascertain software quality. What tests are done.)
%
%Known Bugs:
%
%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]: Krim H., Viberg M.:"Two Decades of Array Signal Processing Research",
%    IEEE Signal Processing Magazine, July 1996, pp. 67-94.
%  [3]: Carlsson B.D.: "Covariance Estimation Errors and Diagonal Loading
%    in Adaptive Arrays", IEEE Trans. AES-24, No. 4, July 1988, pp. 397-401.
%  [4]: Ganz M.W., Moses R.L., Wilson S.L.: "Convergence of the SMI and the
%    Diagonally Loaded SMI Algorithms with Weak Interference", IEEE Trans.
%    AP-38, No. 3, March 1990, pp. 394-399.
%  [5]: Nickel U.: "Comparison of some criteria to determine the dimension
%    of the signal subspace for small sample size", Report No. 457, FGAN-FFM
%    1996, Germany, p. 13.
%
%See Also:
%  eccorrm, basecorrm


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

  T = [];
  Rxx = RIn;
  noElem = size(Rxx,1);

  % ----------------------------------------------------------------------- %
  % Maximum Likelihood
  % ----------------------------------------------------------------------- %
  if (strcmp(method,'ml'))
    % Do nothing more.
    %disp('ml')

  % ----------------------------------------------------------------------- %
  % Forward-Backward Averaging
  % ----------------------------------------------------------------------- %
  elseif (strcmp(method,'forback'))
    %disp('forback')
    J   = fliplr(eye(size(Rxx)));
    Rxx = .5*(Rxx+J*conj(Rxx)*J);

  % ----------------------------------------------------------------------- %
  % Spatial Smoothing
  % ----------------------------------------------------------------------- %
  elseif (strcmp(method,'spatsmooth'))
    %disp('spatsmooth')
    noElemSub  = in5;
    channelIx = in6;
    if isempty(noElemSub)
      error('DBT-Error: Input argument noElemSub is undefined.');
    end%if
    if isempty(channelIx)
      channelIx = 1:noElemSub;
    end%if
  
    %antenna    = inData.antenna;
    %noElem     = antenna.noElem;
    noSubArray = noElem-noElemSub+1;
    Rss        = zeros(noElemSub,noElemSub);
  
    for n=1:noSubArray
      Rss = Rss + Rxx(n:n+noElemSub-1, n:n+noElemSub-1);
    end

    Rxx = (1/noSubArray)*Rss;

    T = zeros(noElem,noElemSub);
    T(channelIx,1:noElemSub) = eye(noElemSub);
  
    %antNew = defant('beamform',antenna,T);
      % Pick out only noElemSub elements as the new antenna.
    %%antenna.noElem = noElemSub;	% This is done by "defant" above.
    %Rout.antenna = antNew; 
  
    changeAntFlag = 1;


  % ----------------------------------------------------------------------- %
  % Diagonal Loading
  % ----------------------------------------------------------------------- %
  elseif (strcmp(method,'diagload'))
    %disp('diagload')
    diagload = in5;
    Rxx = Rxx + diagload*eye(size(Rxx));


  % ----------------------------------------------------------------------- %
  % Diagonal Loading relative the noise eigenvalues.
  % ----------------------------------------------------------------------- %
  elseif (strcmp(method,'diagloadeig'))
    %disp('diagloadeig')
    loadRelNoise = in5;
    %noChannels = inData.antenna.noElem;
    noChannels = noElem;
    noSrc = in6;
    if isempty(noSrc)
      noSrc = noChannels-1;
        % Default value: use only the smallest eigenvalue.
    end%if
    [U,S,V] = svd(Rxx);
    eigVals = diag(S);
      % Calculate eigenvalues of the spatial correlation matrix.
    noisePower = mean(eigVals(noSrc+1:noChannels));
    %fprintf('Noise power in ecorrm: %g dB. ', 10*log10(noisePower))
    diagload = loadRelNoise * noisePower;
    Rxx = Rxx + diagload*eye(size(Rxx));


  % ----------------------------------------------------------------------- %
  % Equal noise eigenvalues.
  % ----------------------------------------------------------------------- %
  elseif (strcmp(method,'equeig'))
    %disp('equeig')
    %noChannels = inData.antenna.noElem;
    noChannels = noElem;
    noSrc = in6;
    if isempty(noSrc)
      noSrc = noChannels-1;
        % Default value: use only the smallest eigenvalue.
    end%if
    [U,S,V] = svd(Rxx);
    eigVals = diag(S);
      % Calculate eigenvalues of the spatial correlation matrix.
    noisePower = mean(eigVals(noSrc+1:noChannels));
    eigVals(noSrc+1:noChannels)  = noisePower;
    Rxx = U*diag(eigVals)*V';


  % ----------------------------------------------------------------------- %
  % Norming relative the noise eigenvalues.
  % ----------------------------------------------------------------------- %
  elseif (strcmp(method,'normeig'))
    %disp('normeig')
    %noChannels = inData.antenna.noElem;
    noChannels = noElem;
    noSrc = in6;
    if isempty(noSrc)
      noSrc = noChannels-1;
        % Default value: use only the smallest eigenvalue.
    end%if
    [U,S,V] = svd(Rxx);
      % Calculate eigenvalues of the spatial correlation matrix.
    eigVals = diag(S);
    noisePower = mean(eigVals(noSrc+1:noChannels));
      % Choose the mean value of the smallest eigenvalues as an estimate of
      % the noise power.
    eigVals = eigVals / noisePower;
    Rxx = Rxx ./ noisePower;
      % Normalize the correlation matrix to the noise power.


  % ----------------------------------------------------------------------- %
  % The else branch
  % ----------------------------------------------------------------------- %
  else
    error('DBT-Error: Selected method is not implemented.')
  end%if

  ROut = Rxx;

%endfunction modcorrmkernel