0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 % Compute Fisher matrix
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 % Parameters are:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 % i1 - input 1st channel (ao)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 % n - noise both channels (matrix 1x1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 % mdl - model (matrix or ssm)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 % params - parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 % numparams - numerical value of parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 % freqs - frequnecies being evaluated
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % N - number of fft frequencies
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 % pl - plist
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 % M Nofrarias 20-09-11
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 % $Id: fisher_1x1.m,v 1.1 2011/10/07 08:17:52 miquel Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 function FisMat = fisher_1x1(i1,n,mdl,params,numparams,freqs,N,pl,inNames,outNames)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 import utils.const.*
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 % Compute psd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 n1 = psd(n.getObjectAtIndex(1,1), pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 % interpolate to given frequencies
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 % noise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 S11 = interp(n1,plist('vertices',freqs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 % get some parameters used below
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 fs = S11.fs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 if ~isempty(mdl) && all(strcmp(class(mdl),'matrix'))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 % compute built-in matrix
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 for i = 1:numel(mdl.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 % set Xvals
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 h(i) = mdl.getObjectAtIndex(i).setXvals(freqs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 % set alias
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 h(i).assignalias(mdl.objs(i),plist('xvals',freqs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 % set paramaters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 h(i).setParams(params,numparams);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 % differentiate and eval
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 for i = 1:length(params)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 utils.helper.msg(msg.IMPORTANT, sprintf('computing symbolic differentiation with respect %s',params{i}), mfilename('class'), mfilename);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 % differentiate symbolically
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 dH11 = diff(h(1),params{i});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 % evaluate
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 d11(i) = eval(dH11);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 elseif ~isempty(mdl) && all(strcmp(class(mdl),'ssm'))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 meval = copy(mdl,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 % set parameter values
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 % meval.doSetParameters(params, numparams);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 meval.setParameters(params, numparams);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 % get the differentiation step
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 step = find(pl,'step');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 % case no diff. step introduced
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 if isempty(step)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 utils.helper.msg(msg.IMPORTANT, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 sprintf('computing optimal differentiation steps'), mfilename('class'), mfilename);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 ranges = find(pl,'stepRanges');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 if isempty(ranges)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 error('### Please input upper and lower ranges for the parameters: ''ranges''')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 ngrid = find(pl,'ngrid');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 if isempty(ngrid)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 error('### Please input a number of points for the grid to compute the diff. step : ''ngrid''')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 % look for numerical differentiation step
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 step = utils.math.diffStepFish_1x1(i1,S11,N,meval,params,numparams,ngrid,ranges,freqs,inNames,outNames);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 % differentiate and eval
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 for i = 1:length(params)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 utils.helper.msg(msg.IMPORTANT, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 sprintf('computing numerical differentiation with respect %s, Step:%4.2d ',params{i},step(i)), mfilename('class'), mfilename);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 % differentiate numerically
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 dH = meval.parameterDiff(plist('names', params(i),'values',step(i)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 % create plist with correct outNames (since parameterDiff change them)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 out1 = strrep(outNames{1},'.', sprintf('_DIFF_%s.',params{i})); % 2x2 case
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 spl = plist('set', 'for bode', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 'outputs', out1, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 'inputs', inNames, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 'reorganize', true,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 'f', freqs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 % do bode
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 d = bode(dH, spl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 % 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)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 d11(i) = d.objs(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 error('### please introduce models for the transfer functions')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 % scaling of PSD
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 % PSD = 2/(N*fs) * FFT *conj(FFT)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 C11 = N*fs/2.*S11.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 % compute elements of inverse cross-spectrum matrix
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105 InvS11 = 1./C11;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 % compute Fisher Matrix
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 for i =1:length(params)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 for j =1:length(params)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 v1v1 = conj(d11(i).y.*i1.y).*(d11(j).y.*i1.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113 FisMat(i,j) = sum(real(InvS11.*v1v1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117 end |