基于SqlServer2008的MATLAB GUI应用程序开发（6）:Matlab Griddata 插值：Duplicat...

haobang008

源地址：http://blog.163.com/6_mao/blog/static/63271315201191603845663/

在matlab里面，当你用griddata函数插值的时候，有时候会出现这样的警告：
Warning: Duplicate x-y data points detected: using average of the z values.
大概意思就是，你的（x,y） 坐标值出现了重复的点，计算时会取z的平均值来计算。

我的数据大概有两百万对，一个个检查有没有重复可不是什么好办法。好在看见网上有人出了一个很聪明的方法：先建立一个新的矩阵A，使得A=x+sqrt(-1)*y. 这样就使得xy坐标变成了一个复数。接着借助unique()函数，这个函数可以判断出一个矩阵中所有不重复的值。举个例子，如果A=[1 1 2 3 4 5 5 6], unique(A)=[1 2 3 4 5 6]，因为A中出现了1到6. 这样利用unique(x+sqrt(-1)*y)就可以判断出数组中有没有重复的了。

可是我的数据经过这样的处理，发现并没有重复的xy值。这是怎么回事？

查询matlab的help文件，上面有这么一句话：

"this means that there are two or more data points where the input values (x, y) are the same or very close together."


这就是说，如果两个值差得很近，也许也会被当做一样来处理？关键是，这个小，究竟是多小？我们来看看griddata()的源代码是怎么说。这一段如下:
____________________________________________________________________________
% Sort x and y so duplicate points can be averaged before passing to delaunay need x,y and z to be column vectors
sz = numel(x);
x = reshape(x,sz,1);y = reshape(y,sz,1);z = reshape(z,sz,1);sxyz = sortrows([x y z],[2 1]);
x = sxyz(:,1);y = sxyz(:,2);z = sxyz(:,3);
myepsx = eps(0.5 * (max(x) - min(x)))^(1/3);
myepsy = eps(0.5 * (max(y) - min(y)))^(1/3);
ind = [0; ((abs(diff(y)) < myepsy) & (abs(diff(x)) < myepsx)); 0];
if sum(ind) > 0
  warning('MATLAB:griddata:DuplicateDataPoints',['Duplicate x-y data points detected: using average of the z values.']);

______________________________________________________________________________
前面的浅灰色的代码没什么意义，主要是转换一下矩阵的格式以方便计算。下面我们来看粗体的部分，含义如下：

sxyz = sortrows([x y z],[2 1]);
将需要插值的（x,y,z）放到一个大的矩阵中，x,y,z分别一列。然后以y为主导，对y从小到大排列，对应的x和z的值也跟着移动。


myepsx = eps(0.5 * (max(x) - min(x)))^(1/3);
myepsy = eps(0.5 * (max(y) - min(y)))^(1/3);

这两行代码，是分别计算x和y的表示的精度值，也可以理解为在你的具体应用中，x和y的值的最小精度是多少。这里计算的方法是计算二分之一x的阈值的开三次方的精度值，应该也可以有别的方法。

ind = [0; ((abs(diff(y)) < myepsy) & (abs(diff(x)) < myepsx)); 0];
diff(y)的意思是[y2-y1,y3-y2,......,yn-y(n-1)]形成的一个新的矩阵，如果这个新的矩阵的每一个元素的值小于之前计算的y的最小精度，那么就认为这个值是0，或者说认为这两个值相等。如果小于的条件是成立的，就返回1（True）；反之则返回0（False）。后面的那个是一样的，就是x的值的对比。
那么我们想想，两个对应的列，如果y的差值还小于精度值（就是认为两个值相等，结果返回1），如果这时候同时X的差值还小于精度值（就是x的两个值也相等，结果返回1），那个1&1=1，结果就返回1.

if sum(ind) > 0
只要结果中有一个1，就说明有一对(x,y)是相等的（或者非常接近），那个总和就必定大于1，就会返回这个警告。

比如说，经过整理以后的结果如下
X=[1 3.001 5 8 0 2 3 3.002 7]
Y=[1 1.001 2 3 5 6 7 7.001 8]
假设精度是0.5（我没有计算，随手举的），那么就意味着在你的计算中最高精度是0.5，低于0.5都可能是没有意义的误差。虽然从这上面看X，Y没有相同的配对，可是如果用上面的计算：
diff(X)=[2.001 1.999 3 -8 2 1 0.002 3.998]
diff(Y)=[0.001 0.999 1  2 1 1 0.001 0.999]
首先看Y，虽然没有0，但是有两个0.001，比0.5小，就可以认为是一样的值，返回两个真值。再检查这两个真值的X的值，其中有一个2.001，大于0.5，说明不同；另外一个，X差值是0.002，小于0.5，说明精度已经低于最小的计算值了，应该认为是相等的。我们回到原来的数据，[3 7]和[3.002 7.001]，如果EPS（精度）是大于0.002的，说明后面的尾数很有可能只是超过精度的计算误差，是没有意义的，这个时候做插值，就会返回之前的那个错误，说你有重复的xy值。

不过，总的来说，这个报错并不影响你的结果；他都说了会取z的平均值计算。如果你真的觉得平均值不准想要取特定的值，可以用unique函数找到特定的行列（或者用griddata里面的这个判断也行），然后更改你的z值达到你的目的。

嗯，我觉得这样的处理还是合理的，不能只一根筋地看是不是数据都一样，我们都知道很多在计算机里面的数据，超过精度的后面的若干位都是没有意义的数据，比如说double精度的计算最小数应该是2.22*e-16吧，如果你的数字比这个还小，或者小数点后位数比这个还多的，那后面的计算都完全没有意义，都是扯淡。就跟我前两天在bbs上看到说matlab算tan(10^100)，出来的结果根本不准，大概就是这个关系。不过tan(10^100)到底怎么算才是准的，那就是另外一个问题了。当然这也提醒了我，如果进行这种太大的数或者位数太多的数的计算，一定要留个心眼，会不会得到一个无效的结果。

kmdqlb

好东西，共分享。。。。。。