Суть вопроса: есть массив данных (40 отсчетов), требуется выполнить интерполяцию таким образом, чтобы повысит частоту дискретизации в 5 раз, т.е. между имеющимися отсчетами поместить еще по 4 отсчета. Внимание привлек вот этот рекуррентный алгоритм. Здесь не требуется решать систему уравнений, вычислений немного, можно вычислять потоково. Для параболической сплайн интерполяции алгоритм заработал (проверял на MATLAB). А вот интерполяция кубическими полиномами работать не хочет. Верней работает неправильно - функции в конечной точке отрезка интерполяции терпит разрыв. Не могу найти причину.
Вот график
Нажмите для просмотра прикрепленного файла
Может кто-нибудь подсказать, что не так с алгоритмом, и что можно сделать?
На всякий случай прикрепляю код скрипта (может в нем накосячил?). Формулы в исходнике немного отличаются, но с ними вообще не работал. Поэтому проверил математику и переписал формулы.
Код
% интерполяция кубическим сплайном
x = [1 0 1 1 1 1 1 1 2 3 7 10 3 1 1 1 4 5 5 2 1 0 0 0 0 1 3 1 1 1];
x1 = zeros(1,5*length(x));
for i = 0:length(x)-1
x1(5*i+1) = x(i+1); % известные значения оставляем на своих местах
end
y = x1;
df = -0.1;
d2f = 0;
for i = 1:length(x)-1
if i == 1
a = 0.04*(0.2*(x(2)-x(1))-df-d2f*5)*(1/6);
b = 0;
c = df;
d = x(1);
for j = 1:4
x1(j+1) = a*j^3 + b*j^2 + c*j + d;
end
else
x_pos = (5*i-5); % позиция текущего интервала анализа в новом массиве x1
beta = 3*a*x_pos^2 + 2*b*x_pos + c;
gamma = 6*a*x_pos + 2*b;
a = 0.04*(0.2*(x(i+1)-x(i))-5*gamma-beta)*(1/6);
b = 0.5*gamma - 3*a*x_pos;
c = beta - 3*a*x_pos^2 - 2*b*x_pos;
d = x(i)-a*x_pos^3-b*x_pos^2-c*x_pos;
for j = x_pos+1:x_pos+4 % известные точки остаются, рассчитанные помещаем между ними
x1(j+1) = a*j^3 + b*j^2 + c*j + d;
end
end
end
plot(0:length(x1)-1,y,0:length(x1)-1,x1)
x = [1 0 1 1 1 1 1 1 2 3 7 10 3 1 1 1 4 5 5 2 1 0 0 0 0 1 3 1 1 1];
x1 = zeros(1,5*length(x));
for i = 0:length(x)-1
x1(5*i+1) = x(i+1); % известные значения оставляем на своих местах
end
y = x1;
df = -0.1;
d2f = 0;
for i = 1:length(x)-1
if i == 1
a = 0.04*(0.2*(x(2)-x(1))-df-d2f*5)*(1/6);
b = 0;
c = df;
d = x(1);
for j = 1:4
x1(j+1) = a*j^3 + b*j^2 + c*j + d;
end
else
x_pos = (5*i-5); % позиция текущего интервала анализа в новом массиве x1
beta = 3*a*x_pos^2 + 2*b*x_pos + c;
gamma = 6*a*x_pos + 2*b;
a = 0.04*(0.2*(x(i+1)-x(i))-5*gamma-beta)*(1/6);
b = 0.5*gamma - 3*a*x_pos;
c = beta - 3*a*x_pos^2 - 2*b*x_pos;
d = x(i)-a*x_pos^3-b*x_pos^2-c*x_pos;
for j = x_pos+1:x_pos+4 % известные точки остаются, рассчитанные помещаем между ними
x1(j+1) = a*j^3 + b*j^2 + c*j + d;
end
end
end
plot(0:length(x1)-1,y,0:length(x1)-1,x1)