Может ли кто-нибудь дать мне некоторое направление для репликации функции MATLAB interp1, используя сплайн-интерполяцию? Я попытался тщательно воспроизвести алгоритм на странице wikipedia, но результаты на самом деле не совпадают.
#include <stdio.h>
#include <stdint.h>
#include <iostream>
#include <vector>
//MATLAB: interp1(x,test_array,query_points,'spline')
int main(){
int size = 10;
std::vector<float> test_array(10);
test_array[0] = test_array[4] = test_array[8] = 1;
test_array[1] = test_array[3] = test_array[5] = test_array[7] = test_array[9] = 4;
test_array[2] = test_array[6] = 7;
std::vector<float> query_points;
for (int i = 0; i < 10; i++)
query_points.push_back(i +.05);
int n = (size - 1);
std::vector<float> a(n+1);
std::vector<float> x(n+1); //sample_points vector
for (int i = 0; i < (n+1); i++){
x[i] = i + 1.0;
a[i] = test_array[i];
}
std::vector<float> b(n);
std::vector<float> d(n);
std::vector<float> h(n);
for (int i = 0; i < (n); ++i)
h[i] = x[i+1] - x[i];
std::vector<float> alpha(n);
for (int i = 1; i < n; ++i)
alpha[i] = ((3 / h[i]) * (a[i+1] - a[i])) - ((3 / h[i-1]) * (a[i] - a[i-1]));
std::vector<float> c(n+1);
std::vector<float> l(n+1);
std::vector<float> u(n+1);
std::vector<float> z(n+1);
l[0] = 1.0;
u[0] = z[0] = 0.0;
for (int i = 1; i < n; ++i){
l[i] = (2 * (x[i+1] - x[i-1])) - (h[i-1] * u[i-1]);
u[i] = h[i] / l[i];
z[i] = (alpha[i] - (h[i-1] * z[i-1])) / l[i];
}
l[n] = 1.0;
z[n] = c[n] = 0.0;
for (int j = (n - 1); j >= 0; j--){
c[j] = z[j] - (u[j] * c[j+1]);
b[j] = ((a[j+1] - a[j]) / h[j]) - ((h[j] / 3) * (c[j+1] + (2 * c[j])));
d[j] = (c[j+1] - c[j]) / (3 * h[j]);
}
std::vector<float> output_array(10);
for (int i = 0; i < n-1; i++){
float eval_point = (query_points[i] - x[i]);
output_array[i] = a[i] + (eval_point * b[i]) + ( eval_point * eval_point * c[i]) + (eval_point * eval_point * eval_point * d[i]);
std::cout << output_array[i] << std::endl;
}
system("pause");
return 0;
}
Оглядываясь назад, ваш код, кажется, правильно закодирован со ссылкой на статью в Википедии. Тем не менее, есть что-то, что вам нужно знать о interp1
который я не думаю, что вы учли при использовании его, чтобы проверить свои ответы.
MATLAB interp1
когда вы указываете флаг spline
предполагает, что условия конечной точки не являются узлами. Алгоритм, указанный в Википедии, является кодом для естественного сплайна.
Таким образом, вероятно, поэтому ваши очки не совпадают. FWIW, проконсультируйтесь: http://www.cs.tau.ac.il/~turkel/notes/numeng/spline_note.pdf и посмотрите на диаграмму на последней странице. Вы увидите, что сплайны с не-узлом и естественные сплайны имеют одинаковую форму, но имеют разные значения y, когда ваши данные состоят только из конечных точек вашего сплайна. Однако, если у вас есть точки данных между конечными точками, все разные типы сплайнов (более или менее) имеют одинаковые значения y.
Для полноты, вот цифра, извлеченная из заметок PDF, на которые я ссылался выше:
Если вы хотите использовать естественные сплайны, используйте csape
вместо interp1
. Это обеспечивает кубический сплайн с конечными условиями. Вы вызываете csape
следующим образом:
pp = csape(x,y);
x
и y
- контрольные точки, определенные для вашего сплайна. По умолчанию это возвращает естественный сплайн, который является тем, что вам нужно, и является структурой типа ppform
. Затем вы можете выяснить, что сплайн оценивает с помощью fnval
:
yval = fnval(pp, xval);
xval
и yval
- входная координата x
а результат вычисляется для сплайна при этом конкретном x
.
Используйте это, а затем проверьте, соответствует ли ваш код значениям, предоставленным csape
.
Для использования csape
вам нужна панель инструментов Curve Fitting Toolbox в MATLAB. Если у вас этого нет, то, к сожалению, этот метод не будет работать.
Я думаю, что interp1
поддерживается MATLAB CODER.
Просто используйте CODER, чтобы сгенерировать код C, и у вас есть то, что вам нужно.
main()
зачем вам нужны указатели вместо локальных массивов или векторов?