Главная | Обратная связь | Поможем написать вашу работу!

Интерполяция и аппроксимация данных

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

Результат аппроксимации может не проходить через узловые точки. Напротив, задача интерполяции — найти данные в окрестности узловых точек. Для этого используются подходящие функции, значения которых в узловых точках совпадают с координатами этих точек. Например, при линейной интерполяции зависимости у(х) узловые точки соединяются друг с другом отрезками прямых и считается, что искомые промежуточные точки расположены на этих отрезках.

Для повышения точности интерполяции применяют параболы (квадратичная интерполяция) или полиномы более высокой степени (полиномиальная интерполяция). Для обработки данных MATLAB использует различные функции интерполяции и аппроксимации данных. Набор таких функций вместе с несколькими вспомогательными функциями описан в этом разделе.

Полиномиальная регрессия

Одна из наиболее известных аппроксимаций — полиномиальная. В системе MATLAB определены функции аппроксимации данных полиномами по методу наименьших квадратов — полиномиальной регрессии. Это выполняет функция, приведенная ниже:

ñpolyfit(x.y.n) — возвращает вектор коэффициентов полинома р(х) степени п, который с наименьшей среднеквадратичной погрешностью аппроксимирует функцию у(х). Результатом является вектор-строка длиной n+1, содержащий коэффициенты полинома в порядке уменьшения степеней х и у равно n+1, то реализуется обычная полиномиальная аппроксимация, при которой график полинома точно проходит через узловые точки с координатами (х.у), хранящиеся в векторах х и у. В противном случае точного совпадения графика с узловыми точками не наблюдается;

ñ[p.S] = polyflt(x.y.n) — возвращает коэффициенты полинома р и структуру S для использования вместе с функцией polyval с целью оценивания или предсказания погрешности;

ñ[p.S] = polyf1t(x,y,n,mu) возвращает коэффициенты полинома р и структуру S для использования вместе с функцией polyval с целью оценивания или предска-зния погрешности, но так, что присходит центрирование (нормирование) и масштабирование х, xnorm = (х - mu(l))/mu(2), где mu(l) = mean(x) и mu(2) = std(x). Центрирование и масштабирование не только улучшают свойства степенного многочлена, получаемого при помощи polyval, но и значительно повышают качественные характеристики самого алгоритма аппроксимации.

Рис. 6.10. Пример использования функции polyfit

Пример (полиномиальная регрессия для функции s

» х=(-3:0.2:3)':

y=sin(x);

p=polyflt(x.y,3)

р =

-0.0953 0.0000 0.8651 -0.0000

»x=(-4:0.2:4)';y=sin(x);

» f=polyval(p,x);plot(x.y.'o',x,f)

Рис. 6.10, построенный в этом примере, дает наглядное представление о точности полиномиальной аппроксимации. Следует помнить, что она достаточно точна в небольших окрестностях от точки х = 0, но может иметь большие погрешности за их пределами или в промежутках между узловыми точками.

График аппроксимирующего полинома третьей степени на рис. 6.10 показан сплошной линией, а точки исходной зависимости обозначены кружками. К сожалению, при степени полинома свыше 5 погрешность полиномиальной регрессии (и аппроксимации) сильно возрастает и ее применение без центрирования и масштабирования становится рискованным. Обратите внимание на то, что при полиномиальной регрессии узловые точки не ложатся точно на график полинома, поскольку их приближение к нему является наилучшим в смысле минимального среднеквадратического отклонения. Об этом уже говорилось.

Интерполяция периодических функций рядом Фурье

Под интерполяцией обычно подразумевают вычисление значений функции f(x) в промежутках между узловыми точками. Линейная, квадратичная и полиномиальная интерполяция реализуются при полиномиальной аппроксимации. А вот для периодических (и особенно для гладких периодических) функций хорошие результаты может дать их интерполяция тригонометрическим рядом Фурье. Для этого используется следующая функция:

ñinterpft(x.n) — возвращает вектор у, содержащий значения периодической функции, определенные в п равномерно расположенных точках. Если length(x)=rr; и х имеет интервал дискретизации dx, то интервал дискретизации для у составляет dy=dx*m/n, причем п не может быть меньше, чем т. Если X — матрица, interpft оперирует столбцами X, возвращая матрицу Y с таким же числом столбцов, как и у X, но с п строками. Функция y=interpft(x.n.dim) работает либо со строками, либо со столбцами в зависимости от значения параметра dim.

Рис. 6.11. Пример использования функции interpft

Пример:

» x=0:10:y-sin(x).^3:

» xl-0:0.1:10;yl-interpft(y,101);

» x2=0:0.01:10:y2-sin(x2).^3;

» plot(xl,yl, '--').hold on.plotCx.y. 'o'.x2.y2)

Рис. 6.11 иллюстрирует эффективность данного вида интерполяции на примере функции sin(x).^3, которая представляет собой сильно искаженную синусоиду.

Исходная функция на рис. 6.11 представлена сплошной линией с кружками, а интерполирующая функция — штрих-пунктирной линией.

Интерполяция на неравномерной сетке

Для интерполяции на неравномерной сетке используется функция griddata:

ñZI = griddata(x.y.z.XI.YI) — преобразует поверхность вида z = f(x. у), которая определяется векторами (x.y.z) с (обычно) неравномерно распределенными элементами. Функция griddata аппроксимирует эту поверхность в точках, определенных векторами (XI.YI) в виде значений ZI. Поверхность всегда проходит через заданные точки. XI и YI обычно формируют однородную сетку (созданную с помощью функции meshgrid).

XI может быть вектором-строкой, в этом случае он определяет матрицу с постоянными столбцами. Точно так же YI может быть вектором-столбцом, тогда он определяет матрицу с постоянными строками.

ñ[XI.YI.ZI] = griddata(x,y,z,xi,yi) — возвращает аппроксимирующую матрицу ZI, как описано выше, а также возвращает матрицы XI и YI, сформированные из вектора-столбца xi и вектора-строки yi. Последние аналогичны матрицам, возвращаемым функцией meshgrid;

ñ[...] = griddata (....method) — использует определенный метод интерполяции:

ñ'nearest' — ступенчатая интерполяция;

ñ'linear' — линейная интерполяция (принята по умолчанию);

ñ'cubic' — кубическая интерполяция;

ñ' v4 ' — метод, используемый в МATLAB 4.

Метод определяет тип аппроксимирующей поверхности. Метод 'cubic1 формирует гладкие поверхности, в то время как 'linear' и 'nearest' имеют разрывы первых и нулевых производных соответственно. Все методы, за исключением v4, основаны на триангуляции Делоне. Метод ' v4 ' включен для обеспечения совместимости с версией 4 системы MATLAB. Пример:

» x=rand(120.1)*4-2;y=rand(120.1)*4-2;

z=x.*y.*exp(-x.^2-y.^2);

» t=-2:0.1:2:[X,Y]=meshgrid(t,t);

Z=griddata(x.y.z.X.Y):

» mesh(X.Y.Z).hold on.plot3(x.y,z, 'ok')

Функции griddataS и griddatan работают аналогично griddata, но для для трехмерного и n-мерного случая — с использованием алгоритма qhul 1. Используются, в частности, при трехмерной и n-мерной триангуляции.

Рис. 6.12 иллюстрирует применение функции griddata.

Рис 6.12. Пример использования функции griddata

Одномерная табличная интерполяция

В ряде случаев очень удобна сплайновая интерполяция и аппроксимация таблично заданных функций. При ней промежуточные точки ищутся по отрезкам полиномов третьей степени — это кубическая сплайновая интерполяция. При этом обычно такие полиномы вычисляются так, чтобы не только их значения совпадали с координатами узловых точек, но также чтобы в узловых точках были непрерывны производные первого и второго порядков. Такое поведение характерно для гибкой линейки, закрепленной в узловых точках, откуда и происходит название spline (сплайн) для этого вида интерполяции (аппроксимации). Для одномерной табличной интерполяции используется функция interpl:

ñyi = Interpl(x.Y.xi) — возвращает вектор yi, содержащий элементы, соответствующие элементам xi и полученные интерполяцией векторов х и Y. Вектор х определяет точки, в которых задано значение Y. Если Y — матрица, то интерполяция выполняется для каждого столбца Y и у1 имеет длину length (xi) - by- size (Y. 2);

ñyi = interpl (x.Y.xi.method) — позволяет с помощью параметра method задать метод интерполяции:

ñ'nearest' — ступенчатая интерполяция;

ñ'linear' — линейная интерполяция (принята по умолчанию);

ñ'spline' — кубическая сплайн-интерполяция;

ñ'cubic' или 'pchip' — интерполяция многочленами Эрмита;

ñ'v5cubic' — кубическая интерполяция MATLAB 5.

ñ yi = interpl (x.Y.xi.method, значение величин вне пределов изменения х) —позволяет отобразить особенные точки на графике;

ñyi = i nterpl(х, Y, xi.method.' сообщение') — позволяет изменить сообщение об особенных точках на графике.

Все методы интерполяции требуют, чтобы значения х изменялись монотонно. Когда х — вектор равномерно распределенных точек, для более быстрой интерполяции лучше использовать методы '*1inear', '*cubic', '*nearest' или '*spline'. Обратите внимание, что в данном случае наименованию метода предшествует знак звездочки.

Пример (интерполяция функции косинуса):

» x=0:10:y=cos(x);

» xi=0:0.1:10:

» yi=interpl(x.y.xi);

» plot(x,y.'x',xi,yi,'g'),hold on

» yi=interpl(x,y.xi.'spline'):

» plot(x.y,'o',xi,yi,'m').grid,hold off

Двумерная табличная интерполяция

Двумерная интерполяция существенно сложнее, чем одномерная, рассмотренная выше, хотя смысл ее тот же — найти промежуточные точки некоторой зависимости z(x, у) вблизи расположенных в пространстве узловых точек. Для двумерной табличной интерполяции используется функция interp2:

ñZI = interp2(X,Y.Z,XI.YI) — возвращает матрицу ZI, содержащую значения функции в точках, заданных аргументами XI и YI, полученные путем интерполяции двумерной зависимости, заданной матрицами X, Y и Z. При этом X и Y должны быть монотонными и иметь тот же формат, как если бы они были получены с помощью функции meshgrid (строки матрицы X являются идентичными; то же можно сказать о столбцах массива Y). Матрицы X и Y определяют точки, в которых задано значение Z. Параметры XI и YI могут быть матрицами, в этом случае interp2 возвращает значения Z, соответствующие точкам (XI(i,j),YI(i.j)). В качестве альтернативы можно передать в качестве параметров вектор-строку xi и вектор-столбец yi. В этом случае interp2 представляет эти векторы так, как если бы использовалась команда mesh-grid(xi.yi);

ñZI = interp2(Z,XJ.YI) — подразумевает, что Х=1:n и Y=l:m, где [m.n]=size(Z);

ñZI = interp2(Z,ntimes) — осуществляет интерполяцию рекурсивным методом с числом шагов ntimes;

ñZI = interp2(X,Y,Z.XI,YI.method) — позволяет с помощью опции method задать метод интерполяции:

ñ'nearest' — интерполяция по соседним точкам;

ñ'linear' — линейная интерполяция;

ñ'cubic' — кубическая интерполяция (полиномами Эрмита);

ñ'spline' — интерполяция сплайнами.

Все методы интерполяции требуют, чтобы X и Y изменялись монотонно и имели такой же формат, как если бы они были получены с помощью функции meshgrid. Когда X и Y — векторы равномерно распределенных точек, для более быстрой интерполяции лучше использовать методы '*1inear', '*cubic', или '*nearest'.

Пример:

» [X.Y]=meshgrid(-3:0.25:3);

Z=peaks(X/2.Y*2):

» [Xl,Yl]=meshgrid(-3:0.1:3);

Zl=interp2(X,Y.Z.Xl.Yl):

» mesh(X.Y,Z).hold on.mesh(Xl.Yl,Zl+15).hold off

Рис. 6.13. Применение функции interpZ

Рис. 6.13 иллюстрирует применение функции interp2 для двумерной интерполяции (на примере функции peaks).

В данном случае поверхность снизу — двумерная линейная интерполяция, которая реализуется по умолчанию, когда не указан параметр method.

Трехмерная табличная интерполяция

Для трехмерной табличной интерполяции используется функция interp3:

ñVI = interp3(X.Y.Z.V.XI,YI.ZI) — интерполирует, чтобы найти VI, значение основной трехмерной функции V в точках матриц XI, YI и ZI. Матрицы X, Y и Z определяют точки, в которых задано значение V. XI, YI и ZI могут быть матрицами, в этом случае InterpS возвращает значения Z, соответствующие точкам (XI (i,j),YI(i. j), ZI (i. j)). В качестве альтернативы можно передать векторы xi, yl и zi. Векторы аргументы, имеющие неодинаковый размер, представляются, как если бы использовалась команда meshgrid;

ñVI = interp3(V.XI.YI.ZI) - подразумевает X=1:N, Y=1:M, Z=1:P, где [M,N.P]=size(V);

ñVI = interpS(V.ntimes) — осуществляет интерполяцию рекурсивным методом с числом шагов ntimes;

ñVI = interp3(....method) — позволяет задать метод интерполяции:

ñ'nearest' — ступенчатая интерполяция;

ñ'linear' — линейная интерполяция;

ñ'cubic' — кубическая интерполяция (полиномами Эрмита);

ñ'spline' — интерполяция сплайнами.

Все методы интерполяции требуют, чтобы X, Y и Z изменялись монотонно и имели такой же формат, как если бы они были получены с помощью функции meshgrid. Когда X и Y и Z — векторы равномерно распределенных в пространстве узловых точек, для более быстрой интерполяции лучше использовать методы '*li'near', '*cubic' или '*nearest'.

N-мерная табличная интерполяция

MATLAB позволяет выполнить даже n-мерную табличную интерполяцию. Для этого используется функция interpn:

ñVI = interpn(Xl.X2,X3,...,V,Y1.Y2.Y3....)- интерполирует, чтобы найти VI, значение основной многомерной функции V в точках массивов Yl, Y2, Y3,.... Функции interpn должно передаваться 2ЛГ+1 аргументов, где N — размерность интерполируемой функции. Массивы XI, Х2, ХЗ,... определяют точки, в которых задано значение V. Параметры Yl, Y2, Y3,... могут быть матрицами, в этом случае interpn возвращает значения VI, соответствующие точкам (YKi, j),Y2(i, j), Y3(i, j),...). В качестве альтернативы можно передать векторы yl, y2, уЗ,... В этом случае interpn интерпретирует их, как если бы использовалась команда ndgrid(yl. У2.у3....);

ñVI = interpn(V.Yl,Y2,Y3,...) - подразумевает Xl=1size(V.l), X2=l:size(V,2), X3=l:size(V,3) и т. д.;

ñVI = Interpn(V.ntimes) — осуществляет интерполяцию рекурсивным методом с числом шагов ntimes;

ñVI = interpn(...method) — позволяет указать метод интерполяции:

ñ'nearest' — ступенчатая интерполяция;

ñ'linear' — линейная интерполяция;

ñ'cubic' — кубическая интерполяция.

В связи с редким применением такого вида интерполяции, наглядная трактовка которой отсутствует, примеры ее использования не приводятся.

Интерполяция кубическим сплайном

Сплайн-интерполяция используется для представления данных отрезками полиномов невысокой степени — чаще всего третьей. При этом кубическая интерполяция обеспечивает непрерывность первой и второй производных результата интерполяции в узловых точках. Из этого вытекают следующие свойства кубической сплайн-интерполяции:

ñграфик кусочно-полиномиальной аппроксимирующей функции проходит точно через узловые точки;

ñв узловых точках нет разрывов и резких перегибов функции;

ñблагодаря низкой степени полиномов погрешность между узловыми точками обычно достаточно мала;

ñсвязь между числом узловых точек и степенью полинома отсутствует;

ñпоскольку используется множество полиномов, появляется возможность аппроксимации функций с множеством пиков и впадин.

Как отмечалось, в переводе spline означает «гибкая линейка». График интерполирующей функции при этом виде интерполяции можно уподобить кривой, по которой изгибается гибкая линейка, закрепленная в узловых точках. Реализуется сплайн-интерполяция следующей функцией:

ñyi = spline(x,y,xi) — использует векторы х и у, содержащие аргументы функции и ее значения, и вектор xi, задающий новые точки; для нахождения элементов вектора yi используется кубическая сплайн-интерполяция;

ñрр = spline(x.y) — возвращает рр-форму сплайна, используемую в функции ppval и других сплайн-функциях.

Пример:

» х=0:10: y=3*cos(x);

» xl=0:0.1:11;

» y1=spline(x,y.xl);

» plot(x.y,'о'.xl.yl.'--')

Сплайн-интерполяция дает неплохие результаты для функций, не имеющих разрывов и резких перегибов. Особенно хорошие результаты получаются для монотонных функций.

Результат интерполяции показан на рис. 6.14.

Рис. 6.14. Пример применения функции spline

Ввиду важности сплайн-интерполяции и аппроксимации в обработке и представлении сложных данных в состав системы MATLAB входит пакет расширения Spline Toolbox, содержащий около 70 дополнительных функций, относящихся к реализации сплайн-интерполяции и аппроксимации, а также графического представления сплайнами их результатов. Для вызова данных об этом пакете (если он установлен) используйте команду help splines.

⇐ Предыдущая 1 2 345 Следующая ⇒