% NOISEGEN1D generates colored noise from white noise.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: noisegen1D can work in two different modes:
% 
% ------------------------------------------------------------------------
% 1) Generates colored noise from white noise with a given spectrum. The
% function constructs a coloring filter through a fitting procedure to the
% model provided. If no model is provided an error is prompted. The colored
% noise provided has one-sided psd corresponding to the input model.
% 
% This mode correspond to the 'Default' set for the method (see the list of
% parameters).
%
% ALGORITHM:
%            1) Fit a set of partial fraction z-domain filters using
%               utils.math.psd2tf
%            2) Convert to array of MIIR filters
%            3) Filter time-series in parallel
%
% CALL:         b = noisegen1D(a, pl)
% 
% INPUT:
%               - a is a white noise time-series analysis object or a
%               vector of analysis objects
%               - pl is a plist with the input parameters
% 
% OUTPUT:
% 
%               - b Colored time-series AOs. The coloring filters used
%               are stored in the objects procinfo field under the
%               parameter 'Filt'.
% 
% ------------------------------------------------------------------------
% 
% 2) Generates noise coloring filters for given input psd models. 
% 
% This mode correspond to the 'Filter' set for the method (see the list of
% parameters).
%
% ALGORITHM:
%            1) Fit a set of partial fraction z-domain filters
%            2) Convert to array of MIIR filters
%
% CALL:         fil = noisegen1D(psd, pl)
% 
% INPUT:
%               - psd is a fsdata analysis object representing the desired
%               model for the power spectral density of the colored noise
%               - pl is a plist with the input parameters
% 
% OUTPUT:
% 
%               - fil is a filterbank parallel object which elements are
%               miir filters. Filters are initialized to avoid startup
%               transients.
%
% ------------------------------------------------------------------------
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'noisegen1D')">Parameters Description</a>
%
% VERSION:     $Id: noisegen1D.m,v 1.30 2011/04/08 08:56:15 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = noisegen1D(varargin)
  
  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end
  
  import utils.const.*
  utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
  
  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end
  
  % Collect all AOs and plists
  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
  
  % Decide on a deep copy or a modify
  bs = copy(as, nargout);
  inhists = [as.hist];
  
  % combine plists
  % define input/output combinations. Different combination are
  % 1) input tsdata and csd# into the plist, output are colored tsdata
  % 2) input fsdata, output is a coloring filter (in a matrix)
  if isempty(pl) % no model input, get model from input
    setpar = 'Filter';
  elseif isa(as(1).data,'fsdata') % get model from input
    setpar = 'Filter';
  else % get model from plist, output tsdata, back compatibility mode
    setpar = 'Default';
  end
  
  pl = parse(pl, getDefaultPlist(setpar));
  pl.getSetRandState();
  
  switch lower(setpar)
    case 'default'
      % Extract necessary parameters
      model = find(pl, 'model');
    case 'filter'
      fs     = find(pl,'fs');
      Iunits = find(pl, 'Iunits');
      Ounits = find(pl, 'Ounits');
      
      % init filter objects
      fil = filterbank.initObjectWithSize(size(bs,1),size(bs,2));
  end
  
  % start building input structure for psd2tf
  params.idtp = 1;
  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');
  params.spy = find(pl, 'Disp');
  
  % 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.usesym = 0;
  params.dterm = 0; % it is better to fit without dterm
  
  % Loop over input AOs
  for jj=1:numel(bs)
    
    switch lower(setpar)
      case 'default'
        if ~isa(bs(jj).data, 'tsdata')
          warning('!!! %s expects ao/tsdata objects. Skipping AO %s', mfilename, ao_invars{jj});
        else      
          %-------------- Colored noise from white noise

          % 1) If we have no model gives an error
          if(isempty(model))
            error('!!! Input a model for the desired PSD')
          end

          % 2) Noise Generation

          params.fs = bs(jj).fs;

          % call psd2tf
          [res, poles, dterm, mresp, rdl] = ...
            utils.math.psd2tf(model.y,[],[],[],model.x,params);
          
          % get init state
          Zi = utils.math.getinitstate(res,poles,1,'mtd','svd');

          % 3) Convert to MIIR filters
          % add yunits, taking them from plist or, if empty, from input object
          Iunits = as(jj).yunits;
          Ounits = find(pl, 'yunits');
          if eq(unit(Ounits), unit(''))
            Ounits = as(jj).yunits;
          end
          pfilts = [];
          for kk=1:numel(res)
            ft = miir(res(kk), [ 1 -poles(kk)], bs(jj).fs);
            ft.setIunits(Iunits);
            ft.setOunits(Ounits);
            ft.setHistout(Zi(kk)); 
            % build parallel bank of filters
            pfilts = [pfilts ft];
          end

          % 4) Filter data
          bs(jj).filter(pfilts);

          % 5) Output data            
          % add history
          bs(jj).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
          % name for this object
          bs(jj).setName(sprintf('noisegen1D(%s)', ao_invars{jj}));
          % Collect the filters into procinfo
          bs(jj).procinfo = plist('filter', pfilts);
      
        end
      case 'filter'
        
        params.fs = fs;

        % call psd2tf
        [res, poles, dterm, mresp, rdl] = ...
          utils.math.psd2tf(bs(jj).y,[],[],[],bs(jj).x,params);

        % get init state
        Zi = utils.math.getinitstate(res,poles,1,'mtd','svd');
        
        % build miir
        pfilts = [];
        for kk=1:numel(res)
          ft = miir(res(kk), [ 1 -poles(kk)], fs);
          ft.setIunits(Iunits);
          ft.setOunits(Ounits);
          ft.setHistout(Zi(kk)); 
          % build parallel bank of filters
          pfilts = [pfilts ft];
        end
        
        % build output filterbanks
        fil(jj) = filterbank(plist('filters',pfilts,'type','parallel'));
        fil(jj).setName(sprintf('noisegen1D(%s)',ao_invars{jj}));
        fil(jj).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
    end
  end
  
  % switch output
  switch lower(setpar)
    
    case 'default'
      % Set output
      if nargout > 1 && nargout == numel(bs)
        % List of outputs
        for ii = 1:numel(bs)
          varargout{ii} = bs.index(ii);
        end
      else
        % Single output
        varargout{1} = bs;
      end
      
    case 'filter'
      
      % Set output
      if nargout > 1 && nargout == numel(bs)
        % List of outputs
        for ii = 1:numel(fil)
          varargout{ii} = fil.index(ii);
        end
      else
        % Single output
        varargout{1} = fil;
      end
      
  end
end

%--------------------------------------------------------------------------
% Get Info Object
%--------------------------------------------------------------------------
function ii = getInfo(varargin)
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pl   = [];
  elseif nargin == 1 && ~isempty(varargin{1}) && ischar(varargin{1})
    sets{1} = varargin{1};
    pl = getDefaultPlist(sets{1});
  else
    sets = SETS();
    % get plists
    pl(size(sets)) = plist;
    for kk = 1:numel(sets)
      pl(kk) =  getDefaultPlist(sets{kk});
    end
  end
  % Build info object
  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: noisegen1D.m,v 1.30 2011/04/08 08:56:15 hewitson Exp $', sets, pl);

end


%--------------------------------------------------------------------------
% Defintion of Sets
%--------------------------------------------------------------------------

function out = SETS()
  out = {...
    'Default', ...
    'Filter'   ...
    };
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------
function plout = getDefaultPlist(set)
  persistent pl;
  persistent lastset;
  if exist('pl', 'var')==0 || isempty(pl) || ~strcmp(lastset, set)
    pl = buildplist(set);
    lastset = set;
  end
  plout = pl;
end

function pl = buildplist(set)
  
  pl = plist();
  
  switch lower(set)
    case 'default'
      % Yunits
      p = param({'yunits',['Unit on Y axis. <br>' ...
        'If left empty, it will take the y-units from the input object']},  paramValue.STRING_VALUE(''));
      pl.append(p);

      % Model
      p = param({'model', 'A frequency-series AO describing the model psd'}, paramValue.EMPTY_DOUBLE);
      pl.append(p);
      
    case 'filter'
      % Fs
      p = param({'fs','The sampling frequency to design for.'}, paramValue.DOUBLE_VALUE(1));
      pl.append(p);
      
      % Iunits
      p = param({'iunits','The input units of the filter.'}, paramValue.EMPTY_STRING);
      pl.append(p);
      
      % Ounits
      p = param({'ounits','The output units of the filter.'}, paramValue.EMPTY_STRING);
      pl.append(p);
  end
  
  % MaxIter
  p = param({'MaxIter', 'Maximum number of iterations in fit routine.'}, paramValue.DOUBLE_VALUE(30));
  pl.append(p);
  
  % PoleType
  p = param({'PoleType', ['Choose the pole type for fitting:<ol>'...
                      '<li>use real starting poles</li>' ...
                      '<li>generates complex conjugate poles of the<br>'...
                      'type <tt>a.*exp(theta*pi*j)</tt><br>'...
                      'with <tt>theta = linspace(0,pi,N/2+1)</tt></li>'...
                      '<li>generates complex conjugate poles of the type<br>'...
                      '<tt>a.*exp(theta*pi*j)</tt><br>'...
                      'with <tt>theta = linspace(0,pi,N/2+2)</tt></li></ol>']}, {3, {1,2,3}, paramValue.SINGLE});
  pl.append(p);
  
  % MinOrder
  p = param({'MinOrder', 'Minimum order to fit with.'}, paramValue.DOUBLE_VALUE(2));
  pl.append(p);
  
  % MaxOrder
  p = param({'MaxOrder', 'Maximum order to fit with.'}, paramValue.DOUBLE_VALUE(25));
  pl.append(p);
  
  % Weights
  p = param({'Weights', ['Choose weighting for the fit:<ol>'...
                         '<li>equal weights for each point</li>'...
                         '<li>weight with <tt>1/abs(model)</tt></li>'...
                         '<li>weight with <tt>1/abs(model).^2</tt></li>'...
                         '<li>weight with inverse of the square mean spread<br>'...
                         'of the model</li></ol>']}, paramValue.DOUBLE_VALUE(3));
  pl.append(p);
  
  % Plot
  p = param({'Plot', 'Plot results of each fitting step.'}, paramValue.FALSE_TRUE);  
  pl.append(p);
  
  % Disp
  p = param({'Disp', 'Display the progress of the fitting iteration.'}, paramValue.FALSE_TRUE);  
  pl.append(p);
  
  % MSEVARTOL
  p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the<br>'...
                           'realtive variation of the mean squared error is<br>'...
                           'smaller than the value specified. This<br>'...
                           'option is useful for finding the minimum of Chi-squared.']}, ...
                           paramValue.DOUBLE_VALUE(1e-2));
  pl.append(p);
  
  % FITTOL
  p = param({'FITTOL', ['Mean Squared Error Value - Check if the mean<br>'...
                        'squared error value is lower than the value<br>'...
                        'specified.']}, paramValue.DOUBLE_VALUE(1e-2));
  pl.append(p);
    
end