В этом примере показано, как использовать команды из Curve Fitting Toolbox™ для подгонки сплайна к тестовым данным на титан с ручным и автоматическим выбором узлов.
Вот некоторые данные, которые регистрируют определенное свойство титана, измеренное как функция от температуры. Мы будем использовать его, чтобы проиллюстрировать некоторые проблемы с сплайн интерполяцией.
[xx,yy] = titanium;
График данных показывает довольно резкий пик.
plot(xx,yy,'bx');
frame = [-10 10 -.1 .3]+[min(xx),max(xx),min(yy),max(yy)];
axis(frame);
Мы выбираем несколько точек данных из этих несколько приблизительных данных, так как мы хотим интерполировать. Вот изображение данных с отмеченными выбранными точками данных.
pick = [1 5 11 21 27 29 31 33 35 40 45 49]; tau = xx(pick); y = yy(pick); hold on plot(tau,y,'ro'); hold off
Поскольку сплайн порядка k с n + k узлами имеет n степеней свободы, и у нас есть 12 точек данных, для подгонки со сплайном четвертого порядка требуется 12 + 4 = 16 узлов. Кроме того, эта последовательность узлов t должна быть такой, чтобы i-й сайт данных лежал в поддержке i-го B-сплайна. Мы достигаем этого, используя сайты данных как узлы, но добавляем два простых узла на любом конце.
dl = tau(2) - tau(1); dr = tau(end) - tau(end-1); t = [tau(1)-dl*[2 1] tau tau(end)+dr*[1 2]]; % construct the knot sequence plot(tau,y,'ro'); hold on axis(frame+[-2*dl 2*dr 0 0]) plot(t,repmat(frame(3)+.03,size(t)),'kx') hold off legend({'Data Values' 'Knots'},'location','NW')
Мы используем эту последовательность узлов, чтобы создать интерполирующий кубический сплайн.
sp = spapi(t,tau,y);
Теперь для графика. Поскольку мы не заботимся о части сплайна вне интервала данных, мы ограничиваем график этим интервалом.
plot(tau,y,'ro') axis(frame) hold on fnplt(sp,[tau(1) tau(end)], 'k') hold off
Более близкий взгляд на левую часть сплайна подгонки показывает некоторые волнистости.
xxx = linspace(tau(1),tau(5),41); plot(xxx, fnval(sp, xxx), 'k', tau, y, 'ro'); axis([tau(1) tau(5) 0.6 1.2]);
Неразумный отбойник в первом интервале проистекает из того, что наш сплайн плавно переходит в нуль в своём первом узле. Чтобы увидеть это, вот изображение всего сплайна, вместе с его последовательностью узлов и точками данных.
fnplt(sp,'k'); hold on plot(tau,y,'ro', t,repmat(.1,size(t)),'kx'); hold off legend({'Spline Interpolant' 'Data Values' 'Knots'},'location','NW')
Вот простой способ обеспечить более разумное граничное поведение. Мы добавляем еще две точки данных за пределами заданного интервала данных и выбираем в качестве наших данных значения прямой линии через первые две точки данных.
tt = [tau(1)-[4 3 2 1]*dl tau tau(end)+[1 2 3 4]*dr]; xx = [tau(1)-[2 1]*dl tau tau(end)+[1 2]*dr]; yy = [y(1)-[2 1]*(y(2)-y(1)) y y(end)+[1 2]*(y(end)-y(end-1))]; sp2 = spapi(tt,xx,yy); plot(tau,y,'ro', xx([1 2 end-1 end]),yy([1 2 end-1 end]),'bo'); axis(frame+[-2*dl 2*dr 0 0]); hold on fnplt(sp2,'b',tau([1 end])) hold off legend({'Original Data' 'Data Added for End Conditions' ... 'Fit with Added Data'},'location','NW')
Вот сравнение двух сплайнов подгонки, чтобы показать уменьшение волнообразования в первом и последнем интервале.
hold on fnplt(sp,'k',tau([1 end])) hold off legend({'Original Data' 'Data Added for End Conditions' ... 'Fit with Added Data' 'Original Fit'},'location','NW')
Наконец, вот более пристальный взгляд на первые четыре интервала данных, который более четко показывает уменьшение волнистости около левого конца.
plot(tau,y,'ro',xxx,fnval(sp2,xxx),'b',xxx,fnval(sp,xxx),'k'); axis([tau(1) tau(5) .6 1.2]); legend({'Original Data' 'Fit with Added Data' ... 'Original Fit'},'location','NW')
Если вся эта деталь отключает вас, пусть Curve Fitting Toolbox выберет узлы для вас. Задайте требуемый порядок интерполяции в качестве первого входного параметра в команду сплайн интерполяции spapi, а не последовательность узлов.
autosp = spapi(4, tau, y); knots = fnbrk(autosp,'knots'); plot(tau, y, 'ro') hold on fnplt(autosp,'g') plot(knots, repmat(.5,size(knots)),'gx') hold off legend({'Data Values' 'Fit With Knots Chosen by SPAPI' ... 'Knots Chosen by SPAPI'}, 'location','NW')
Ниже приведен результат гораздо лучшего выбора узла, полученного смещением узла в 842 слегка вправо и узла в 985 слева.
knots([7 12]) = [851, 971]; adjsp = spapi(knots, tau, y); hold on fnplt(adjsp,'r',2) plot(knots, repmat(.54,size(knots)),'rx') hold off legend({'Data Values' 'Fit With Knots Chosen by SPAPI' ... 'Knots Chosen by SPAPI' 'Fit With Knots Adjusted' ... 'Adjusted Knots'}, 'location','NW')
Или же просто попробуйте стандартную кубическую сплайн интерполяцию, поставляемую csapi. Это сводится к несколько другому выбору узлов.
autocs = csapi(tau, y); plot(tau, y, 'ro') hold on fnplt(autocs,'c') hold off
С такими быстро изменяющимися данными трудно договориться между всеми разумными интерполянтами, даже если каждая из них является кубическим сплайном. На графике ниже показаны все пять интерполяций, для сравнения.
plot(tau, y, 'ro') hold on fnplt(sp,'k',tau([1 end])) % black: original fnplt(sp2,'b',tau([1 end])) % blue: with special end conditions fnplt(autosp,'g') % green: automatic knot choice by SPAPI fnplt(autocs,'c') % cyan: automatic knot choice by CSAPI fnplt(adjsp,'r',2) % red: knot choice by SPAPI slightly changed hold off legend({'Data Values' 'Original Fit' 'Special End Conditions' ... 'With Knots Chosen by SPAPI' 'With Knots Chosen by CSAPI' ... 'With Adjusted Knots'},'location','NW')
