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