% FROMCSD Construct a matrix filter from cross-spectral density matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION:    fromCSD
%
% DESCRIPTION: Construct matrix filter from cross-spectral density. Such a
% filter can be used for multichannel noise generation in combination with
% the MultiChannelNoise constructor of the matrix class.
%
%
% CALL:        a = fromCSD(a, pl)
%
% PARAMETER:   pl: plist containing 'csd'
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = fromCSD(a, pli)
  
  VERSION = '$Id: fromCSD.m,v 1.8 2010/10/29 16:09:12 ingo Exp $';
  
  % get AO info
  iobj = matrix.getInfo('matrix', 'From CSD');
  
  % Set the method version string in the minfo object
  iobj.setMversion([VERSION '-->' iobj.mversion]);
  
  % Add default values
  pl = parse(pli, iobj.plists);
  
  % Get parameters and set params for fit
  csdm        = find(pl, 'csd');
  if isa(csdm,'matrix') % get elements out of the input matrix
    csdao = csdm.objs;
  else
    csdao = csdm;
  end
  fs          = find(pl, 'fs');
  
  target      = lower(find(pl, 'targetobj')); % decide to perform s domain or z domain identification
  % if target is parfrac output a matrix of parfarc objs (s domain
  % identification)
  % if target is miir output a matrix of filterbank parallel miir objects
  % (z domain identification)
  
  usesym = lower(find(pl, 'UseSym'));
  
  if  (fs == 0) && strcmpi(target,'miir')
    error('### Please provide a valid sampling frequency for CSD constructor.');
  elseif isempty(fs) && strcmpi(target,'miir')
    error('### Please provide a valid sampling frequency for CSD constructor.');
  end
  
  % get units for filters
  tgiunit = find(pl,'iunits');
  tgounit = find(pl,'ounits');
  
  params = struct();
  
  params.Nmaxiter = find(pl, 'MaxIter');
  params.minorder = find(pl, 'MinOrder');
  params.maxorder = find(pl, 'MaxOrder');
  params.spolesopt = find(pl, 'PoleType');
  params.weightparam = find(pl, 'Weights');
  
  % set the target output
  if strcmpi(target,'miir')
    params.TargetDomain = 'z';
  elseif strcmpi(target,'parfrac')
    params.TargetDomain = 's';
  else
    error('### Unknown option for ''targetobj''.');
  end
  
  % Tolerance for MSE Value
  lrscond = find(pl, 'FITTOL');
  % give an error for strange values of lrscond
  if lrscond<0
    error('### Negative values for FITTOL are not allowed. ')
  end
  % handling data
  lrscond = -1*log10(lrscond);
  % give a warning for strange values of lrscond
  if lrscond<0
    warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond)))
  end
  params.lrscond = lrscond;
  
  % Tolerance for the MSE relative variation
  msevar = find(pl, 'MSEVARTOL');
  % handling data
  msevar = -1*log10(msevar);
  % give a warning for strange values of msevar
  if msevar<0
    warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar)))
  end
  params.msevar = msevar;
  
  if isempty(params.msevar)
    params.ctp = 'chival';
  else
    params.ctp = 'chivar';
  end
  
  if(find(pl, 'plot'))
    params.plot = 1;
  else
    params.plot = 0;
  end
  
  params.fs = fs;
  params.dterm = 0; % it is better to fit without direct term
  
  % check if symbolic calculation is required
  if strcmpi(usesym,'on')
    params.usesym = 1;
  elseif strcmpi(usesym,'off')
    params.usesym = 0;
  else
    error('### Unknown option for ''UseSym''.');
  end
  
  % extracting csd
  if numel(csdao)==1 % one dimensional psd
    csd = csdao.y;
    freq = csdao.x;
    dim = 'one';
  else % multichannel system
    dim = 'multi';
    [nn,mm] = size(csdao);
    if nn~=mm
      error('### CSD Matrix must be square. ')
    end
    freq = csdao.index(1).x;
    for ii = 1:nn
      for jj = 1:nn
        tcsd = csdao.index(ii,jj).y;
        % willing to work with columns
        [aa,bb] = size(tcsd);
        if aa<bb
          tcsd = tcsd.';
        end
        csd(ii,jj,:) = tcsd;
      end
    end
    
  end
  
  
  % call csd2tf
  % ostruct is a struct array whose fields contain the residues and poles
  % of estimated TFs. Since the fit is porformed on the columns of the TF
  % matrix, each element of the array contains the residues and poles
  % corresponding to the functions on the given column of the TF matrix.
  
  %ostruct = utils.math.csd2tf(csd,freq,params);
  
  ostruct = utils.math.csd2tf2(csd,freq,params);
  
  
  % the filter for each channel is implemented by the rows of the TF matrix
  
  switch dim
    case 'one'
      
      switch target
        case 'miir'
          % --- filter ---
          res = ostruct.res;
          poles = ostruct.poles;
          % construct a struct array of miir filters vectors
          pfilts(numel(res),1) = miir;
          for kk=1:numel(res)
            ft = miir(res(kk), [ 1 -poles(kk)], fs);
            ft.setIunits(unit(tgiunit));
            ft.setOunits(unit(tgounit));
            pfilts(kk,1) = ft;
            clear ft
          end
          filt = filterbank(pfilts,'parallel');
          
          a = matrix(filt);
          
          % Add history
          a.addHistory(iobj, pl, [], [csdao(:).hist]);
          
        case 'parfrac'
          res = ostruct.res;
          poles = ostruct.poles;
          
          fbk = parfrac(res,poles,0);
          fbk.setIunits(unit(tgiunit));
          fbk.setOunits(unit(tgounit));
          
          a = matrix(fbk);
          
          % Add history
          a.addHistory(iobj, pl, [], [csdm(:).hist]);
      end
      
      
    case 'multi'
      
      switch target
        case 'miir'
          % init filters array
          %fbk(nn*nn,1) = filterbank;
          %fbk = filterbank.newarray([nn nn]);
          
          for zz=1:nn*nn % run over system dimension
            % --- get column filter coefficients ---
            % each column of mres\mpoles are the coefficients of a given filter
            clear res poles
            res = ostruct(zz).res;
            poles = ostruct(zz).poles;
            
            % construct a struct array of miir filters vectors
            %ft(numel(res),1) = miir;
            for kk=1:numel(res)
              ft(kk,1) = miir(res(kk), [1 -poles(kk)], fs);
              ft(kk,1).setIunits(unit(tgiunit));
              ft(kk,1).setOunits(unit(tgounit));
            end
            
            fbk(zz,1) = filterbank(ft,'parallel');
            clear ft
            
          end
          
          mfbk = reshape(fbk,nn,nn);
          
          a = matrix(mfbk);
          
          %           % overide default name
          %           if strcmpi(pl.find('name'), 'None')
          %             pl.pset('name', sprintf('fromCSD(%s)', csdao.name));
          %           end
          %           % overide default description
          %           if isempty(pl.find('description'))
          %             pl.pset('description', csdao.description);
          %           end
          
          % Add history
          a.addHistory(iobj, pl, [], [csdao(:).hist]);
          
        case 'parfrac'
          % init filters array
          %fbk(nn*nn,1) = parfrac;
          
          for zz=1:nn*nn % run over system dimension
            % --- get column filter coefficients ---
            % each column of mres\mpoles are the coefficients of a given filter
            clear res poles
            res = ostruct(zz).res;
            poles = ostruct(zz).poles;
            
            fbk(zz,1) = parfrac(res,poles,0);
            fbk(zz,1).setIunits(unit(tgiunit));
            fbk(zz,1).setOunits(unit(tgounit));
            
          end
          
          mfbk = reshape(fbk,nn,nn);
          
          a = matrix(mfbk);
          
          %           % overide default name
          %           if strcmpi(pl.find('name'), 'None')
          %             pl.pset('name', sprintf('fromCSD(%s)', csdao.name));
          %           end
          %           % overide default description
          %           if isempty(pl.find('description'))
          %             pl.pset('description', csdao.description);
          %           end
          
          % Add history
          a.addHistory(iobj, pl, [], [csdao(:).hist]);
      end
      
  end
  
  % Set properties from the plist
  warning('off', utils.const.warnings.METHOD_NOT_FOUND);
  % remove parameters we already used
  pl_set = copy(pl,1);
  if pl_set.isparam('csd')
    pl_set.remove('csd');
  end
  if pl_set.isparam('fs')
    pl_set.remove('fs');
  end
  if pl_set.isparam('targetobj')
    pl_set.remove('targetobj');
  end
  if pl_set.isparam('UseSym')
    pl_set.remove('UseSym');
  end
  if pl_set.isparam('iunits')
    pl_set.remove('iunits');
  end
  if pl_set.isparam('ounits')
    pl_set.remove('ounits');
  end
  if pl_set.isparam('MaxIter')
    pl_set.remove('MaxIter');
  end
  if pl_set.isparam('MinOrder')
    pl_set.remove('MinOrder');
  end
  if pl_set.isparam('MaxOrder')
    pl_set.remove('MaxOrder');
  end
  if pl_set.isparam('PoleType')
    pl_set.remove('PoleType');
  end
  if pl_set.isparam('Weights')
    pl_set.remove('Weights');
  end
  if pl_set.isparam('FITTOL')
    pl_set.remove('FITTOL');
  end
  if pl_set.isparam('MSEVARTOL')
    pl_set.remove('MSEVARTOL');
  end
  if pl_set.isparam('plot')
    pl_set.remove('plot');
  end
  a.setProperties(pl_set);
  warning('on', utils.const.warnings.METHOD_NOT_FOUND);
  
  
end