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

Octave поддерживает несколько методов одномерной интерполяции, большинство из которых описаны в этом разделе. Полиномиальная интерполяция и Интерполяция на рассеянных данных описывают дополнительные методы.

: yi = interp1 (x, y, xi) ¶

: yi = interp1 (y, xi) ¶

: yi = interp1 (…, method) ¶

: yi = interp1 (…, extrap) ¶

: yi = interp1 (…, "left") ¶

: yi = interp1 (…, "right") ¶

: pp = interp1 (…, "pp") ¶

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

Интерполирует входные данные, чтобы определить значение yi в точках xi. Если не указано, x принимается как индексы y (1:length (y)). Если y является матрицей или N-мерным массивом, интерполяция выполняется для каждой колонки y.

Метод интерполяции method может быть одним из:

"nearest"

Возвращает ближайшего соседа.

"previous"

Возвращает предыдущего соседа.

"next"

Возвращает следующего соседа.

"linear" (по умолчанию)

Линейная интерполяция от ближайших соседей.

"pchip"

Кубическая интерполирующая полиномиальная функция Эрмита — сохраняющая форму интерполяция с плавным первым производным.

"cubic"

Кубическая интерполяция (то же самое, что и "pchip").

"spline"

Кубическая сплайн-интерполяция — гладкие первые и вторые производные по всей кривой.

Добавление '*' в начало любого метода выше заставляет interp1 предположить, что x равномерно распределён, и используются только x(1) и x(2). Это обычно быстрее и никогда не медленнее. По умолчанию используется "linear".

Если extrap — строка "extrap", то значения за пределами конечных точек экстраполируются с использованием текущего method. Если extrap — число, то значения за пределами конечных точек заменяются этим числом. По умолчанию extrap — NA.

Если указан строковый аргумент "pp", то xi не должно быть указано, и interp1 возвращает объект кусочно-полиномиальной функции. Этот объект может быть позднее использован с ppval для вычисления интерполяции. Существует эквивалентность, такая что ppval (interp1 (x, y, method, "pp"), xi) == interp1 (x, y, xi, method, "extrap").

Дублирование точек в x задает разрывную интерполянт. Может быть не более 2 последовательных точек с одинаковым значением. Если x возрастает, по умолчанию разрывная интерполянта является правонепрерывной. Если x убывает, по умолчанию разрывная интерполянта является левонепрерывной. Условие непрерывности интерполянта может быть задано с помощью опций "left" или "right" для выбора левонепрерывной или правонепрерывной интерполянты соответственно. Разрывная интерполяция допускается только для методов "nearest" и "linear"; во всех остальных случаях значения x должны быть уникальными.

Пример использования interp1:

xf = [0:0.05:10];
yf = sin (2*pi*xf/5);
xp = [0:10];
yp = sin (2*pi*xp/5);
lin = interp1 (xp, yp, xf);
near = interp1 (xp, yp, xf, "nearest");
pch = interp1 (xp, yp, xf, "pchip");
spl = interp1 (xp, yp, xf, "spline");
plot (xf,yf,"r", xf,near,"g", xf,lin,"b", xf,pch,"c", xf,spl,"m",
      xp,yp,"r*");
legend ("original", "nearest", "linear", "pchip", "spline");

См. также: pchip, spline, interpft, interp2, interp3, interpn.

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

Это можно продемонстрировать кодом:

t = -2:2;
dt = 1;
ti =-2:0.025:2;
dti = 0.025;
y = sign (t);
ys = interp1 (t,y,ti,"spline");
yp = interp1 (t,y,ti,"pchip");
ddys = diff (diff (ys)./dti) ./ dti;
ddyp = diff (diff (yp)./dti) ./ dti;
figure (1);
plot (ti,ys,"r-", ti,yp,"g-");
legend ("spline", "pchip", 4);
figure (2);
plot (ti,ddys,"r+", ti,ddyp,"g*");
legend ("spline", "pchip");

Результат которого можно увидеть на рисунках 29.1 и 29.2.

Интерполяция Фурье — метод ресемплирования, в котором сигнал преобразуется в частотную область, дополняется нулями, а затем преобразуется обратно во временную область.

: interpft (x, n) ¶

: interpft (x, n, dim) ¶

Интерполяция Фурье.

Если x — вектор, то x ресемплируется до n точек. Предполагается, что данные в x равноотстоящие. Если x — матрица или N-мерный массив, интерполяция выполняется для каждой колонки x.

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

interpft предполагает, что интерполируемая функция периодична, поэтому делаются предположения относительно крайних точек интерполяции.

См. также: interp1.

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

t = 0 : 0.3 : pi; dt = t(2)-t(1);
n = length (t); k = 100;
ti = t(1) + [0 : k-1]*dt*n/k;
y = sin (4*t + 0.3) .* cos (3*t - 0.1);
yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1);
plot (ti, yp, "g", ti, interp1 (t, y, ti, "spline"), "b", ...
      ti, interpft (y, k), "c", t, y, "r+");
legend ("sin(4t+0.3)cos(3t-0.1)", "spline", "interpft", "data");

который демонстрирует плохое поведение интерполяции Фурье для непериодических функций, как можно увидеть на рисунке 29.3.

Кроме того, вспомогательные функции spline и lookup, лежащие в основе функции interp1, могут вызываться напрямую.

: pp = spline (x, y) ¶

: yi = spline (x, y, xi) ¶

Возвращает кубический сплайн-интерполянт для точек x и y.

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

При вызове с третьим входным аргументом, spline вычисляет сплайн в точках xi. Третья форма вызова spline (x, y, xi) эквивалентна ppval (spline (x, y), xi).

Переменная x должна быть вектором длины n.

y может быть вектором или массивом. Если y — вектор, его длина должна быть либо n, либо n + 2. Если длина y равна n, то используется граничное условие "not-a-knot". Если длина y равна n + 2, то первые и последние значения вектора y — значения первой производной кубического сплайна в конечных точках.

Если y — массив, то размер y должен иметь вид [s1, s2, …, sk, n] или [s1, s2, …, sk, n + 2]. Массив преобразуется внутри в матрицу, где ведущее измерение задаётся s1 * s2 * … * sk, и затем каждая строка этой матрицы обрабатывается отдельно. Обратите внимание, что это прямо противоположно interp1, но сделано для совместимости с MATLAB.

См. также: pchip, ppval, mkpp, unmkpp.