Home > . > peaklock.m

# peaklock

## PURPOSE

'peaklock': determines peacklocking errors from velocity histograms.

## SYNOPSIS

function [histinter,x,error]=peaklock(nbb,minim,maxim,histu)

## DESCRIPTION

    'peaklock': determines peacklocking errors from velocity histograms.
    -------------------------------------------------------
    first smooth the input histogram 'histu' in such a way that the integral over
    n-n+1 is preserved, then deduce the peaklocking 'error' function of the pixcel displacement 'x'.

     [histinter,x,error]=peaklock(nbb,minim,maxim,histu)
    OUTPUT:
    histinter: smoothed interpolated histogram
     x: vector of displacement values.
     error: vector of estimated errors corresponding to x
    INPUT:
    histu=vector representing the values of histogram  of measured velocity ;
    minim, maxim: extremal values of the measured velocity (absica for histu)
    nbb: number of bins inside each integer interval for the histograms
    SUBROUTINES INCLUDED:
    spline4.m% spline interpolation at 4th order
    splinhist.m: give spline coeff cc for a smooth histo (call spline4)
    histsmooth.m(x,cc): calculate the smooth histo for any value x
    histder.m(x,cc): calculate the derivative of the smooth histo

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- get_field 'get_field': display variables and attributes from a Netcdf file, and plot selected fields

## SUBFUNCTIONS

- function [histsmooth,cc]= splinhist(Integ,mini,nbb)
- function [histsmooth,bb]= spline4(aa,mini,n)
- function histx= histsmooth(chi,cc)
- function histder= histder(chi,cc)

## SOURCE CODE

    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;