source: trunk/src/uvmat_doc/FUNCTIONS_DOC/peaklock.html @ 48

Last change on this file since 48 was 37, checked in by sommeria, 15 years ago

create_grid.fig ,
uvmat_doc and all the included files added
File size: 16.5 KB
Line 
1<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
2                "http://www.w3.org/TR/REC-html40/loose.dtd">
3<html>
4<head>
5  <title>Description of peaklock</title>
6  <meta name="keywords" content="peaklock">
7  <meta name="description" content="'peaklock': determines peacklocking errors from velocity histograms.">
8  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
9  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
10  <meta name="robots" content="index, follow">
11  <link type="text/css" rel="stylesheet" href="../m2html.css">
12</head>
13<body>
14<a name="_top"></a>
15<div><a href="../index.html">Home</a> &gt;  <a href="index.html">.</a> &gt; peaklock.m</div>
16
17<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
18<td align="right"><a href="index.html">Index for .&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->
19
20<h1>peaklock
21</h1>
22
23<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
24<div class="box"><strong>'peaklock': determines peacklocking errors from velocity histograms.</strong></div>
25
26<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
27<div class="box"><strong>function [histinter,x,error]=peaklock(nbb,minim,maxim,histu) </strong></div>
28
29<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
30<div class="fragment"><pre class="comment">'peaklock': determines peacklocking errors from velocity histograms.
31-------------------------------------------------------
32first smooth the input histogram 'histu' in such a way that the integral over
33n-n+1 is preserved, then deduce the peaklocking 'error' function of the pixcel displacement 'x'.
34
35 [histinter,x,error]=peaklock(nbb,minim,maxim,histu)
36OUTPUT:
37histinter: smoothed interpolated histogram
38 x: vector of displacement values.
39 error: vector of estimated errors corresponding to x
40INPUT:
41histu=vector representing the values of histogram  of measured velocity ;
42minim, maxim: extremal values of the measured velocity (absica for histu)
43nbb: number of bins inside each integer interval for the histograms
44SUBROUTINES INCLUDED:
45spline4.m% spline interpolation at 4th order
46splinhist.m: give spline coeff cc for a smooth histo (call spline4)
47histsmooth.m(x,cc): calculate the smooth histo for any value x
48histder.m(x,cc): calculate the derivative of the smooth histo</pre></div>
49
50<!-- crossreference -->
51<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
52This function calls:
53<ul style="list-style-image:url(../matlabicon.gif)">
54</ul>
55This function is called by:
56<ul style="list-style-image:url(../matlabicon.gif)">
57<li><a href="get_field.html" class="code" title="function varargout = get_field(varargin)">get_field</a>        'get_field': display variables and attributes from a Netcdf file, and plot selected fields</li></ul>
58<!-- crossreference -->
59
60<h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
61<ul style="list-style-image:url(../matlabicon.gif)">
62<li><a href="#_sub1" class="code">function [histsmooth,cc]= splinhist(Integ,mini,nbb)</a></li><li><a href="#_sub2" class="code">function [histsmooth,bb]= spline4(aa,mini,n)</a></li><li><a href="#_sub3" class="code">function histx= histsmooth(chi,cc)</a></li><li><a href="#_sub4" class="code">function histder= histder(chi,cc)</a></li></ul>
63<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
64<div class="fragment"><pre>0001 <span class="comment">%'peaklock': determines peacklocking errors from velocity histograms.</span>
650002 <span class="comment">%-------------------------------------------------------</span>
660003 <span class="comment">%first smooth the input histogram 'histu' in such a way that the integral over</span>
670004 <span class="comment">%n-n+1 is preserved, then deduce the peaklocking 'error' function of the pixcel displacement 'x'.</span>
680005 <span class="comment">%</span>
690006 <span class="comment">% [histinter,x,error]=peaklock(nbb,minim,maxim,histu)</span>
700007 <span class="comment">%OUTPUT:</span>
710008 <span class="comment">%histinter: smoothed interpolated histogram</span>
720009 <span class="comment">% x: vector of displacement values.</span>
730010 <span class="comment">% error: vector of estimated errors corresponding to x</span>
740011 <span class="comment">%INPUT:</span>
750012 <span class="comment">%histu=vector representing the values of histogram  of measured velocity ;</span>
760013 <span class="comment">%minim, maxim: extremal values of the measured velocity (absica for histu)</span>
770014 <span class="comment">%nbb: number of bins inside each integer interval for the histograms</span>
780015 <span class="comment">%SUBROUTINES INCLUDED:</span>
790016 <span class="comment">%spline4.m% spline interpolation at 4th order</span>
800017 <span class="comment">%splinhist.m: give spline coeff cc for a smooth histo (call spline4)</span>
810018 <span class="comment">%histsmooth.m(x,cc): calculate the smooth histo for any value x</span>
820019 <span class="comment">%histder.m(x,cc): calculate the derivative of the smooth histo</span>
830020 <a name="_sub0" href="#_subfunctions" class="code">function [histinter,x,error]=peaklock(nbb,minim,maxim,histu)</a>
840021
850022 nint=maxim-minim+1
860023 xfin=[minim-0.5+1/(2*nbb):(1/nbb):maxim+0.5-(1/(2*nbb))];
870024 histo=(reshape(histu,nbb,nint));<span class="comment">%extract values with x between integer -1/2 integer +1/2</span>
880025 Integ=sum(histo)/nbb; <span class="comment">%integral of the pdf on each integer bin</span>
890026 [histinter,cc]=<a href="#_sub1" class="code" title="subfunction [histsmooth,cc]= splinhist(Integ,mini,nbb)">splinhist</a>(Integ,minim,nbb);
900027 histx=reshape(histinter,nbb,nint);
910028 xint=[minim:1:maxim];
920029 x=zeros(nbb,nint);
930030 <span class="comment">%determination of the displacement x(j,:)</span>
940031 <span class="comment">%j=1</span>
950032 delx=histo(1,:)./<a href="#_sub3" class="code" title="subfunction histx= histsmooth(chi,cc)">histsmooth</a>(-0.5*ones(1,nint),cc)/nbb;
960033 <span class="comment">%del(1,:)=delx;</span>
970034 x(1,:)=-0.5+delx-(delx.*delx/2).*<a href="#_sub4" class="code" title="subfunction histder= histder(chi,cc)">histder</a>(-0.5*ones(1,nint),cc);
980035 <span class="comment">%histx(1,:)=histsmooth(x(j-1,:),cc);</span>
990036 <span class="keyword">for</span> j=2:nbb
1000037     delx=histo(j,:)./<a href="#_sub3" class="code" title="subfunction histx= histsmooth(chi,cc)">histsmooth</a>(x(j-1,:),cc)/nbb;
1010038     <span class="comment">%delx=delx.*(delx&lt;3*ones(1,nint)/nbb)+3*ones(1,nint)/nbb.*~(delx &lt;3*ones(1,nint)/nbb)</span>
1020039     x(j,:)=x(j-1,:)+delx-(delx.*delx/2).*<a href="#_sub4" class="code" title="subfunction histder= histder(chi,cc)">histder</a>(x(j-1,:),cc);
1030040 <span class="keyword">end</span>
1040041 <span class="comment">%reshape</span>
1050042 xint=ones(nbb,1)*xint;
1060043 x=x+xint;
1070044 x=reshape(x,1,nbb*nint);
1080045 error=xfin+1/(2*nbb)-x;
1090046
1100047 <span class="comment">%-------------------------------------------------------</span>
1110048 <span class="comment">% --- determine the spline coefficients cc for the interpolated histogram.</span>
1120049 <span class="comment">%-------------------------------------------------</span>
1130050 <a name="_sub1" href="#_subfunctions" class="code">function [histsmooth,cc]= splinhist(Integ,mini,nbb)</a>
1140051 <span class="comment">% provides a smooth histogramm histmooth, which remains always positive,</span>
1150052 <span class="comment">% and is such that its sum over each integer bin [i-1/2 i+1/2] is equal to</span>
1160053 <span class="comment">% Integ(i). The function determines histmooth as the exponential of a 4th</span>
1170054 <span class="comment">% order spline function and adjust the cefficients by a Newton method to</span>
1180055 <span class="comment">% fit the integral conditions Integ</span>
1190056 <span class="comment">% histmooth is determined at the abscissa</span>
1200057 <span class="comment">% xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))] (maxi=mini+size(aa)-1)</span>
1210058 <span class="comment">%cc(1-5,i) provides the spline coefficients</span>
1220059
1230060 <span class="comment">% order 0</span>
1240061 siz=size(Integ);
1250062 nint=siz(2);
1260063 izero=find(Integ==0); <span class="comment">%indices of zero elements</span>
1270064 inonzero=find(Integ);
1280065 Integ(izero)=min(Integ(inonzero));
1290066 aa=log(Integ);<span class="comment">%initial guess for a coeff</span>
1300067 spli=<a href="#_sub2" class="code" title="subfunction [histsmooth,bb]= spline4(aa,mini,n)">spline4</a>(aa,mini,nbb);  <span class="comment">%appel à la fonction spline4</span>
1310068 <a href="#_sub3" class="code" title="subfunction histx= histsmooth(chi,cc)">histsmooth</a>=exp(spli);
1320069
1330070 S=(sum(reshape(<a href="#_sub3" class="code" title="subfunction histx= histsmooth(chi,cc)">histsmooth</a>,nbb,nint)))/nbb;<span class="comment">% integral of the fit histsmooth on ]i-1/2 i+1/2[</span>
1340071 epsilon=max(abs(Integ-S));
1350072 iter=0;
1360073 <span class="keyword">while</span> epsilon &gt; 0.000001 &amp; iter&lt;10
1370074 ident=eye(nint);
1380075 dSda=ones(nint);
1390076 <span class="keyword">for</span> j=1:nint<span class="comment">% determination of the jacobian matrix dSda</span>
1400077 dhistda=<a href="#_sub2" class="code" title="subfunction [histsmooth,bb]= spline4(aa,mini,n)">spline4</a>(ident(j,:),mini,nbb);
1410078 expdhistda=dhistda.*<a href="#_sub3" class="code" title="subfunction histx= histsmooth(chi,cc)">histsmooth</a>;
1420079 dSda(j,:)=(sum(reshape(expdhistda,nbb,nint)))/nbb;
1430080 <span class="keyword">end</span>
1440081 aa=aa+(Integ-S)*inv(dSda);<span class="comment">%new estimate of coefficients aa by linear interpolation</span>
1450082 [spli,bb]=<a href="#_sub2" class="code" title="subfunction [histsmooth,bb]= spline4(aa,mini,n)">spline4</a>(aa,mini,nbb);<span class="comment">% new fit histsmooth</span>
1460083 <a href="#_sub3" class="code" title="subfunction histx= histsmooth(chi,cc)">histsmooth</a>=exp(spli);
1470084 S=(sum(reshape(<a href="#_sub3" class="code" title="subfunction histx= histsmooth(chi,cc)">histsmooth</a>,nbb,nint)))/nbb;<span class="comment">% integral of the fit histsmooth on ]i-1/2 i+1/2[</span>
1480085 epsilon=max(abs(Integ-S));
1490086 iter=iter+1;
1500087 <span class="keyword">end</span>
1510088 <span class="keyword">if</span> iter==10, errordlg(<span class="string">'splinhist did not converge after 10 iterations'</span>),<span class="keyword">end</span>
1520089 cc(1,:)=aa;
1530090 cc(2,:)=bb(1,:);
1540091 cc(3,:)=bb(2,:);
1550092 cc(4,:)=bb(3,:);
1560093 cc(5,:)=bb(4,:);
1570094
1580095 <span class="comment">%-------------------------------------------------------</span>
1590096 <span class="comment">% --- determine the 4th order spline coefficients from the function values aa.</span>
1600097 <span class="comment">%-------------------------------------------------</span>
1610098 <a name="_sub2" href="#_subfunctions" class="code">function [histsmooth,bb]= spline4(aa,mini,n)</a>
1620099 <span class="comment">% spline interpolation at 4th order</span>
1630100 <span class="comment">%aa=vector of values of a function at integer abscissa, starting at mini</span>
1640101 <span class="comment">%n=number of subdivisions for the interpolated function</span>
1650102 <span class="comment">% histmooth =interpolated values at absissa</span>
1660103 <span class="comment">% xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))] (maxi=mini+size(aa)-1)</span>
1670104 <span class="comment">%bb=[b(i);c(i);d(i); e(i)] matrix of spline coeff</span>
1680105 L1=[1/2 1/4 1/8 1/16;1 1 3/4 1/2;0 2 3 3;0 0 6 12];
1690106 L2=[-1/2 1/4 -1/8 1/16;1 -1 3/4 -1/2;0 2 -3 3;0 0 6 -12];
1700107 M=inv(L2)*L1;
1710108 [V,D]=eig(M);
1720109 F=-inv(V)*inv(L2)*[1 ;0 ;0;0];
1730110 a1rev=[1 -1/D(1,1)];
1740111 b1rev=[F(1)/D(1,1)];
1750112 a2rev=[1 -1/D(2,2)];
1760113 b2rev=[F(2)/D(2,2)];
1770114 a3=[1 -D(3,3)];
1780115 b3=[F(3)];
1790116 a4=[1 -D(4,4)];
1800117 b4=[F(4)];
1810118
1820119 <span class="comment">%data</span>
1830120 <span class="comment">% n=10;% résolution de la pdf: nbre de points par unite de u</span>
1840121 <span class="comment">% mini=-10.0;%general mini=uint16(min(values)-1 CHOOSE maxi-mini+1 EVEN</span>
1850122 <span class="comment">% maxi=9.0; % general maxi=uint16(max(values))+1</span>
1860123 <span class="comment">%nint=double(maxi-mini+1); % nombre d'intervals entiers EVEN!</span>
1870124 siz=size(aa);
1880125 nint=siz(2);
1890126 maxi=mini+nint-1;
1900127 npdf=nint*n;<span class="comment">% nbre total d'intervals à introduire dans la pdf: hist(u,npdf)</span>
1910128 <span class="comment">%simulation de pdf</span>
1920129 xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))];<span class="comment">% valeurs d'interpolation: we take n values in each integer interval</span>
1930130 <span class="comment">%histolin=exp(-(xfin-1).*(xfin-1)).*(2+cos(10*(xfin-1)));% simulation d'une pdf</span>
1940131 <span class="comment">%histo=log(histolin);</span>
1950132 <span class="comment">%histo=sin(2*pi*xfin);</span>
1960133 <span class="comment">%histextract=(reshape(histo,n,nint));</span>
1970134 <span class="comment">%aa=sum(histextract)/n %integral of the pdf on each integer bin</span>
1980135 IP=[0 diff(aa)];
1990136 Irev=zeros(size(aa));
2000137 <span class="keyword">for</span> i=1:nint
2010138     Irev(i)=aa(end-i+1);
2020139 <span class="keyword">end</span>
2030140 IPrev=[0 diff(Irev)];
2040141
2050142 <span class="comment">%get the spline coelfficients a_d, using filter on the eigen vectors A,B,C</span>
2060143 Arev=filter(b1rev,a1rev,IPrev);
2070144 Brev=filter(b2rev,a2rev,IPrev);
2080145 C=filter(b3,a3,IP);
2090146 D=filter(b4,a4,IP);
2100147 A=zeros(size(Arev));
2110148 B=zeros(size(Brev));
2120149 <span class="keyword">for</span> i=1:nint
2130150     A(i)=Arev(end-i+1);
2140151     B(i)=Brev(end-i+1);
2150152 <span class="keyword">end</span>
2160153 <span class="comment">%Matr=V*[A;B;C;D];</span>
2170154 bb=V*[A;B;C;D];
2180155 <span class="comment">%b=Matr(1,:);</span>
2190156 <span class="comment">%c=Matr(2,:);</span>
2200157 <span class="comment">%d=Matr(3,:);</span>
2210158 <span class="comment">%e=Matr(4,:);</span>
2220159 <span class="comment">%a=aa;</span>
2230160
2240161 <span class="comment">%calculate the interpolation using the spline coefficients a-d</span>
2250162 <span class="comment">%xextract=(reshape(xfin,n,nint));%</span>
2260163 chi=xfin+1/(2*n)-min(xfin)-double(int16(xfin+(1/(2*n))-min(xfin)))-0.5;<span class="comment">% decimal part</span>
2270164 chi2=chi.*chi;
2280165 chi3=chi2.*chi;
2290166 chi4=chi3.*chi;
2300167 avec=reshape(ones(n,1)*aa,1,n*nint);
2310168 bvec=reshape(ones(n,1)*bb(1,:),1,n*nint);
2320169 cvec=reshape(ones(n,1)*bb(2,:),1,n*nint);
2330170 dvec=reshape(ones(n,1)*bb(3,:),1,n*nint);
2340171 evec=reshape(ones(n,1)*bb(4,:),1,n*nint);
2350172 <a href="#_sub3" class="code" title="subfunction histx= histsmooth(chi,cc)">histsmooth</a>=avec+bvec.*chi+cvec.*chi2+dvec.*chi3+evec.*chi4;
2360173
2370174 <span class="comment">%-------------------------------------------------------</span>
2380175 <span class="comment">% --- determine the interpolated histogram at points chi from the spline ceff cc.</span>
2390176 <span class="comment">%-------------------------------------------------</span>
2400177 <a name="_sub3" href="#_subfunctions" class="code">function histx= histsmooth(chi,cc)</a>
2410178 <span class="comment">% provides the value of the interpolated histogram at values chi=x-i</span>
2420179 <span class="comment">%(difference with the mnearest integer)</span>
2430180 <span class="comment">% cc(5,size(chi)) is the set of spline coefficients obtained by splinhist</span>
2440181 chi2=chi.*chi;
2450182 chi3=chi2.*chi;
2460183 chi4=chi3.*chi;
2470184 histx=exp(cc(1,:)+cc(2,:).*chi+cc(3,:).*chi2+cc(4,:).*chi3+cc(5,:).*chi4);
2480185
2490186 <span class="comment">%-------------------------------------------------------</span>
2500187 <span class="comment">% --- determine the derivative p'/p of the interpolated histogram at points chi from the spline ceff cc.</span>
2510188 <span class="comment">%-------------------------------------------------</span>
2520189 <a name="_sub4" href="#_subfunctions" class="code">function histder= histder(chi,cc)</a>
2530190 <span class="comment">% provides the logarithmique derivative p'/p of the interpolated histogram</span>
2540191 <span class="comment">%at values chi=x-i</span>
2550192 <span class="comment">%(difference with the nearest integer)</span>
2560193 <span class="comment">% cc(5,size(chi)) is the set of spline coefficients obtained by splinhist</span>
2570194 chi2=chi.*chi;
2580195 chi3=chi2.*chi;
2590196 chi4=chi3.*chi;
2600197 <a href="#_sub4" class="code" title="subfunction histder= histder(chi,cc)">histder</a>=cc(2,:)+2*cc(3,:).*chi+3*cc(4,:).*chi2+4*cc(5,:).*chi3;</pre></div>
261<hr><address>Generated on Fri 13-Nov-2009 11:17:03 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
262</body>
263</html>
Note: See TracBrowser for help on using the repository browser.