--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/noisegen1D.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,413 @@
+% 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
+
+