calcfiltproc.m

阵列信号处理的工具箱

function [sigMatOut, basisMat, freqPos] = calcfiltproc(sigMatIn, method, weightMatrix,  taperCoeff, orthBeamFlag, extraParam, normType, normLevel, mcorrmParams)

%CALCFILTPROC Linear filtering.
%
%--------
%Synopsis:
%  [sigMatOut, basisMat, freqPos] = calcfiltproc(sigMatIn, method, weightMatrix,
%    taperCoeff, interCorrMat)
%
%Description:
%  Linear filtering of several signal vectors. 
%
%  This function must be used in conjunction with the function "calcfiltinit".
%  See the help text of "calcfiltinit" for more information.
%
%Output and Input:
%  sigMatOut (CxMatrixT): Output signal matrix. Rows and columns as for 
%    "sigMatIn".
%  basisMat (CxMatrixT): After the filtering, the signal consists of the 
%    coefficients in the basis, which is made up of the columns of this
%    output matrix. "basisMat" is the memory of what is done with the signal 
%    during the filtering. 
%    For filtering methods where the basis matrix is constant during the 
%    filtering the function "calcfiltinit" delivers the basis mattrix.
%    But for methods where the basis matrix is different in different parts
%    of the whole signal, this function "calcfiltproc" deliveres the 
%    previously used local basis matrix. 
%  freqPos (RealVectorT.'): Normalized frequencies of the output channels.
%    Not used.
%
%  sigMatIn (CxMatrixT): Input signal matrix. Different columns means
%    different signal vector samples or snapshots. Different rows means 
%    different elements in the same signal vector. For example with
%    spacial filtering (beamforming), different rows are different
%    antenna channels and different columns are different samples or
%    snapshots.
%  method (StringT):Filtering method. See the help text of the function
%    "calcfiltinit" for sections about each method.
%  weightMatrix (CxMatrixT): Filter coefficients which are calculated
%    by the function "calcfiltinit".
%  taperCoeff (CxVectorT): Weights to weight the filter coefficients with.
%  normType [D](StringT): Type of normalization. Not implemented.
%  mcorrmParams [D](CellArrayT): Parameters to send to the function
%    "" to modify the correlation matrix in methods, e.g. 'capon',
%    that uses a correlation matrix in "calcfiltproc". Not implemented.
%
%--------
%
%--------
%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.)
%  This function is not tested yet.
%
%Known Bugs:
%  There should be input and output parameters for "freqPos", "doaPos", etc.
%  Is there a bug in method 'trans'? Look at the code!
%
%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]: Johnson D.H., Dudgeon D.E.: "Array Signal Processing, Concepts
%    and Techniques", Prentice Hall 1993, ISBN 0-13-048513-6.
%
%See Also:
%  calcfiltinit

%   *  DBT, A Matlab Toolbox for Radar Signal Processing  *
% (c) FOA 1994-2000. See the file dbtright.m for copyright notice.
%
%  Start        : 980913 Svante Bj鰎klund (svabj).
%  Latest change: $Date: 2000/11/17 12:45:10 $ $Author: svabj $.
%  $Revision: 1.16 $
% *****************************************************************************

%  Move the filtering to after each method !?!


% ----------------------------------------------------------------------- %
% Handle input parameters
% ----------------------------------------------------------------------- %

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

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

arginNo=4;
if (nargin < arginNo)
  taperCoeff = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  extraParam = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  orthBeamFlag = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  extraParam = [];
end
arginNo = arginNo +1;
if (nargin < arginNo)
  normType = [];
end
arginNo = arginNo +1;


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

sigVectLen = size(sigMatIn,1);
noSamples = size(sigMatIn,2);

if isempty(taperCoeff)
  taperCoeff = ones(sigVectLen,1);
end%if
if isempty(orthBeamFlag)
  orthBeamFlag = 0;
end%if
if isempty(extraParam)
  extraParam = [];
end%if
if isempty(normType)
  normType = 'noNorm';
end%if


% ************* Pick out some more fields from input parameters. *************



% ****************** Error check input parameters ******************

%chkdtype(sigMatIn, 'CxMatrixT')
%chkdtype(method, 'StringT')
%chkdtype(weightMatrix, 'CxMatrixT')
%chkdtype(taperCoeff, 'CxVectorT')
  % The test is removed to save execution time. 
  % Some of the tests are already done in "calcfiltinit" and some we do
  % not care about.

% ----------------------------------------------------------------------- %
% Loops for several range gates and trials.
% ----------------------------------------------------------------------- %

% ****************** Create output variable. ******************


% ****************** Loops. ******************

%for pulseLoop = 1:length()
%  for rangeLoop = 1:length(rangeIx)
%  for trialLoop = 1:length(trialIx)



    % ----------------------------------------------------------------------- %
    % Code common to all methods.
    % ----------------------------------------------------------------------- %

    % **************** Pick out fields from input parameters. ****************


    % ****************** Precalculations. ******************


    % ----------------------------------------------------------------------- %
    % Methods
    % ----------------------------------------------------------------------- %

    % ----------------------------------------------------------------------- %
    % Adaptive Total Finite Response Filter (Adaptive Conventional Beamforming)
    %
    % Start: 981013 Svante Bj鰎klund (svabj).
    % ----------------------------------------------------------------------- %
    if (strcmp(method,'atfir') | strcmp(method,'tfir') | strcmp(method,'gcaap'))

      %fprintf('calcfiltproc: issparse(weightMatrix) = %d\n', ...
      %  issparse(weightMatrix));
      %if (issparse(weightMatrix))
      %   sigMatIn = sparse(sigMatIn);
      %end%if

      sigMatOut = weightMatrix' * sigMatIn;

      %fprintf('calcfiltproc: issparse(sigMatOut) = %d\n',issparse(sigMatOut));


    % ----------------------------------------------------------------------- %
    % Capons Method.
    %
    % Start: 981013 Svante Bj鰎klund (svabj).
    % ----------------------------------------------------------------------- %
    elseif (strcmp(method,'capon'))
      steMat = weightMatrix;
      sigCorrMat = (1/noSamples) * sigMatIn*sigMatIn';
 
      if (0)
        % This version can be used for execution time measurements.
        %figure
        %s = svd(sigCorrMat);
        %plot(10*log10(s))
        d4 = sigCorrMat \ steMat;
        d1 = steMat' * d4;
          % This (above) line is time consuming. 
        d2 = diag(d1);
        normFactors = 1 ./ d2;
        %d3 = diag(normFactors);
        %W = (d4) * d3;
        W = weightcols(d4,normFactors);
        %difference = sum(sum(abs(W - W2)));
        %limit = 0.5*size(W,1)*size(W,2)*eps;
        %if (difference > limit)
        %  difference
        %  error('Difference to large.')
        %else
        %  disp('  Test Ok.')
        %end%if
      else
        d4 = sigCorrMat \ steMat;
        W = weightcols(d4,(1 ./ diag(steMat' * d4)));
          % This line is time consuming. In a run with ceship.m 74% of
          % the time was spent on this line.
      end%if
      if (orthBeamFlag), W = orthbeam(W); end%if
      sigMatOut = W' * sigMatIn;
      basisMat =  W;
      %freqPos = ...


    % ----------------------------------------------------------------------- %
    % Discrete Fourier Transform Method.
    %
    % Start: 981013 Svante Bj鰎klund (svabj).
    % ----------------------------------------------------------------------- %
    elseif (strcmp(method,'fft'))
      if (0)
        if isempty(extraParam)
          extraParam = size(sigMatIn,1);
        end%if
        noPulses = size(sigMatIn,1);
        noDopChan = extraParam;
      end%if (0)

      noPulses = size(weightMatrix,2);
      noDopChan = size(weightMatrix,1);
      sigMatIn = diag(taperCoeff) * sigMatIn;
      % The beams in FFT are always orthogonal. No need for an explicit
      % orthogonalization.
      sigMatOut = fftshift1(fft(sigMatIn, noDopChan));

      %basisMat =  fftshift1(dftopm(noDopChan,noPulses)');
      %freqPos = ...


    % ----------------------------------------------------------------------- %
    % Inverse Discrete Fourier Transform Method.
    %
    % Start: 98xxxx Svante Bj鰎klund (svabj).
    % ----------------------------------------------------------------------- %
    elseif (strcmp(method,'ifft'))
      if (0)
      if isempty(extraParam)
        extraParam = size(sigMatIn,1);
      end%if
      noPulses = extraParam;
      noDopChan = size(sigMatIn,1);
      end%if (0)

      noPulses = size(weightMatrix,1);
      noDopChan = size(weightMatrix,2);
      sigMatIn = diag(taperCoeff) * sigMatIn;
      % The beams in IFFT are always orthogonal. No need for an explicit
      % orthogonalization.
      sigMatOut = fftshift1(ifft(sigMatIn, noPulses),'inv');

      %basisMat =  fftshift1(dftopm(noDopChan,noPulses)');
      %freqPos = ...



    % ----------------------------------------------------------------------- %
    % Direct Transformation Matrix Method.
    % (Direct beamforming weight multiplication)
    %
    % Start: 981013 Svante Bj鰎klund (svabj).
    % ----------------------------------------------------------------------- %
    elseif (strcmp(method,'trans'))
      transMat = weightMatrix;
      sigMatIn = diag(taperCoeff) * sigMatIn;  % Shouldn't it be "transMat"
        % instead of "sigMatIn" ?
      %if (orthBeamFlag), transMat = orthbeam(transMat); end%if
      sigMatOut = transMat' * sigMatIn;


    % ----------------------------------------------------------------------- %
    % SB Interference Eigenvector Method.
    % (Using the eigenvectors as weights)
    %
    % Implementation is not finnished.
    %
    % Start: 991017 Svante Bj鰎klund (svabj).
    % ----------------------------------------------------------------------- %
    elseif (strcmp(method,'intereigvec'))
      sigMatOut = weightMatrix' * sigMatIn;

  
    % ----------------------------------------------------------------------- %
    % SB Data Eigenvector Method.
    % (Using the eigenvectors as weights)
    %
    % Implementation is not finnished.
    % There should be a special treatment, as with 'fft', of this method
    % in "puldfilt" etc. 
    % because the number of output channels cannot be chosen by "freqChoice".
    %
    % Start: 991017 Svante Bj鰎klund (svabj).
    % ----------------------------------------------------------------------- %
    elseif (strcmp(method,'dataeigvec'))
        % The calculation is the same as in 'intereigvec' but the correlation 
        % matrix is formed from the signal samples instead from the input 
        % interference correlation matrix.
      %steMat = weightMatrix; % Not used.
      sigCorrMat = (1/noSamples) * sigMatIn*sigMatIn';
      [W,S,V]=svd(sigCorrMat);
      if (orthBeamFlag), W = orthbeam(W); end%if
      sigMatOut = W' * sigMatIn;
      basisMat =  W;


    % ----------------------------------------------------------------------- %
    % MUSIC.
    % 
    % Note, this method is non-linear. Other linear processing cannot be 
    % performed after this method. This method should be in the function 
    % "calcspcproc" instead.
    %
    % Start: 990307 Svante Bj鰎klund (svabj).
    % ----------------------------------------------------------------------- %
    elseif (strcmp(method,'music'))
      if isempty(extraParam)
        extraParam = fix(size(sigMatIn,1)/2);
      end%if
      noSrc = extraParam;
      noChan = size(sigMatIn,1);
      steMat = weightMatrix;
      L = size(sigMatIn,2);
      Rxx = 1/L * sigMatIn*sigMatIn';
      [U,S,V]=svd(Rxx);   % Eigendecomp. of correlation matrix using SVD.
      En=U(:,(noSrc+1):noChan);  % Splitting in signal and noise subspace.
      %Pspect = sum((abs(steMat)).^2)./(eps+sum((abs(En.'*conj(steMat))).^2)); 
      Pspect = sqrt(sum((abs(steMat)).^2) ./ ...
        (eps+sum((abs(En.'*conj(steMat))).^2))); 
        % eps to avoid divide by zero.
        % sqrt(.) because the output of this function is supposed to be 
        % amplitudes.
      sigMatOut = repmat(Pspect.',1,L);
      basisMat =  weightMatrix;
        % The output parameter "weightMatrix" is nonsense for this method.
      %freqPos = ...


    % ----------------------------------------------------------------------- %
    %
    % ----------------------------------------------------------------------- %
    elseif (strcmp(method,'dummy'))


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



    % ----------------------------------------------------------------------- %
    % Return data.
    % ----------------------------------------------------------------------- %