Description of peaklock

0001 %'peaklock': determines peacklocking errors from velocity histograms.
0002 %-------------------------------------------------------
0003 %first smooth the input histogram 'histu' in such a way that the integral over
0004 %n-n+1 is preserved, then deduce the peaklocking 'error' function of the pixcel displacement 'x'.
0005 %
0006 % [histinter,x,error]=peaklock(nbb,minim,maxim,histu)
0007 %OUTPUT:
0008 %histinter: smoothed interpolated histogram
0009 % x: vector of displacement values.
0010 % error: vector of estimated errors corresponding to x
0011 %INPUT:
0012 %histu=vector representing the values of histogram  of measured velocity ;
0013 %minim, maxim: extremal values of the measured velocity (absica for histu)
0014 %nbb: number of bins inside each integer interval for the histograms
0015 %SUBROUTINES INCLUDED:
0016 %spline4.m% spline interpolation at 4th order
0017 %splinhist.m: give spline coeff cc for a smooth histo (call spline4)
0018 %histsmooth.m(x,cc): calculate the smooth histo for any value x
0019 %histder.m(x,cc): calculate the derivative of the smooth histo
0020 function [histinter,x,error]=peaklock(nbb,minim,maxim,histu)
0021 
0022 nint=maxim-minim+1
0023 xfin=[minim-0.5+1/(2*nbb):(1/nbb):maxim+0.5-(1/(2*nbb))];
0024 histo=(reshape(histu,nbb,nint));%extract values with x between integer -1/2 integer +1/2
0025 Integ=sum(histo)/nbb; %integral of the pdf on each integer bin
0026 [histinter,cc]=splinhist(Integ,minim,nbb);
0027 histx=reshape(histinter,nbb,nint);
0028 xint=[minim:1:maxim];
0029 x=zeros(nbb,nint);
0030 %determination of the displacement x(j,:)
0031 %j=1
0032 delx=histo(1,:)./histsmooth(-0.5*ones(1,nint),cc)/nbb;
0033 %del(1,:)=delx;
0034 x(1,:)=-0.5+delx-(delx.*delx/2).*histder(-0.5*ones(1,nint),cc);
0035 %histx(1,:)=histsmooth(x(j-1,:),cc);
0036 for j=2:nbb
0037     delx=histo(j,:)./histsmooth(x(j-1,:),cc)/nbb;
0038     %delx=delx.*(delx<3*ones(1,nint)/nbb)+3*ones(1,nint)/nbb.*~(delx <3*ones(1,nint)/nbb)
0039     x(j,:)=x(j-1,:)+delx-(delx.*delx/2).*histder(x(j-1,:),cc);
0040 end
0041 %reshape
0042 xint=ones(nbb,1)*xint;
0043 x=x+xint;
0044 x=reshape(x,1,nbb*nint);
0045 error=xfin+1/(2*nbb)-x;
0046 
0047 %-------------------------------------------------------
0048 % --- determine the spline coefficients cc for the interpolated histogram.
0049 %-------------------------------------------------
0050 function [histsmooth,cc]= splinhist(Integ,mini,nbb)
0051 % provides a smooth histogramm histmooth, which remains always positive,
0052 % and is such that its sum over each integer bin [i-1/2 i+1/2] is equal to
0053 % Integ(i). The function determines histmooth as the exponential of a 4th
0054 % order spline function and adjust the cefficients by a Newton method to
0055 % fit the integral conditions Integ
0056 % histmooth is determined at the abscissa
0057 % xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))] (maxi=mini+size(aa)-1)
0058 %cc(1-5,i) provides the spline coefficients
0059 
0060 % order 0
0061 siz=size(Integ);
0062 nint=siz(2);
0063 izero=find(Integ==0); %indices of zero elements
0064 inonzero=find(Integ);
0065 Integ(izero)=min(Integ(inonzero));
0066 aa=log(Integ);%initial guess for a coeff
0067 spli=spline4(aa,mini,nbb);  %appel à la fonction spline4
0068 histsmooth=exp(spli);
0069 
0070 S=(sum(reshape(histsmooth,nbb,nint)))/nbb;% integral of the fit histsmooth on ]i-1/2 i+1/2[
0071 epsilon=max(abs(Integ-S));
0072 iter=0;
0073 while epsilon > 0.000001 & iter<10
0074 ident=eye(nint);
0075 dSda=ones(nint);
0076 for j=1:nint% determination of the jacobian matrix dSda
0077 dhistda=spline4(ident(j,:),mini,nbb);
0078 expdhistda=dhistda.*histsmooth;
0079 dSda(j,:)=(sum(reshape(expdhistda,nbb,nint)))/nbb;
0080 end
0081 aa=aa+(Integ-S)*inv(dSda);%new estimate of coefficients aa by linear interpolation
0082 [spli,bb]=spline4(aa,mini,nbb);% new fit histsmooth
0083 histsmooth=exp(spli);
0084 S=(sum(reshape(histsmooth,nbb,nint)))/nbb;% integral of the fit histsmooth on ]i-1/2 i+1/2[
0085 epsilon=max(abs(Integ-S));
0086 iter=iter+1;
0087 end
0088 if iter==10, errordlg('splinhist did not converge after 10 iterations'),end
0089 cc(1,:)=aa;
0090 cc(2,:)=bb(1,:);
0091 cc(3,:)=bb(2,:);
0092 cc(4,:)=bb(3,:);
0093 cc(5,:)=bb(4,:);
0094 
0095 %-------------------------------------------------------
0096 % --- determine the 4th order spline coefficients from the function values aa.
0097 %-------------------------------------------------
0098 function [histsmooth,bb]= spline4(aa,mini,n)
0099 % spline interpolation at 4th order
0100 %aa=vector of values of a function at integer abscissa, starting at mini
0101 %n=number of subdivisions for the interpolated function
0102 % histmooth =interpolated values at absissa
0103 % xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))] (maxi=mini+size(aa)-1)
0104 %bb=[b(i);c(i);d(i); e(i)] matrix of spline coeff
0105 L1=[1/2 1/4 1/8 1/16;1 1 3/4 1/2;0 2 3 3;0 0 6 12];
0106 L2=[-1/2 1/4 -1/8 1/16;1 -1 3/4 -1/2;0 2 -3 3;0 0 6 -12];
0107 M=inv(L2)*L1;
0108 [V,D]=eig(M);
0109 F=-inv(V)*inv(L2)*[1 ;0 ;0;0];
0110 a1rev=[1 -1/D(1,1)];
0111 b1rev=[F(1)/D(1,1)];
0112 a2rev=[1 -1/D(2,2)];
0113 b2rev=[F(2)/D(2,2)];
0114 a3=[1 -D(3,3)];
0115 b3=[F(3)];
0116 a4=[1 -D(4,4)];
0117 b4=[F(4)];
0118 
0119 %data
0120 % n=10;% résolution de la pdf: nbre de points par unite de u
0121 % mini=-10.0;%general mini=uint16(min(values)-1 CHOOSE maxi-mini+1 EVEN
0122 % maxi=9.0; % general maxi=uint16(max(values))+1
0123 %nint=double(maxi-mini+1); % nombre d'intervals entiers EVEN!
0124 siz=size(aa);
0125 nint=siz(2);
0126 maxi=mini+nint-1;
0127 npdf=nint*n;% nbre total d'intervals à introduire dans la pdf: hist(u,npdf)
0128 %simulation de pdf
0129 xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))];% valeurs d'interpolation: we take n values in each integer interval
0130 %histolin=exp(-(xfin-1).*(xfin-1)).*(2+cos(10*(xfin-1)));% simulation d'une pdf
0131 %histo=log(histolin);
0132 %histo=sin(2*pi*xfin);
0133 %histextract=(reshape(histo,n,nint));
0134 %aa=sum(histextract)/n %integral of the pdf on each integer bin
0135 IP=[0 diff(aa)];
0136 Irev=zeros(size(aa));
0137 for i=1:nint
0138     Irev(i)=aa(end-i+1);
0139 end
0140 IPrev=[0 diff(Irev)];
0141 
0142 %get the spline coelfficients a_d, using filter on the eigen vectors A,B,C
0143 Arev=filter(b1rev,a1rev,IPrev);
0144 Brev=filter(b2rev,a2rev,IPrev);
0145 C=filter(b3,a3,IP);
0146 D=filter(b4,a4,IP);
0147 A=zeros(size(Arev));
0148 B=zeros(size(Brev));
0149 for i=1:nint
0150     A(i)=Arev(end-i+1);
0151     B(i)=Brev(end-i+1);
0152 end
0153 %Matr=V*[A;B;C;D];
0154 bb=V*[A;B;C;D];
0155 %b=Matr(1,:);
0156 %c=Matr(2,:);
0157 %d=Matr(3,:);
0158 %e=Matr(4,:);
0159 %a=aa;
0160 
0161 %calculate the interpolation using the spline coefficients a-d
0162 %xextract=(reshape(xfin,n,nint));%
0163 chi=xfin+1/(2*n)-min(xfin)-double(int16(xfin+(1/(2*n))-min(xfin)))-0.5;% decimal part
0164 chi2=chi.*chi;
0165 chi3=chi2.*chi;
0166 chi4=chi3.*chi;
0167 avec=reshape(ones(n,1)*aa,1,n*nint);
0168 bvec=reshape(ones(n,1)*bb(1,:),1,n*nint);
0169 cvec=reshape(ones(n,1)*bb(2,:),1,n*nint);
0170 dvec=reshape(ones(n,1)*bb(3,:),1,n*nint);
0171 evec=reshape(ones(n,1)*bb(4,:),1,n*nint);
0172 histsmooth=avec+bvec.*chi+cvec.*chi2+dvec.*chi3+evec.*chi4;
0173 
0174 %-------------------------------------------------------
0175 % --- determine the interpolated histogram at points chi from the spline ceff cc.
0176 %-------------------------------------------------
0177 function histx= histsmooth(chi,cc)
0178 % provides the value of the interpolated histogram at values chi=x-i
0179 %(difference with the mnearest integer)
0180 % cc(5,size(chi)) is the set of spline coefficients obtained by splinhist
0181 chi2=chi.*chi;
0182 chi3=chi2.*chi;
0183 chi4=chi3.*chi;
0184 histx=exp(cc(1,:)+cc(2,:).*chi+cc(3,:).*chi2+cc(4,:).*chi3+cc(5,:).*chi4);
0185 
0186 %-------------------------------------------------------
0187 % --- determine the derivative p'/p of the interpolated histogram at points chi from the spline ceff cc.
0188 %-------------------------------------------------
0189 function histder= histder(chi,cc)
0190 % provides the logarithmique derivative p'/p of the interpolated histogram
0191 %at values chi=x-i
0192 %(difference with the nearest integer)
0193 % cc(5,size(chi)) is the set of spline coefficients obtained by splinhist
0194 chi2=chi.*chi;
0195 chi3=chi2.*chi;
0196 chi4=chi3.*chi;
0197 histder=cc(2,:)+2*cc(3,:).*chi+3*cc(4,:).*chi2+4*cc(5,:).*chi3;
peaklock

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
