Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/diffStepFish_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 % Look for differentiation step for a given parameter and
4 %
5 % Parameters are:
6
7 %
8 % $Id: diffStepFish_1x1.m,v 1.1 2011/10/07 08:17:52 miquel Exp $
9 %
10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
11 function best = diffStepFish_1x1(i1,S11,N,meval,params,numparams,ngrid,ranges,freqs,inNames,outNames)
12
13 import utils.const.*
14
15 % remove aux file if existing
16 if exist('diffStepFish.txt') == 2
17 ! rm diffStepFish.txt
18 end
19
20 step = ones(ngrid,numel(params));
21
22 % initialize matrix of steps
23 for ii = 1:numel(params)
24 step(:,ii) = ranges(1)*numparams(ii);
25 end
26
27 % step(:,1) = logspace(ranges(1,1),ranges(2,1),ngrid);
28
29 for kk = 1:length(params)
30 step(:,kk) = numparams(kk)*logspace(log10(ranges(1)),log10(ranges(2)),ngrid);
31 Rmat = [];
32 for jj = 1:ngrid
33 for ii = 1:length(params)
34 % differentiate numerically
35 dH = meval.parameterDiff(plist('names', params(ii),'values',step(jj,ii)));
36 % create plist with correct outNames (since parameterDiff change them)
37 out1 = strrep(outNames{1},'.', sprintf('_DIFF_%s.',params{ii})); % 2x2 case
38 spl = plist('set', 'for bode', ...
39 'outputs', {out1}, ...
40 'inputs', inNames, ...
41 'reorganize', true,...
42 'f', freqs);
43 % do bode
44 d = bode(dH, spl);
45 % assign according matlab's matrix notation:
46 % H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
47 d11(ii) = d.objs(1);
48 end
49
50 fs = S11.fs;
51 % scaling of PSD
52 % PSD = 2/(N*fs) * FFT *conj(FFT)
53 C11 = N*fs/2.*S11.y;
54
55 % compute elements of inverse cross-spectrum matrix
56 InvS11 = 1./C11;
57
58 % compute Fisher Matrix
59 for i =1:length(params)
60 for j =1:length(params)
61
62 v1v1 = conj(d11(i).y.*i1.y).*(d11(j).y.*i1.y);
63 FisMat(i,j) = sum(real(InvS11.*v1v1));
64 end
65 end
66
67 detFisMat = det(FisMat);
68 R = [step(jj,:) detFisMat];
69 % only file diffStepFish.txt stores all iterations. Rmat is
70 % initialized for each loop
71 save('diffStepFish.txt','R','-ascii','-append');
72 Rmat = [Rmat; R];
73 end
74
75 % look for the stable step: compute diff and
76 % look for the smallest one in absolute value
77 % The smallest slope marks the plateau
78 diffDetFisMat = abs(diff(Rmat(:,end)));
79 lowdet = diffDetFisMat(1);
80 ind = 2;
81 for k = 1:numel(diffDetFisMat)
82 if diffDetFisMat(k) < lowdet
83 lowdet = diffDetFisMat(k);
84 ind = k+1; % index give by diff = x(2) - x(1). We take the step corresponding to x(2)
85 end
86 end
87 % display message
88 utils.helper.msg(msg.IMPORTANT, ...
89 sprintf('Best numerical diff. step with respect %s: %d',params{kk}, step(ind,kk)), mfilename('class'), mfilename);
90 % reassing all current column to the best step
91 step(:,kk) = step(ind,kk)*ones(ngrid,1);
92
93 figure
94 diffDetFisMat(diffDetFisMat == 0) = 1e-20; % to avoid zeros in loglog plot
95 loglog(Rmat(1:end-1,kk)/numparams(kk),diffDetFisMat,'--ks','LineWidth',2,'MarkerSize',10)
96 title(sprintf('Parameter: %s',params{kk}))
97 ylabel('\Delta FisMat / \Delta\theta')
98 xlabel('Normalised \Delta\theta')
99
100 end
101 best = step(1,:);
102 end
103