前几天在进行数据仿真的时候,对于将表格离散数据转化成连续数据一直是一件十分棘手的事情，在网站上找了许多资源最后才找到可以利用二维线性插值的方法将数据进行转化。
是要将
m
×
n
m\times n
m×n的二维数据如下图所示进行拟合操作找出
C
=
f
(
M
a
,
α
)
C=f(Ma,\alpha)
C=f(Ma,α)的插值式。在这里我们假设对表格中的一个
2
×
2
2 \times 2
2×2的数据进行插值操作,其基本系数描述如下:
f
(
M
a
,
α
)
=
a
0
+
a
1
M
a
+
a
2
α
+
a
3
M
a
×
α
f(Ma,\alpha)=a_0+a_1Ma+a_2\alpha+a_3Ma\times\alpha
f(Ma,α)=a0+a1Ma+a2α+a3Ma×α
在已知四个坐标
(
x
1
,
y
1
,
f
1
)
,
(
x
2
,
y
2
,
f
2
)
,
(
x
3
,
y
3
,
f
3
)
,
(
x
4
,
y
4
,
f
4
)
(x_1,y_1,f_1),(x_2,y_2,f_2),(x_3,y_3,f_3),(x_4,y_4,f_4)
(x1,y1,f1),(x2,y2,f2),(x3,y3,f3),(x4,y4,f4)为插值点的情况下可以得到这四个坐标点区域内的拟合曲面的系数
(
α
0
,
α
1
,
α
2
,
α
3
)
T
(\alpha_0,\alpha_1,\alpha_2,\alpha_3)^T
(α0,α1,α2,α3)T:
(
1
x
1
y
1
x
1
y
1
1
x
2
y
2
x
2
y
2
1
x
3
y
3
x
3
y
3
1
x
4
y
4
x
4
y
4
)
(
a
0
a
1
a
2
a
3
)
=
(
f
1
f
2
f
3
f
4
)
A
x
=
b
→
x
=
A
−
1
b
\begin{pmatrix} 1 & x_1&y_1&x_1y_1\\1 & x_2&y_2&x_2y_2\\1 & x_3&y_3&x_3y_3\\1 & x_4&y_4&x_4y_4\\\end{pmatrix}\begin{pmatrix} a_0\\a_1\\a_2\\a_3\\\end{pmatrix}=\begin{pmatrix} f_1\\f_2\\f_3\\f_4\\\end{pmatrix}\\ Ax=b\rightarrow x = A^{-1}b
⎝⎜⎜⎛1111x1x2x3x4y1y2y3y4x1y1x2y2x3y3x4y4⎠⎟⎟⎞⎝⎜⎜⎛a0a1a2a3⎠⎟⎟⎞=⎝⎜⎜⎛f1f2f3f4⎠⎟⎟⎞Ax=b→x=A−1b
例如对下面的表格中的数据进行插值:
主函数:
clc,clear;close all;
X = 0.9;Y = 9*pi/180;
CX = [0.4171 0.3858 0.3779 0.3785 0.3787 0.3829 0.3855 0.4082 0.4947;...
0.4404 0.4086 0.4007 0.4015 0.4018 0.4062 0.4091 0.4321 0.5192;...
0.5219 0.4903 0.4827 0.4838 0.4846 0.4897 0.4934 0.5175 0.6073;...
0.6603 0.6290 0.6218 0.6234 0.6249 0.6310 0.6363 0.6621 0.7571;...
0.8534 0.8826 0.8160 0.8184 0.8209 0.8284 0.8358 0.8641 0.9672;...
1.1023 1.0723 1.0666 1.0700 1.0738 1.0835 1.0938 1.1254 1.2392];
%X,Y是想求的点处的位置
y = linear_cha(CX,X,Y);
%下面对表格中的数据进行线性插值函数
function y = linear_cha(CX,X,Y)
ma = 0.1:0.1:0.9;alpha = (0:2:10)*pi/180;
%寻找横坐标
for i = 1:8
if (ma(i)<X)&&(ma(i+1)>=X)
X_ = ma(i);
index_i = i;
break;
end
end
%寻找纵坐标
for i = 1:5
if (alpha(i)<Y)&&(alpha(i+1)>=Y)
Y_ = alpha(i);
index_j = i;
break;
end
end
x1 = ma(index_i);y1 = alpha(index_j);
x2 = ma(index_i);y2 = alpha(index_j + 1);
x3 = ma(index_i + 1);y3 = alpha(index_j + 1);
x4 = ma(index_i + 1);y4 = alpha(index_j);
%构造多项式矩阵求解系数
A = [1 x1 y1 x1*y1;...
1 x2 y2 x2*y2;...
1 x3 y3 x3*y3;...
1 x4 y4 x4*y4];
b = [CX(index_j,index_i);...
CX(index_j+1,index_i);...
CX(index_j+1,index_i+1);...
CX(index_j,index_i+1)];
x = A\b;
y = x'*[1;X;Y;X*Y];
end