В этом примере показано, как использовать команды из Toolbox™ Фитинг кривой (Curve Fitting) для подгонки сплайна к данным титановых испытаний с ручным и автоматическим выбором узлов.
Вот некоторые данные, которые записывают определенное свойство титана, измеренное как функция температуры. Мы используем его для иллюстрации некоторых проблем с интерполяцией сплайнов.
[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')
