source: trunk/src/peaklock.m @ 446

Last change on this file since 446 was 19, checked in by gostiaux, 15 years ago

the private files have been moved down to the root folder

File size: 7.3 KB
Line 
1%'peaklock': determines peacklocking errors from velocity histograms.
2%-------------------------------------------------------
3%first smooth the input histogram 'histu' in such a way that the integral over
4%n-n+1 is preserved, then deduce the peaklocking 'error' function of the pixcel displacement 'x'.
5%
6% [histinter,x,error]=peaklock(nbb,minim,maxim,histu)
7%OUTPUT:
8%histinter: smoothed interpolated histogram
9% x: vector of displacement values.
10% error: vector of estimated errors corresponding to x
11%INPUT:
12%histu=vector representing the values of histogram  of measured velocity ;
13%minim, maxim: extremal values of the measured velocity (absica for histu)
14%nbb: number of bins inside each integer interval for the histograms
15%SUBROUTINES INCLUDED:
16%spline4.m% spline interpolation at 4th order
17%splinhist.m: give spline coeff cc for a smooth histo (call spline4)
18%histsmooth.m(x,cc): calculate the smooth histo for any value x
19%histder.m(x,cc): calculate the derivative of the smooth histo
20function [histinter,x,error]=peaklock(nbb,minim,maxim,histu)
21
22nint=maxim-minim+1
23xfin=[minim-0.5+1/(2*nbb):(1/nbb):maxim+0.5-(1/(2*nbb))];
24histo=(reshape(histu,nbb,nint));%extract values with x between integer -1/2 integer +1/2
25Integ=sum(histo)/nbb; %integral of the pdf on each integer bin
26[histinter,cc]=splinhist(Integ,minim,nbb);
27histx=reshape(histinter,nbb,nint);
28xint=[minim:1:maxim];
29x=zeros(nbb,nint);
30%determination of the displacement x(j,:)
31%j=1
32delx=histo(1,:)./histsmooth(-0.5*ones(1,nint),cc)/nbb;
33%del(1,:)=delx;
34x(1,:)=-0.5+delx-(delx.*delx/2).*histder(-0.5*ones(1,nint),cc);
35%histx(1,:)=histsmooth(x(j-1,:),cc);
36for j=2:nbb
37    delx=histo(j,:)./histsmooth(x(j-1,:),cc)/nbb;
38    %delx=delx.*(delx<3*ones(1,nint)/nbb)+3*ones(1,nint)/nbb.*~(delx <3*ones(1,nint)/nbb)
39    x(j,:)=x(j-1,:)+delx-(delx.*delx/2).*histder(x(j-1,:),cc);
40end
41%reshape
42xint=ones(nbb,1)*xint;
43x=x+xint;
44x=reshape(x,1,nbb*nint);
45error=xfin+1/(2*nbb)-x;
46
47%-------------------------------------------------------
48% --- determine the spline coefficients cc for the interpolated histogram.
49%-------------------------------------------------
50function [histsmooth,cc]= splinhist(Integ,mini,nbb)
51% provides a smooth histogramm histmooth, which remains always positive,
52% and is such that its sum over each integer bin [i-1/2 i+1/2] is equal to
53% Integ(i). The function determines histmooth as the exponential of a 4th
54% order spline function and adjust the cefficients by a Newton method to
55% fit the integral conditions Integ
56% histmooth is determined at the abscissa
57% xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))] (maxi=mini+size(aa)-1)
58%cc(1-5,i) provides the spline coefficients
59
60% order 0
61siz=size(Integ);
62nint=siz(2);
63izero=find(Integ==0); %indices of zero elements
64inonzero=find(Integ);
65Integ(izero)=min(Integ(inonzero));
66aa=log(Integ);%initial guess for a coeff
67spli=spline4(aa,mini,nbb);  %appel à la fonction spline4
68histsmooth=exp(spli);
69
70S=(sum(reshape(histsmooth,nbb,nint)))/nbb;% integral of the fit histsmooth on ]i-1/2 i+1/2[
71epsilon=max(abs(Integ-S));
72iter=0;
73while epsilon > 0.000001 & iter<10
74ident=eye(nint);
75dSda=ones(nint);
76for j=1:nint% determination of the jacobian matrix dSda
77dhistda=spline4(ident(j,:),mini,nbb);
78expdhistda=dhistda.*histsmooth;
79dSda(j,:)=(sum(reshape(expdhistda,nbb,nint)))/nbb;
80end
81aa=aa+(Integ-S)*inv(dSda);%new estimate of coefficients aa by linear interpolation
82[spli,bb]=spline4(aa,mini,nbb);% new fit histsmooth
83histsmooth=exp(spli);
84S=(sum(reshape(histsmooth,nbb,nint)))/nbb;% integral of the fit histsmooth on ]i-1/2 i+1/2[
85epsilon=max(abs(Integ-S));
86iter=iter+1;
87end
88if iter==10, errordlg('splinhist did not converge after 10 iterations'),end
89cc(1,:)=aa;
90cc(2,:)=bb(1,:);
91cc(3,:)=bb(2,:);
92cc(4,:)=bb(3,:);
93cc(5,:)=bb(4,:);
94
95%-------------------------------------------------------
96% --- determine the 4th order spline coefficients from the function values aa.
97%-------------------------------------------------
98function [histsmooth,bb]= spline4(aa,mini,n)
99% spline interpolation at 4th order
100%aa=vector of values of a function at integer abscissa, starting at mini
101%n=number of subdivisions for the interpolated function
102% histmooth =interpolated values at absissa
103% xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))] (maxi=mini+size(aa)-1)
104%bb=[b(i);c(i);d(i); e(i)] matrix of spline coeff
105L1=[1/2 1/4 1/8 1/16;1 1 3/4 1/2;0 2 3 3;0 0 6 12];
106L2=[-1/2 1/4 -1/8 1/16;1 -1 3/4 -1/2;0 2 -3 3;0 0 6 -12];
107M=inv(L2)*L1;
108[V,D]=eig(M);
109F=-inv(V)*inv(L2)*[1 ;0 ;0;0];
110a1rev=[1 -1/D(1,1)];
111b1rev=[F(1)/D(1,1)];
112a2rev=[1 -1/D(2,2)];
113b2rev=[F(2)/D(2,2)];
114a3=[1 -D(3,3)];
115b3=[F(3)];
116a4=[1 -D(4,4)];
117b4=[F(4)];
118
119%data
120% n=10;% résolution de la pdf: nbre de points par unite de u
121% mini=-10.0;%general mini=uint16(min(values)-1 CHOOSE maxi-mini+1 EVEN
122% maxi=9.0; % general maxi=uint16(max(values))+1
123%nint=double(maxi-mini+1); % nombre d'intervals entiers EVEN!
124siz=size(aa);
125nint=siz(2);
126maxi=mini+nint-1;
127npdf=nint*n;% nbre total d'intervals à introduire dans la pdf: hist(u,npdf)
128%simulation de pdf
129xfin=[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
130%histolin=exp(-(xfin-1).*(xfin-1)).*(2+cos(10*(xfin-1)));% simulation d'une pdf
131%histo=log(histolin);
132%histo=sin(2*pi*xfin);
133%histextract=(reshape(histo,n,nint));
134%aa=sum(histextract)/n %integral of the pdf on each integer bin
135IP=[0 diff(aa)];
136Irev=zeros(size(aa));
137for i=1:nint
138    Irev(i)=aa(end-i+1);
139end
140IPrev=[0 diff(Irev)];
141
142%get the spline coelfficients a_d, using filter on the eigen vectors A,B,C
143Arev=filter(b1rev,a1rev,IPrev);
144Brev=filter(b2rev,a2rev,IPrev);
145C=filter(b3,a3,IP);
146D=filter(b4,a4,IP);
147A=zeros(size(Arev));
148B=zeros(size(Brev));
149for i=1:nint
150    A(i)=Arev(end-i+1);
151    B(i)=Brev(end-i+1);
152end
153%Matr=V*[A;B;C;D];
154bb=V*[A;B;C;D];
155%b=Matr(1,:);
156%c=Matr(2,:);
157%d=Matr(3,:);
158%e=Matr(4,:);
159%a=aa;
160
161%calculate the interpolation using the spline coefficients a-d
162%xextract=(reshape(xfin,n,nint));%
163chi=xfin+1/(2*n)-min(xfin)-double(int16(xfin+(1/(2*n))-min(xfin)))-0.5;% decimal part
164chi2=chi.*chi;
165chi3=chi2.*chi;
166chi4=chi3.*chi;
167avec=reshape(ones(n,1)*aa,1,n*nint);
168bvec=reshape(ones(n,1)*bb(1,:),1,n*nint);
169cvec=reshape(ones(n,1)*bb(2,:),1,n*nint);
170dvec=reshape(ones(n,1)*bb(3,:),1,n*nint);
171evec=reshape(ones(n,1)*bb(4,:),1,n*nint);
172histsmooth=avec+bvec.*chi+cvec.*chi2+dvec.*chi3+evec.*chi4;
173
174%-------------------------------------------------------
175% --- determine the interpolated histogram at points chi from the spline ceff cc.
176%-------------------------------------------------
177function histx= histsmooth(chi,cc)
178% provides the value of the interpolated histogram at values chi=x-i
179%(difference with the mnearest integer)
180% cc(5,size(chi)) is the set of spline coefficients obtained by splinhist
181chi2=chi.*chi;
182chi3=chi2.*chi;
183chi4=chi3.*chi;
184histx=exp(cc(1,:)+cc(2,:).*chi+cc(3,:).*chi2+cc(4,:).*chi3+cc(5,:).*chi4);
185
186%-------------------------------------------------------
187% --- determine the derivative p'/p of the interpolated histogram at points chi from the spline ceff cc.
188%-------------------------------------------------
189function histder= histder(chi,cc)
190% provides the logarithmique derivative p'/p of the interpolated histogram
191%at values chi=x-i
192%(difference with the nearest integer)
193% cc(5,size(chi)) is the set of spline coefficients obtained by splinhist
194chi2=chi.*chi;
195chi3=chi2.*chi;
196chi4=chi3.*chi;
197histder=cc(2,:)+2*cc(3,:).*chi+3*cc(4,:).*chi2+4*cc(5,:).*chi3;
Note: See TracBrowser for help on using the repository browser.