+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−%
+
−%  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
+
−