Mercurial > hg > ltpda
view m-toolbox/classes/@ao/crbound.m @ 5:5a49956df427 database-connection-manager
LTPDAPreferences panel for new LTPDADatabaseConnectionManager
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
% 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