karjalainen_decay2_fit.m 2.14 KB
Newer Older
1
function [v,norm]=decay2_fit(data,expf,w,figno)
% fit decay_model to dB-scaled energy-time decay curve% data(:,1) = decay data (dB scaled)% data(:,2) = optional time sample points% expf = optional weighting factor% w = weighting vector, default equal weights% figno = optional figure number, 0=no figure% v = model coefficients% norm = error norm from curvefit
sizd = size(data);                      % data dimensionalitysn   = sizd(1);sn2  = ceil(sn/2);                      % size, size/2sn10 = ceil(sn/10);                     % size/10
if (nargin<4)                         % figure to draw to or none=0
    figno = 0; endif nargin<3 || isempty(w)             % weight vector
	w=ones(sn,1);                    endif (nargin<2)                         % expf for model function
    expf=0.4;endif isempty(expf)                      % expf for model function
    expf=0.4; end %y=data(:,1);                          % dB scaled magnitudesif (sizd(2)<2)                        % time value vector default
    x = (1:sizd(1))';else                                  % optionally given time points
    x = data(:,2); endydata   = (10.^(y./20)).^expf;               % magnitude data scalednz      = find(w); bi=nz(1);                 % start index (bi) of nonzero weighty1      = mean(ydata(bi:sn10+bi));           % average of first 10% (after bi)yn      = mean(ydata(sn-sn10+1:sn));         % average of last 10%tmat    = [ones(sn2,1) x(bi:sn2+bi-1)];p       = tmat\y(bi:sn2+bi-1);               % regression line 50% after bitau0    = p(2)/8.7;                          % decay parameter initial valuev0      = [y1 tau0 yn];                      % vector of variables initial valueydata   = w.*ydata;                          % ydata weighted%[v,norm] = lsqcurvefit('karjalainen_decay_model',v0,x,ydata,[],[],[],expf,w); % optimize%                                     % v = model coefficientsif (figno>0)   % plot target data if figno~=0
    plot(x,y); 
    hold on; endz = 20*log10(karjalainen_decay_model(v,x,1));
if (figno>0) % plot fitted data if figno~=0
    plot(x,z,'r');
    hold off;end 
% Call this way: decay2_fit(vp(:,40),[],[],1)