スプラインの作成方法
この例では、Curve Fitting Toolbox™ でスプライン関数を使用してさまざまな方法でスプラインを作成する方法を説明します。
内挿
csapi
コマンドを使用して、次のサイト x
で余弦関数に一致する "3 次スプライン内挿" を作成できます。
x = 2*pi*[0 1 .1:.2:.9]; y = cos(x); cs = csapi(x,y);
これで、fnplt
を使用して内挿スプラインを表示できます。
fnplt(cs,2); axis([-1 7 -1.2 1.2]) hold on plot(x,y,'o') hold off
内挿のチェック
余弦関数は 2*pi 周期です。その点について、この 3 次スプライン内挿はどの程度の近似になっているでしょうか。これを確認する 1 つの方法は、2 つの端点における 1 階微分の差異を計算することです。
diff( fnval( fnder(cs), [0 2*pi] ) )
ans = -0.1375
周期性を適用するには、csapi
の代わりに csape
を使用します。
csp = csape( x, y, 'periodic' ); hold on fnplt(csp,'g') hold off
この確認により、次が得られます。
diff( fnval( fnder(csp), [0 2*pi] ) )
ans = -2.2806e-17
2 階微分でも端点で一致しています。
diff( fnval( fnder(csp, 2), [0 2*pi] ) )
ans = -2.2204e-16
同じデータへの "区分的線形内挿" は spapi
で使用できます。ここでは、以前のプロットにこの内挿 (赤色) を追加します。
pl = spapi(2, x, y); hold on fnplt(pl, 'r', 2) hold off
平滑化
データのノイズが多い場合、通常は内挿ではなく近似を行う必要があります。たとえば、次のようなデータがあるとします。
x = linspace(0,2*pi,51);
noisy_y = cos(x) + .2*(rand(size(x))-.5);
plot(x,noisy_y,'x')
axis([-1 7 -1.2 1.2])
内挿すると、次の青色で示される、うねりのある内挿が得られます。
hold on fnplt( csapi(x, noisy_y) ) hold off
それに対し、次のような適切な許容誤差で平滑化を行うと
tol = (.05)^2*(2*pi)
tol = 0.0157
次の赤色で示す平滑化された近似が得られます。
hold on fnplt( spaps(x, noisy_y, tol), 'r', 2 ) hold off
この近似は区間の両端付近で非常に悪化していて、周期性がありません。周期性を適用するには、周期的に拡張されたデータに近似を行ってから、その近似を元の区間に制限します。
noisy_y([1 end]) = mean( noisy_y([1 end]) ); lx = length(x); lx2 = round(lx/2); range = [lx2:lx 2:lx 2:lx2]; sps = spaps([x(lx2:lx)-2*pi x(2:lx) x(2:lx2)+2*pi],noisy_y(range),2*tol);
これにより、ほぼ周期的な近似が得られます (黒で表示)。
hold on fnplt(sps, [0 2*pi], 'k', 2) hold off
最小二乗近似
あるいは、自由度がほとんどないスプラインで、ノイズを含むデータに最小二乗近似を使用することもできます。
たとえば、4 つの区分のみから成る 3 次スプラインを試してみるとします。
spl2 = spap2(4, 4, x, noisy_y); fnplt(spl2,'b',2); axis([-1 7 -1.2 1.2]) hold on plot(x,noisy_y,'x') hold off
節点の選択
spapi
または spap2
を使用する場合、通常は特定のスプライン空間を指定する必要があります。このために、"節点シーケンス" および "次数" を指定しますが、これが少し問題になる場合があります。ただし、次数 k
のスプラインを使用して x,y
データにスプライン内挿を行う場合は、関数 optknt
を使用して適切な節点シーケンスを指定できます (次の例を参照)。
k = 5; % order 5, i.e., we are working with quartic splines x = 2*pi*sort([0 1 rand(1,10)]); y = cos(x); sp = spapi( optknt(x,k), x, y ); fnplt(sp,2,'g'); hold on plot(x,y,'o') hold off axis([-1 7 -1.1 1.1])
最小二乗近似を行う場合は、現在の近似とともに newknt
を用いることでより適切な節点の選択を決定できます。たとえば、次の赤色でプロットされた指数関数の近似は、その誤差からわかるように、あまり良いものではありません。
x = linspace(0,10,101); y = exp(x); sp0 = spap2( augknt(0:2:10,4), 4, x, y ); plot(x,y-fnval(sp0,x),'r','LineWidth',2)
しかし、この元の近似を利用して、節点の数は "同じ" でも、分散が改善された別の近似を作成できます。その誤差を黒でプロットしています。
sp1 = spap2( newknt(sp0), 4, x, y ); hold on plot(x,y-fnval(sp1,x),'k','LineWidth',2) hold off
グリッド データ
Curve Fitting Toolbox のスプライン内挿および近似のコマンドはすべて、任意の数の変数においてグリッド データを扱うこともできます。
例として、次に Mexican Hat 関数に対する双三次スプライン内挿を示します。
x =.0001+(-4:.2:4); y = -3:.2:3; [yy,xx] = meshgrid(y,x); r = pi*sqrt(xx.^2+yy.^2); z = sin(r)./r; bcs = csapi({x,y}, z); fnplt(bcs) axis([-5 5 -5 5 -.5 1])
次に、同じグリッド上の同じ関数で、ノイズが含まれる値への "最小二乗" 近似を示します。
knotsx = augknt(linspace(x(1), x(end), 21), 4); knotsy = augknt(linspace(y(1), y(end), 15), 4); bsp2 = spap2({knotsx,knotsy},[4 4], {x,y},z+.02*(rand(size(z))-.5)); fnplt(bsp2) axis([-5 5 -5 5 -.5 1])
曲線
Curve Fitting Toolbox で "ベクトル値の" スプラインを扱うことができるため、グリッド データを容易に処理できます。また、そのおかげでパラメトリック曲線も簡単に操作できます。
例として、以下の Figure に付けられた点を 3 次スプライン曲線で処理して得られた無限大の近似を示します。
t = 0:8; xy = [0 0;1 1; 1.7 0;1 -1;0 0; -1 1; -1.7 0; -1 -1; 0 0].'; infty = csape(t, xy, 'periodic'); fnplt(infty, 2) axis([-2 2 -1.1 1.1]) hold on plot(xy(1,:),xy(2,:),'o') hold off
次に、同じ曲線を 3 次元の運きで表したものを示します。
roller = csape( t , [ xy ;0 1/2 1 1/2 0 1/2 1 1/2 0], 'periodic'); fnplt( roller , 2, [0 4],'b' ) hold on fnplt( roller, 2, [4 8], 'r') plot3(0,0,0,'o') hold off
2 枚の羽のような空間曲線を可視化しやすくするために曲線を二等分して別々の色でプロットし、原点に印を付けてあります。
曲面
R^3 の値をもつ二変量テンソル積スプラインからは曲面が得られます。例として、次にトーラスの適切な近似を示します。
x = 0:4; y = -2:2; R = 4; r = 2; v = zeros(3,5,5); v(3,:,:) = [0 (R-r)/2 0 (r-R)/2 0].'*[1 1 1 1 1]; v(2,:,:) = [R (r+R)/2 r (r+R)/2 R].'*[0 1 0 -1 0]; v(1,:,:) = [R (r+R)/2 r (r+R)/2 R].'*[1 0 -1 0 1]; dough0 = csape({x,y},v,'periodic'); fnplt(dough0) axis equal, axis off
次に、その曲面に対する法線から成る冠を示します。
nx = 43; xy = [ones(1,nx); linspace(2,-2,nx)]; points = fnval(dough0,xy)'; ders = fnval(fndir(dough0,eye(2)),xy); normals = cross(ders(4:6,:),ders(1:3,:)); normals = (normals./repmat(sqrt(sum(normals.*normals)),3,1))'; pn = [points;points+normals]; hold on for j=1:nx plot3(pn([j,j+nx],1),pn([j,j+nx],2),pn([j,j+nx],3)) end hold off
最後に、(x,y) 面への投影を示します。
fnplt(fncmb(dough0, [1 0 0; 0 1 0]))
axis([-5.25 5.25 -4.14 4.14]), axis off
散布データ
平面のグリッドなしのデータ サイトで与えられた値に内挿することもできます。例として、単位正方形を単位円板に滑らかにマッピングする作業を考えてみます。データ値を作成して円でマークし、それに対応するデータ サイトを x でマークします。各データ サイトをそれに関連する値に矢印で接続します。
n = 64; t = linspace(0,2*pi,n+1); t(end) = []; values = [cos(t); sin(t)]; plot(values(1,:),values(2,:),'or') axis equal, axis off sites = values./repmat(max(abs(values)),2,1); hold on plot(sites(1,:),sites(2,:),'xk') quiver(sites(1,:),sites(2,:), ... values(1,:)-sites(1,:), values(2,:)-sites(2,:)) hold off
次に、tpaps
を使用して二変量で内挿したベクトル値の薄板スプラインを作成します。
st = tpaps(sites, values, 1);
fnplt
によるプロットが示すように、このスプラインは実際に単位正方形を単位円板に滑らかに (近似的に) マップします。このプロットは、st
内のスプライン マップにおける等間隔の正方形グリッドのイメージを示しています。
hold on fnplt(st) hold off