关于matlab的v4插值

matlab的v4插值

griddata(x,y,z,linspace(min(x),max(x))’,linspace(min(y),max(y))

1. 根据帮助文件，v4插值为双调和样条插值（biharmonic spline interpolation）
   

2. 关于双调和样条插值，参考的此网址中的文献
   http://www.fon.hum.uva.nl/praat/manual/biharmonic_spline_interpolation.html

3. 文献（doi：10.1007/s11004-011-9346-5）中给出的方法：
   

4. 根据文献所给方法得到的matlab代码（参数仅适用于此案例，具体情况应按文献修改）

   x=[27.5	77.5	111.25	0	0	0	-27.5	-77.5	-111.25	0	0	0];
    y=[0	0	0	27.5	77.5	111.25	0	0	0	-27.5	-77.5	-111.25];
    z=[0.256700000000000,0.0720900000000000,0.0357500000000000,0.246100000000000,0.273500000000000,0.0820800000000000,0.257700000000000,0.0710900000000000,0.0357500000000000,0.236100000000000,0.283500000000000,0.0820800000000000];
   %[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y)),'v4');%插值
%Z(find(Z<0))=0;
%以下内容为根据文献得出代码，作用相当于上方两行代码
    X=linspace(min(x),max(x));%将最小值与最大值划分为100块
    Y=linspace(min(y),max(y));
 %以下求权重矩阵   
i=1;
while i~=13
    j=1;
    while j~=13
        r(i,j)=((x(i)-x(j))^2+(y(i)-y(j))^2)^0.5;
        if (r(i,j)==0)
            d(i,j)=0;
        else
            d(i,j)=(log(r(i,j))-1)*r(i,j)^2;%格林函数矩阵，即文献中d(i,j)，即G
        end
        j=j+1;
    end
    i=i+1;
end
w=inv(d)*z';%权重矩阵

%以下求XY的100*100矩阵
n=1;
while n~=101
    m=1;
    while m~=101
        X(m,n)=X(1,n);
        
        m=m+1;
    end
    n=n+1;
end

Y=X';

%以下求Z矩阵
n=1;
while n~=101
    m=1;
    while m~=101
        i=1;
        while i~=13
            v(1,i)=((X(m,n)-x(i))^2+(Y(m,n)-y(i))^2)^0.5;%即文献中r(i,j)
            if v(1,i)==0
                u(1,i)=0;
            else
                u(1,i)=(log(v(1,i))-1)*v(1,i)^2;%即文献中d(i,j)，即Gp
            end
            i=i+1;
        end
        Z(m,n)=u*w;
        m=m+1;
    end
    n=n+1;
end
Z(find(Z<0))=0;%将Z中小于0的数取0
figure,surf(X,Y,Z); 

5. v4插值代码（此结果与标题4的计算结果相同）

x=[27.5	77.5	111.25	0	0	0	-27.5	-77.5	-111.25	0	0	0];
    y=[0	0	0	27.5	77.5	111.25	0	0	0	-27.5	-77.5	-111.25];
    z=[0.256700000000000,0.0720900000000000,0.0357500000000000,0.246100000000000,0.273500000000000,0.0820800000000000,0.257700000000000,0.0710900000000000,0.0357500000000000,0.236100000000000,0.283500000000000,0.0820800000000000];
   [X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y)),'v4');%插值
Z(find(Z<0))=0;
figure,surf(X,Y,Z); 

6.本人matlab初学者，欢迎批评指正