--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/crbound.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,410 @@
+% CRBOUND computes the inverse of the Fisher Matrix
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: CRBOUND computes the Cramer-Rao lowe bound for parametric
+%              model given the input signals and the noise. The method
+%              accepts 2D (2 input/2 output) models and 1D models
+%              (1 input/1 output) and (2 input/1 output)
+%
+% CALL:        bs = crbound(in,noise,pl)
+%
+% INPUTS:      in      - analysis objects with input signals to the system
+%                        (x1) if 1 input 
+%                        (x2) if 2 input 
+%              noise   - analysis objects with measured noise 
+%                        (x1) if 1 input: S11 
+%                        (x4) if 2 input: (S11,S12,S21,S22) 
+%              model   - symbolic model (smodel) containing the evaluated
+%                        transfer function models
+%                        (x1) if 1 input/1 output
+%                        (x2) if 2 input/1 output
+%                        (x4) if 2 input/2 output
+%              pl      - parameter list
+%
+% OUTPUTS:     bs   - covariance matrix AO
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'crbound')">Parameters Description</a>
+%
+% VERSION:    $Id: crbound.m,v 1.13 2011/04/08 08:56:17 hewitson Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = crbound(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);
+  
+  % Method can not be used as a modifier
+  if nargout == 0
+    error('### crbound cannot be used as a modifier. Please give an output variable.');
+  end
+  
+  % Collect input variable names
+  in_names = cell(size(varargin));
+  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+  
+  % Collect all AOs smodels and plists
+  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+  % Combine plists
+  pl = parse(pl, getDefaultPlist);
+  
+  % get params
+  params = find(pl,'params');
+  numparams = find(pl,'paramsValues');
+  mdl2 = find(pl,'models');
+  bmdl = find(pl,'built-in');
+  f1 = find(pl,'f1');
+  f2 = find(pl,'f2');
+  freqs = find(pl,'frequencies');
+  pseudoinv = find(pl,'pinv');
+  tol = find(pl,'tol');
+  %   ninputs = find(pl,'Ninputs');
+  
+  
+  % Decide on a deep copy or a modify
+  bs = copy(as, nargout);
+  if ~isempty(mdl2)
+    mdl = copy(mdl2, nargout);
+  end
+  
+  if numel(bs) == 2
+    % implementing here for 1D
+    in(1) = bs(1);
+    n(1) = bs(2);
+    
+    
+    
+    % fft
+    i1 = fft(in(1));
+    Nsamples = length(i1.x);
+    if ~isempty(f1) &&  ~isempty(f2)
+      i1 = split(i1,plist('frequencies',[f1 f2]));
+    elseif ~isempty(freqs)
+      i1 = split(i1,plist('frequencies',freqs));
+      i1 = join(i1);
+    end
+    
+    % get frequency vector
+    f = i1.x;
+    
+    if ~isempty(mdl)
+      % get model
+      h(1) = mdl;
+      h(1).setXvals(f);
+    else
+      error('### One model should be introduced, please check the help.')
+    end
+    
+  elseif numel(bs) == 3
+    if numel(mdl) == 2
+      in(1) = bs(1);
+      in(2) = bs(2);
+      n(1) = bs(3);
+      
+      % fft
+      i1 = fft(in(1));
+      Nsamples = length(i1.x);
+      i2 = fft(in(2));
+      
+      % Get rid of fft f =0, reduce frequency range if needed
+      if ~isempty(f1) &&  ~isempty(f2)
+        i1 = split(i1,plist('frequencies',[f1 f2]));
+        i2 = split(i2,plist('frequencies',[f1 f2]));
+      elseif ~isempty(freqs)
+        i1 = split(i1,plist('frequencies',freqs));
+        i1 = join(i1);
+        i2 = split(i2,plist('frequencies',freqs));
+        i2 = join(i2);
+      end
+      
+      % get frequency vector
+      f = i1.x;
+      
+      if ~isempty(mdl)
+        % compute built-in matrix
+        h(1) = mdl(1);
+        h(1).setXvals(f);
+        h(2) = mdl(2);
+        h(2).setXvals(f);
+      else
+        error('### One model should be introduced, please check the help.')
+      end
+    elseif numel(mdl) == 1
+      error('Check the model or use crbound with the matrix class')
+    end
+    % Get input objects
+  elseif numel(bs) == 4
+    in(1) = bs(1);
+    in(2) = bs(2);
+    n(1) = bs(3);
+    n(2) = bs(4);
+    
+    % fft
+    i1 = fft(in(1));
+    Nsamples = length(i1.x);
+    i2 = fft(in(2));
+    
+    % Get rid of fft f =0, reduce frequency range if needed
+    if ~isempty(f1) &&  ~isempty(f2)
+      i1 = split(i1,plist('frequencies',[f1 f2]));
+      i2 = split(i2,plist('frequencies',[f1 f2]));
+    elseif ~isempty(freqs)
+      i1 = split(i1,plist('frequencies',freqs));
+      i1 = join(i1);
+      i2 = split(i2,plist('frequencies',freqs));
+      i2 = join(i2);
+    end
+    
+    % get frequency vector
+    f = i1.x;
+    
+    if ~isempty(bmdl)
+      % compute built-in smodels
+      for i = 1:4
+        if strcmp(bmdl{i},'0');
+          h(i) = smodel('0');
+          h(i).setXvar('f');
+          h(i).setXvals(f);
+        else
+          h(i) = smodel(plist('built-in',bmdl{i},'f',f));
+          % set all params to all models. It is not true but harmless
+          for k = 1:numel(params)
+            vecparams(k) = {numparams(k)*ones(size(f))};
+          end
+          h(i).setParams(params,vecparams);
+        end
+      end
+    elseif ~isempty(mdl)
+      % compute built-in matrix
+      h(1) = mdl{1};
+      h(1).setXvals(f);
+      h(2) = mdl{2};
+      h(2).setXvals(f);
+      h(3) = mdl{3};
+      h(3).setXvals(f);
+      h(4) = mdl{4};
+      h(4).setXvals(f);
+    end
+    
+  else
+    error('### The number of input in crbound objects is not correct, please check the help.')
+  end
+  
+  % Set params
+  for i = 1:numel(h)
+    h(i).setParams(params,numparams);
+  end
+  
+  if numel(bs) == 2
+    % Compute psd
+    n1  = psd(n(1), pl);
+    
+    % interpolate to fft frequencies
+    S11 = interp(n1,plist('vertices',f));
+    
+    for i = 1:length(params)
+      utils.helper.msg(msg.IMPORTANT, sprintf('computing symbolic differentiation with respect %s',params{i}), mfilename('class'), mfilename);
+      % differentiate symbolically
+      dH11 = diff(h(1),params{i});
+      % evaluate
+      d11(i) = eval(dH11);
+    end
+    
+    % get some parameters used below
+    fs = S11.fs;
+%     N = len(S11);
+    Navs = find(pl,'Navs');
+    
+    % scaling of PSD
+    % PSD = 2/(N*fs) * FFT *conj(FFT)
+    C11 = Nsamples*fs/2.*S11.y;
+
+    InvS11 = 1./C11;
+
+    % compute Fisher Matrix
+    for i =1:length(params)
+      for j =1:length(params)
+        v1v1 = conj(d11(i).y.*i1.y).*(d11(j).y.*i1.y);
+
+        FisMat(i,j) = sum(real(InvS11.*v1v1));
+      end
+    end
+    
+  elseif numel(bs) == 3
+  
+    if numel(mdl) == 2
+      % Compute psd
+      n1  = psd(n(1), pl);
+      
+      % interpolate to fft frequencies
+      S11 = interp(n1,plist('vertices',f));
+      
+      for i = 1:length(params)
+        utils.helper.msg(msg.IMPORTANT, sprintf('computing symbolic differentiation with respect %s',params{i}), mfilename('class'), mfilename);
+        % differentiate symbolically
+        dH11 = diff(h(1),params{i});
+        dH12 = diff(h(2),params{i});
+        
+        % evaluate
+        d11(i) = eval(dH11);
+        d12(i) = eval(dH12);    
+      end
+      
+      % get some parameters used below
+      fs = S11.fs;
+%       N = len(S11);
+      Navs = find(pl,'Navs');
+      
+      % scaling of PSD
+      % PSD = 2/(N*fs) * FFT *conj(FFT)
+      C11 = Nsamples*fs/2.*S11.y;
+      
+      % compute elements of inverse cross-spectrum matrix
+      InvS11 = 1./C11;
+      
+      % compute Fisher Matrix
+      for i =1:length(params)
+        for j =1:length(params)
+          v1v1 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
+          
+          FisMat(i,j) = sum(real(InvS11.*v1v1));
+        end
+      end
+      
+    elseif numel(mdl) == 1
+      error('Please check model sizes')
+    end
+    
+    
+  elseif numel(bs) == 4
+    
+    % Compute psd
+    n1  = psd(n(1), pl);
+    n2  = psd(n(2), pl);
+    n12 = cpsd(n(1),n(2), pl);
+    
+    % interpolate to fft frequencies
+    S11 = interp(n1,plist('vertices',f));
+    S12 = interp(n12,plist('vertices',f));
+    S22 = interp(n2,plist('vertices',f));
+    S21 = conj(S12);
+    
+    for i = 1:length(params)
+      utils.helper.msg(msg.IMPORTANT, sprintf('computing symbolic differentiation with respect %s',params{i}), mfilename('class'), mfilename);
+      % differentiate symbolically
+      dH11 = diff(h(1),params{i});
+      dH12 = diff(h(2),params{i});
+      dH21 = diff(h(3),params{i});
+      dH22 = diff(h(4),params{i});
+      % evaluate
+      d11(i) = eval(dH11);
+      d12(i) = eval(dH12);
+      d21(i) = eval(dH21);
+      d22(i) = eval(dH22);
+    end
+    
+    % get some parameters used below
+    fs = S11.fs;
+%     N = len(S11);
+    Navs = find(pl,'Navs');
+    
+    % scaling of PSD
+    % PSD = 2/(N*fs) * FFT *conj(FFT)
+    C11 = Nsamples*fs/2.*S11.y;
+    C22 = Nsamples*fs/2.*S22.y;
+    C12 = Nsamples*fs/2.*S12.y;
+    C21 = Nsamples*fs/2.*S21.y;
+    
+    % compute elements of inverse cross-spectrum matrix
+    InvS11 = (C22./(C11.*C22 - C12.*C21));
+    InvS22 = (C11./(C11.*C22 - C12.*C21));
+    InvS12 = (C21./(C11.*C22 - C12.*C21));
+    InvS21 = (C12./(C11.*C22 - C12.*C21));
+    
+    % compute Fisher Matrix
+    for i =1:length(params)
+      for j =1:length(params)
+        
+        v1v1 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
+        v2v2 = conj(d21(i).y.*i1.y + d22(i).y.*i2.y).*(d21(j).y.*i1.y + d22(j).y.*i2.y);
+        v1v2 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d21(j).y.*i1.y + d22(j).y.*i2.y);
+        v2v1 = conj(d21(i).y.*i1.y + d22(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
+        
+        FisMat(i,j) = sum(real(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1));
+      end
+    end
+    
+  end
+  
+  % inverse is the optimal covariance matrix
+  if pseudoinv && isempty(tol)
+    cov = pinv(FisMat);
+  elseif pseudoinv
+    cov = pinv(FisMat,tol);
+  else
+    cov = FisMat\eye(size(FisMat));
+  end
+  
+  % create AO
+  out = ao(cov);
+  % Fisher Matrix in the procinfo
+  out.setProcinfo(plist('FisMat',FisMat));
+  
+  varargout{1} = out;
+end
+
+
+%--------------------------------------------------------------------------
+% Get Info Object
+%--------------------------------------------------------------------------
+function ii = getInfo(varargin)
+  if nargin == 1 && strcmpi(varargin{1}, 'None')
+    sets = {};
+    pls  = [];
+  else
+    sets = {'Default'};
+    pls  = getDefaultPlist;
+  end
+  % Build info object
+  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: crbound.m,v 1.13 2011/04/08 08:56:17 hewitson Exp $', sets, pls);
+end
+
+%--------------------------------------------------------------------------
+% Get Default Plist
+%--------------------------------------------------------------------------
+
+function plout = getDefaultPlist()
+  persistent pl;  
+  if exist('pl', 'var')==0 || isempty(pl)
+    pl = buildplist();
+  end
+  plout = pl;  
+end
+
+function pl = buildplist()
+  pl = plist.WELCH_PLIST;
+  pset(pl,'Navs',1)
+  
+  p = plist({'params', 'Parameters of the model'}, paramValue.EMPTY_STRING);
+  pl.append(p);
+  
+  p = plist({'models','Symbolic models of the system'}, paramValue.EMPTY_STRING);
+  pl.append(p);
+  
+  p = plist({'frequencies','Array of start/sop frequencies where the analysis is performed'}, paramValue.EMPTY_STRING);
+  pl.append(p);
+  
+  p = plist({'pinv','Use the Penrose-Moore pseudoinverse'}, paramValue.TRUE_FALSE);
+  pl.append(p);
+  
+  p = plist({'tol','Tolerance for the Penrose-Moore pseudoinverse'}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+end