Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/fisher_1x1.m @ 0:f0afece42f48

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