Mercurial > hg > ltpda

annotate m-toolbox/classes/@matrix/dispersion.m @ 42:f90d4f666cc7 database-connection-manager

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleanup
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Mon, 05 Dec 2011 18:04:34 +0100
parents f0afece42f48
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
1 % DISPERSION computes the dispersion function
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
3 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
4 % DESCRIPTION: DISPERSION computes the dispersion function
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
5 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
6 % CALL: bs = dipersion(in,pl)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
7 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
8 % INPUTS: in - matrix objects with input signals to the system
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
9 % model - symbolic models containing the transfer function model
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
10 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
11 % pl - parameter list
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
12 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
13 % OUTPUTS: bs - dispersion function AO
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
14 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
15 % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'dispersion')">Parameters Description</a>
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
16 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
17 % VERSION: $Id: dispersion.m,v 1.1 2011/06/22 09:54:19 miquel Exp $
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
18 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
19 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
20
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
21 function varargout = dispersion(varargin)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
22
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
23 % Check if this is a call for parameters
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
24 if utils.helper.isinfocall(varargin{:})
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
25 varargout{1} = getInfo(varargin{3});
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
26 return
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
27 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
28
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
29 import utils.const.*
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
30 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
31
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
32 % Method can not be used as a modifier
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
33 if nargout == 0
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
34 error('### crb cannot be used as a modifier. Please give an output variable.');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
35 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
36
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
37 % Collect input variable names
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
38 in_names = cell(size(varargin));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
39 for ii = 1:nargin,in_names{ii} = inputname(ii);end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
40
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
41 % Collect all AOs smodels and plists
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
42 [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
43 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
44
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
45 % Combine plists
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
46 pl = parse(pl, getDefaultPlist);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
47
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
48 % get params
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
49 params = find(pl,'FitParams');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
50 numparams = find(pl,'paramsValues');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
51 mmdl = find(pl,'model');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
52 channel = find(pl,'channel');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
53 mtxns = find(pl,'noise');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
54 outModel = find(pl,'outModel');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
55 bmdl = find(pl,'built-in');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
56 f1 = find(pl,'f1');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
57 f2 = find(pl,'f2');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
58 pseudoinv = find(pl,'pinv');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
59 tol = find(pl,'tol');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
60 outNames = find(pl,'outNames');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
61 inNames = find(pl,'inNames');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
62
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
63 % Decide on a deep copy or a modify
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
64 fin = copy(mtxs, nargout);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
65 n = copy(mtxns, nargout);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
66 mdl = copy(mmdl,1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
67
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
68 % Get number of experiments
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
69 nexp = numel(fin);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
70
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
71 % fft
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
72 % fin = fft(in);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
73
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
74 % N should get before spliting, in order to convert correctly from psd to
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
75 % fft
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
76 N = length(fin(1).getObjectAtIndex(1).x);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
77
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
78 % Get rid of fft f =0, reduce frequency range if needed
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
79 if ~isempty(f1) && ~isempty(f2)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
80 fin = split(fin,plist('frequencies',[f1 f2]));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
81 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
82
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
83 FMall = zeros(numel(params),numel(params));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
84 % loop over experiments
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
85 for k = 1:nexp
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
86
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
87 utils.helper.msg(msg.IMPORTANT, sprintf('Analysis of experiment #%d',k), mfilename('class'), mfilename);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
88
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
89 % if (((numel(n(1).objs)) == 2) && (numel(fin(1).objs) == 2))
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
90
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
91 % use signal fft to get frequency vector. Take into account signal
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
92 % could be empty or set to zero
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
93 % 1st channel
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
94 if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
95 i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
96 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
97 i1 = fin(k).getObjectAtIndex(1,1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
98 freqs = i1.x;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
99 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
100 % 2nd channel
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
101 if all(fin(k).getObjectAtIndex(2,1).y == 0) || isempty(fin(k).getObjectAtIndex(2,1).y)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
102 i2 = ao(plist('type','fsdata','xvals',0,'yvals',0));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
103 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
104 i2 = fin(k).getObjectAtIndex(2,1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
105 freqs = i2.x;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
106 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
107
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
108 % Compute psd
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
109 n1 = lpsd(n(k).getObjectAtIndex(1,1));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
110 n2 = lpsd(n(k).getObjectAtIndex(2,1));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
111 n12 = lcpsd(n(k).getObjectAtIndex(1,1),n(k).getObjectAtIndex(2,1));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
112
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
113 % interpolate to given frequencies
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
114 % noise
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
115 S11 = interp(n1,plist('vertices',freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
116 S12 = interp(n12,plist('vertices',freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
117 S22 = interp(n2,plist('vertices',freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
118 S21 = conj(S12);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
119
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
120 % get some parameters used below
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
121 fs = S11.fs;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
122
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
123
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
124 % if (isempty(outModel))
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
125
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
126 % if ~isempty(bmdl)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
127 % compute built-in smodels
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
128 % for i = 1:4
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
129 % if strcmp(bmdl{i},'0');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
130 % h(i) = smodel('0');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
131 % h(i).setXvar('f');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
132 % h(i).setXvals(freqs);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
133 % h(i).setParams(params,numparams);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
134 % else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
135 % h(i) = smodel(plist('built-in',bmdl{i},'f',freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
136 % % set all params to all models. It is not true but harmless
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
137 % for ii = 1:numel(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
138 % vecparams(ii) = {numparams(ii)*ones(size(freqs))};
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
139 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
140 % h(i).setParams(params,vecparams);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
141 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
142 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
143 % elseif ~isempty(mdl) && all(strcmp(class(mdl),'matrix'))
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
144
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
145 for i = 1:numel(mdl.objs)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
146 % set Xvals
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
147 h(i) = mdl.getObjectAtIndex(i).setXvals(freqs);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
148 % set alias
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
149 h(i).assignalias(mdl.objs(i),plist('xvals',freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
150 % set paramaters
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
151 h(i).setParams(params,numparams);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
152 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
153 % differentiate and eval
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
154 for i = 1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
155 utils.helper.msg(msg.IMPORTANT, sprintf('computing symbolic differentiation with respect %s',params{i}), mfilename('class'), mfilename);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
156 % differentiate symbolically
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
157 dH11 = diff(h(1),params{i});
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
158 dH12 = diff(h(3),params{i}); % taking into account matrix index convention h(2) > H(2,1)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
159 dH21 = diff(h(2),params{i});
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
160 dH22 = diff(h(4),params{i});
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
161 % evaluate
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
162 d11(i) = eval(dH11,plist('output type','fsdata','output x',freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
163 d12(i) = eval(dH12,plist('output type','fsdata','output x',freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
164 d21(i) = eval(dH21,plist('output type','fsdata','output x',freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
165 d22(i) = eval(dH22,plist('output type','fsdata','output x',freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
166 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
167
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
168 % elseif ~isempty(mdl) && all(strcmp(class(mdl),'ssm'))
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
169 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
170 % meval = copy(mdl,1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
171 % % set parameter values
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
172 % meval.doSetParameters(params, numparams);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
173 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
174 % % make numeric
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
175 % % meval.doSubsParameters(params, true);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
176 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
177 % % get the differentiation step
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
178 % step = find(pl,'diffStep');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
179 % if isempty(step)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
180 % error('### Please input a step for the numerical differentiation')
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
181 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
182 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
183 % % differentiate and eval
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
184 % for i = 1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
185 % utils.helper.msg(msg.IMPORTANT, sprintf('computing numerical differentiation with respect %s, Step:%4.2d ',params{i},step(i)), mfilename('class'), mfilename);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
186 % % differentiate numerically
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
187 % dH = meval.parameterDiff(plist('names', params(i),'values',step(i)));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
188 % % create plist with correct outNames (since parameterDiff change them)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
189 % out1 = strrep(outNames{1},'.', sprintf('_DIFF_%s.',params{i})); % 2x2 case
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
190 % out2 =strrep(outNames{2},'.', sprintf('_DIFF_%s.',params{i}));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
191 % spl = plist('set', 'for bode', ...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
192 % 'outputs', {out1,out2}, ...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
193 % 'inputs', inNames, ...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
194 % 'reorganize', true,...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
195 % 'f', freqs);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
196 % % do bode
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
197 % d = bode(dH, spl);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
198 % % assign according matlab's matrix notation: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
199 % d11(i) = d(1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
200 % d21(i) = d(2);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
201 % d12(i) = d(3);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
202 % d22(i) = d(4);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
203 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
204 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
205 % else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
206 % error('### please introduce models for the transfer functions')
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
207 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
208
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
209 % elseif (~isempty(outModel))
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
210 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
211 % % if(~isempty(outModel))
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
212 % for lll=1:size(outModel,1)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
213 % for kkk=1:size(outModel,2)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
214 % outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[f1 f2]));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
215 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
216 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
217 % %end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
218 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
219 % % Avoid numerical differentiation (faster for the magnetic case)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
220 % Param{1} = [1 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
221 % 0 0;];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
222 % Param{2} = [0 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
223 % 0 1;];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
224 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
225 % for pp = 1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
226 % for ll = 1:size(outModel,1)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
227 % for kk = 1:size(Param{pp},2)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
228 % % index convention: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
229 % tmp = 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
230 % for innerIndex = 1:size(outModel,2)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
231 % tmp = tmp + outModel(ll,innerIndex).y * Param{pp}(innerIndex,kk);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
232 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
233 % h{pp}(:,(kk-1)*size(outModel,1) + ll) = tmp;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
234 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
235 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
236 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
237 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
238 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
239 % for i = 1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
240 % d11(i).y = h{i}(:,1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
241 % d21(i).y = h{i}(:,2);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
242 % d12(i).y = h{i}(:,3);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
243 % d22(i).y = h{i}(:,4);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
244 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
245 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
246 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
247
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
248 % scaling of PSD
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
249 % PSD = 2/(N*fs) * FFT *conj(FFT)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
250 C11 = N*fs/2.*S11.y;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
251 C22 = N*fs/2.*S22.y;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
252 C12 = N*fs/2.*S12.y;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
253 C21 = N*fs/2.*S21.y;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
254
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
255 % compute elements of inverse cross-spectrum matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
256 InvS11 = (C22./(C11.*C22 - C12.*C21));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
257 InvS22 = (C11./(C11.*C22 - C12.*C21));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
258 InvS12 = (C21./(C11.*C22 - C12.*C21));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
259 InvS21 = (C12./(C11.*C22 - C12.*C21));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
260
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
261
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
262 % compute Fisher Matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
263 for i =1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
264 for j =1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
265
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
266 v1v1 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
267 v2v2 = conj(d21(i).y.*i1.y + d22(i).y.*i2.y).*(d21(j).y.*i1.y + d22(j).y.*i2.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
268 v1v2 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d21(j).y.*i1.y + d22(j).y.*i2.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
269 v2v1 = conj(d21(i).y.*i1.y + d22(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
270
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
271 FisMat(i,j) = sum(real(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
272 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
273 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
274 utils.helper.msg(msg.IMPORTANT, sprintf('rank(FisMat) = %d', rank(FisMat)));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
275
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
276 for kk = 1:numel(freqs)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
277 % create input signal with power at single freq.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
278 % depending on input channel
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
279 p = zeros(1,numel(freqs));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
280 if channel == 1
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
281 p(kk) = sum(i1.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
282 % create aos
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
283 i1single = ao(plist('Xvals',freqs,'Yvals',p));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
284 i2single = ao(plist('Xvals',freqs,'Yvals',zeros(1,numel(freqs))));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
285 elseif channel == 2
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
286 p(kk) = sum(i2.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
287 % create aos
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
288 i1single = ao(plist('Xvals',freqs,'Yvals',zeros(1,numel(freqs))));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
289 i2single = ao(plist('Xvals',freqs,'Yvals',p));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
290 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
291 error('### wrong channel')
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
292 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
293
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
294 % compute Fisher Matrix for single frequencies
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
295 for i =1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
296 for j =1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
297
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
298 v1v1 = conj(d11(i).y.*i1single.y + d12(i).y.*i2single.y).*(d11(j).y.*i1single.y + d12(j).y.*i2single.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
299 v2v2 = conj(d21(i).y.*i1single.y + d22(i).y.*i2single.y).*(d21(j).y.*i1single.y + d22(j).y.*i2single.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
300 v1v2 = conj(d11(i).y.*i1single.y + d12(i).y.*i2single.y).*(d21(j).y.*i1single.y + d22(j).y.*i2single.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
301 v2v1 = conj(d21(i).y.*i1single.y + d22(i).y.*i2single.y).*(d11(j).y.*i1single.y + d12(j).y.*i2single.y);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
302
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
303 FisMatsingle(i,j) = sum(real(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
304 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
305 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
306 d(kk) = trace(FisMat\FisMatsingle)/numel(i1.x); % had to divide for num. freqs
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
307 % d(kk) = trace(pinv(FisMat)*FisMatsingle)/numel(i1.x); % had to divide for num. freqs
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
308
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
309 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
310
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
311 % % store Fisher Matrix for this run
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
312 % FM{k} = FisMat;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
313 % % adding up
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
314 % FMall = FMall + FisMat;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
315
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
316 % elseif ((numel(n(1).objs) == 3) && (numel(in.objs) == 4) && ~isempty(outModel))
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
317 % % this is only valid for the magnetic model, where we have 4 inputs
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
318 % % (corresponding to the 4 conformator waveforms) and 3 outputs
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
319 % % (corresponding to IFO.x12, IFO.eta1 and IFO.phi1). And there is a
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
320 % % contribution of an outModel converting the conformator waveforms
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
321 % % into forces and torques.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
322 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
323 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
324 % % For other cases not implemented yet.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
325 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
326 % % use signal fft to get frequency vector. Take into account signal
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
327 % % could be empty or set to zero
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
328 % % 1st channel
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
329 % freqs = fin.getObjectAtIndex(1,1).x;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
330 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
331 % for ii = 1:numel(n.objs)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
332 % for jj = ii:numel(n.objs)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
333 % % Compute psd
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
334 % if (ii==jj)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
335 % spec(ii,jj) = psd(n(k).getObjectAtIndex(ii), pl);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
336 % S2(ii,jj) = interp(spec(ii,jj),plist('vertices',freqs,'method','linear'));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
337 % else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
338 % spec(ii,jj) = cpsd(n(k).getObjectAtIndex(ii),n(k).getObjectAtIndex(jj),pl);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
339 % S2(ii,jj) = interp(spec(ii,jj),plist('vertices',freqs,'method','linear'));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
340 % S2(jj,ii) = conj(S2(ii,jj));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
341 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
342 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
343 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
344 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
345 % S = matrix(S2,plist('shape',[numel(n.objs) numel(n.objs)]));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
346 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
347 % % get some parameters used below
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
348 % fs = S.getObjectAtIndex(1,1).fs;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
349 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
350 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
351 % if(~isempty(outModel))
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
352 % for lll=1:size(outModel,1)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
353 % for kkk=1:size(outModel,2)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
354 % outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[f1 f2]));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
355 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
356 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
357 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
358 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
359 % % Avoid numerical differentiation (faster for the magnetic case)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
360 % Param{1} = [ 1 0 0 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
361 % 0 0 0 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
362 % 0 0 0 0;];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
363 % Param{2} = [ 0 1 0 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
364 % 0 0 0 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
365 % 0 0 0 0;];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
366 % Param{3} = [ 0 0 0 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
367 % 0 0 1 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
368 % 0 0 0 0;];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
369 % Param{4} = [ 0 0 0 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
370 % 0 0 0 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
371 % 0 0 0 1;];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
372 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
373 % % scaling of PSD
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
374 % % PSD = 2/(N*fs) * FFT *conj(FFT)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
375 % for j = 1: numel(S.objs)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
376 % % spectra to variance
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
377 % C(:,j) = (N*fs/2)*S.objs(j).data.getY;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
378 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
379 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
380 % detm = (C(:,1).*C(:,5).*C(:,9) + ...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
381 % C(:,2).*C(:,6).*C(:,7) + ...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
382 % C(:,3).*C(:,4).*C(:,8) -...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
383 % C(:,7).*C(:,5).*C(:,3) -...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
384 % C(:,8).*C(:,6).*C(:,1) -...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
385 % C(:,9).*C(:,4).*C(:,2));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
386 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
387 % InvS11 = (C(:,5).*C(:,9) - C(:,8).*C(:,6))./detm;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
388 % InvS12 = -(C(:,4).*C(:,9) - C(:,7).*C(:,6))./detm;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
389 % InvS13 = (C(:,4).*C(:,8) - C(:,7).*C(:,5))./detm;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
390 % InvS21 = -(C(:,2).*C(:,9) - C(:,8).*C(:,3))./detm;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
391 % InvS22 = (C(:,1).*C(:,9) - C(:,7).*C(:,3))./detm;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
392 % InvS23 = -(C(:,1).*C(:,8) - C(:,7).*C(:,2))./detm;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
393 % InvS31 = (C(:,2).*C(:,6) - C(:,5).*C(:,3))./detm;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
394 % InvS32 = -(C(:,1).*C(:,6) - C(:,4).*C(:,3))./detm;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
395 % InvS33 = (C(:,1).*C(:,5) - C(:,4).*C(:,2))./detm;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
396 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
397 % for pp = 1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
398 % for ll = 1:size(outModel,1)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
399 % for kk = 1:size(Param{pp},2)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
400 % % index convention: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
401 % tmp = 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
402 % for innerIndex = 1:size(outModel,2)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
403 % tmp = tmp + outModel(ll,innerIndex).y * Param{pp}(innerIndex,kk);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
404 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
405 % h{pp}(:,(kk-1)*size(outModel,1) + ll) = tmp;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
406 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
407 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
408 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
409 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
410 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
411 % for kk = 1:numel(in.objs)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
412 % inV(:,kk) = fin.objs(kk).data.getY;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
413 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
414 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
415 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
416 % % compute Fisher Matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
417 % for i =1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
418 % for j =1:length(params)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
419 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
420 % for ll = 1:size(outModel,1)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
421 % tmp = 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
422 % for kk = 1:size(Param{1},2)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
423 % tmp = tmp + h{i}(:,(kk-1)*size(outModel,1) + ll).*inV(:,kk);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
424 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
425 % v{i}(:,ll) = tmp;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
426 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
427 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
428 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
429 % for ll = 1:size(outModel,1)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
430 % tmp = 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
431 % for kk = 1:size(Param{1},2)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
432 % tmp = tmp + h{j}(:,(kk-1)*size(outModel,1) + ll).*inV(:,kk);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
433 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
434 % v{j}(:,ll) = tmp;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
435 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
436 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
437 % v1v1 = conj(v{i}(:,1)).*v{j}(:,1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
438 % v1v2 = conj(v{i}(:,1)).*v{j}(:,2);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
439 % v1v3 = conj(v{i}(:,1)).*v{j}(:,3);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
440 % v2v1 = conj(v{i}(:,2)).*v{j}(:,1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
441 % v2v2 = conj(v{i}(:,2)).*v{j}(:,2);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
442 % v2v3 = conj(v{i}(:,2)).*v{j}(:,3);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
443 % v3v1 = conj(v{i}(:,3)).*v{j}(:,1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
444 % v3v2 = conj(v{i}(:,3)).*v{j}(:,2);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
445 % v3v3 = conj(v{i}(:,3)).*v{j}(:,3);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
446 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
447 % FisMat(i,j) = sum(real(InvS11.*v1v1 +...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
448 % InvS12.*v1v2 +...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
449 % InvS13.*v1v3 +...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
450 % InvS21.*v2v1 +...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
451 % InvS22.*v2v2 +...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
452 % InvS23.*v2v3 +...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
453 % InvS31.*v3v1 +...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
454 % InvS32.*v3v2 +...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
455 % InvS33.*v3v3));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
456 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
457 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
458 % % store Fisher Matrix for this run
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
459 % FM{k} = FisMat;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
460 % % adding up
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
461 % FMall = FMall + FisMat;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
462 % else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
463 % error('Implemented cases: 2 inputs / 2outputs (TN3045 analysis), and 4 inputs / 3 outpus (magnetic complete analysis model. Other cases have not been implemented yet. Sorry for the inconvenience)');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
464 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
465
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
466
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
467 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
468
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
469 % inverse is the optimal covariance matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
470 % if pseudoinv && isempty(tol)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
471 % cov = pinv(FMall);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
472 % elseif pseudoinv
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
473 % cov = pinv(FMall,tol);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
474 % else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
475 % cov = FMall\eye(size(FMall));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
476 % end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
477
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
478
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
479 % create AO
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
480 out = ao(plist('Xvals',freqs,'Yvals',d,'type','fsdata','fs',fs,'name',''));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
481 % spectrum in the procinfo
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
482 if channel == 1
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
483 out.setProcinfo(plist('S11',S11));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
484 elseif channel == 2
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
485 out.setProcinfo(plist('S22',S22));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
486 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
487 % Fisher Matrix in the procinfo
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
488 out.setProcinfo(plist('FisMat',FisMat));
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
489
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
490 varargout{1} = out;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
491 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
492
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
493 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
494
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
495
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
496 %--------------------------------------------------------------------------
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
497 % Get Info Object
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
498 %--------------------------------------------------------------------------
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
499 function ii = getInfo(varargin)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
500 if nargin == 1 && strcmpi(varargin{1}, 'None')
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
501 sets = {};
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
502 pls = [];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
503 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
504 sets = {'Default'};
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
505 pls = getDefaultPlist;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
506 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
507 % Build info object
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
508 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: dispersion.m,v 1.1 2011/06/22 09:54:19 miquel Exp $', sets, pls);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
509 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
510
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
511 %--------------------------------------------------------------------------
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
512 % Get Default Plist
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
513 %--------------------------------------------------------------------------
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
514 function plout = getDefaultPlist()
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
515 persistent pl;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
516 if exist('pl', 'var')==0 || isempty(pl)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
517 pl = buildplist();
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
518 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
519 plout = pl;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
520 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
521
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
522 function pl = buildplist()
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
523 pl = plist.LPSD_PLIST;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
524 pset(pl,'Navs',1)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
525
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
526 p = plist({'f1', 'Initial frequency for the analysis'}, paramValue.EMPTY_DOUBLE);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
527 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
528
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
529 p = plist({'f2', 'Final frequency for the analysis'}, paramValue.EMPTY_DOUBLE);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
530 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
531
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
532 p = plist({'FitParamas', 'Parameters of the model'}, paramValue.EMPTY_STRING);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
533 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
534
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
535 p = plist({'model','An array of matrix models'}, paramValue.EMPTY_STRING);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
536 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
537
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
538 p = plist({'noise','An array of matrices with the cross-spectrum matrices'}, paramValue.EMPTY_STRING);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
539 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
540
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
541 p = plist({'built-in','Symbolic models of the system as a string of built-in models'}, paramValue.EMPTY_STRING);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
542 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
543
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
544 p = plist({'frequencies','Array of start/sop frequencies where the analysis is performed'}, paramValue.EMPTY_STRING);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
545 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
546
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
547 p = plist({'pinv','Use the Penrose-Moore pseudoinverse'}, paramValue.TRUE_FALSE);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
548 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
549
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
550 p = plist({'tol','Tolerance for the Penrose-Moore pseudoinverse'}, paramValue.EMPTY_DOUBLE);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
551 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
552
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
553 p = plist({'diffStep','Numerical differentiation step for ssm models'}, paramValue.EMPTY_DOUBLE);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
554 pl.append(p);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
555
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
556 end