Интерполяция

Можно создать кубический сплайновый интерполятор, соответствующий косинусной функции на следующих площадках x, с использованием csapi команда.

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-периодической. Насколько хорошо в этом отношении работает наш кубический сплайн-интерполятор? Один из способов проверки - вычислить разницу в первой производной на двух конечных точках.

diff( fnval( fnder(cs), [0 2*pi] ) )

ans = -0.1375

Для обеспечения периодичности использования csape вместо csapi.

csp = csape( x, y, 'periodic' );
hold on
fnplt(csp,'g')
hold off

Теперь чек дает

diff( fnval( fnder(csp), [0 2*pi] ) )

ans = -2.2806e-17

Даже вторая производная теперь совпадает в конечных точках.

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])

интерполяция дает интерполяцию wiggly, показанную ниже синим цветом.

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

Аппроксимация методом наименьших квадратов

Можно также использовать аппроксимацию наименьших квадратов к шумным данным сплайном с несколькими степенями свободы.

Например, можно попробовать кубический сплайн всего с четырьмя участками.

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обычно необходимо указать определенное сплайновое пространство. Это делается путем задания последовательности узлов и порядка, и это может быть немного проблемой. Однако при выполнении интерполяции сплайна в x,y данные с использованием сплайна порядка k, вы можете использовать функцию 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

Данные с привязкой к сетке

Все команды интерполяции и аппроксимации сплайна на панели инструментов фитинга кривой также могут обрабатывать данные с сеткой в любом числе переменных.

Например, вот бикубический сплайн, интерполированный с функцией мексиканской шляпы.

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) может работать с векторными сплайнами. Это также упрощает работу с параметрическими кривыми.

Вот, например, аппроксимация до бесконечности, полученная подведением кубической сплайновой кривой через точки, отмеченные на следующем рисунке.

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

Вот та же кривая, но с движением в третьем измерении.

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

Две половины кривой отображаются различными цветами, а начало координат обозначается как средство визуализации этой двукрылой кривой пространства.

Поверхности

Двумерные сплайны тензор-произведение со значениями в 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