Mercurial > hg > ltpda
view m-toolbox/classes/+utils/@math/diffStepFish.m @ 46:ca0b8d4dcdb6 database-connection-manager
Fix
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Tue, 06 Dec 2011 19:07:27 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Look for differentiation step for a given parameter and % % Parameters are: % % $Id: diffStepFish.m,v 1.2 2011/09/19 06:17:45 miquel Exp $ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function best = diffStepFish(i1,i2,S11,S12,S21,S22,N,meval,params,ngrid,ranges,freqs,inNames,outNames) % remove aux file if existing if exist('diffStepFish.txt') == 2 ! rm diffStepFish.txt end step = ones(ngrid,numel(params)); % build matrix of steps % for ii = 1:length(params) % step(:,ii) = [] logspace(ranges(1,ii),ranges(2,ii),ngrid); % end for ii = 1:ngrid step(ii,:) = ranges(1,:); end % step(:,1) = logspace(ranges(1,1),ranges(2,1),ngrid); for kk = 1:length(params) step(:,kk) = logspace(log10(ranges(1,kk)),log10(ranges(2,kk)),ngrid); Rmat = []; for jj = 1:ngrid for ii = 1:length(params) tic % differentiate numerically dH = meval.parameterDiff(plist('names', params(ii),'values',step(jj,ii))); % create plist with correct outNames (since parameterDiff change them) out1 = strrep(outNames{1},'.', sprintf('_DIFF_%s.',params{ii})); % 2x2 case out2 =strrep(outNames{2},'.', sprintf('_DIFF_%s.',params{ii})); spl = plist('set', 'for bode', ... 'outputs', {out1,out2}, ... 'inputs', inNames, ... 'reorganize', true,... 'f', freqs); % do bode d = bode(dH, spl); % 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) d11(ii) = d.objs(1); d21(ii) = d.objs(2); d12(ii) = d.objs(3); d22(ii) = d.objs(4); end fs = S11.fs; % scaling of PSD % PSD = 2/(N*fs) * FFT *conj(FFT) C11 = N*fs/2.*S11.y; C22 = N*fs/2.*S22.y; C12 = N*fs/2.*S12.y; C21 = N*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 detFisMat = det(FisMat); R = [step(jj,:) detFisMat]; save('diffStepFish.txt','R','-ascii','-append'); Rmat = [Rmat; R]; toc end % look for the stable step: compute diff and % look for the smallest one in absolute value % The smallest slope marks the plateau diffDetFisMat = abs(diff(Rmat(:,end))); lowdet = diffDetFisMat(1); ind = 2; for k = 1:numel(diffDetFisMat) if diffDetFisMat(k) < lowdet lowdet = diffDetFisMat(k); ind = k+1; % index give by diff = x(2) - x(1). We take the step corresponding to x(2) end end step(:,kk) = step(jj,kk)*ones(ngrid,1); end step(:,end) = logspace(log10(ranges(1,end)),log10(ranges(2,end)),ngrid); best = step(1,:); end