Победить Фарроу - 2

[_Ivana]

Небольшая предыстория. Есть один популярный вид локальной интерполяции, когда по 4 точкам строится полином Лагранжа 3 степени, проходящий через эти 4 точки, и значения этого полинома берутся в качестве интерполирующей функции на среднем интервале, задаваемом второй и третьей точкой. Для соседнего интервала строится свой такой же полином по 4 точкам, и т.д. Этот вид интерполяции широко применяется во многих практических приложениях, особенно в ЦОС (сразу оговоримся, что далее везде подразумевается равномерная сетка), он настолько популярен, что ему посвящается много статей, например вот, где приведена даже модификация Фарроу для оптимизации расчета коэффициентов полинома: требуется 1 деление на 6, 2 умножения/деления на 2 (сдвига при реализации алгоритма в машинном коде) и 8 сложений/вычитаний. Некоторое время назад я в определенном смысле "победил" этот алгоритм - нашел другой метод интерполяции, который требует для расчета коэффициентов столько же операций, но при этом на порядок (в 10 раз) точнее и как бонус имеет непрерывную первую производную. Даже написал маленькую статейку в инете по этому поводу, но это никого не заинтересовало

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

ЗЫ те 2 участника этого форума, которые прочли ту мою статейку: если вы там обнаружили подсказки к решению этой задачи, не спешите пока их обнародовать - я лелею надежду, что другие участники могут заинтересоваться и захотеть подумать самостоятельно

[_Ivana]

В теге оффтопа код демо примера на Матлабе, если кто хочет подсмотреть ответ. В коде только конечные формулы расчета и построение графика, без теоретических выкладок и доказательств алгоритма.
ЗЫ при выборе синтаксиса Матлаба или С текст кода отображается почему-то некорректно, поэтому выбрал свободный синтаксис.

(Оффтоп)

1. clf reset 
2. N = 10; h = 0.01; x = 1:N; y = rand(1, N); 
3. plot(x, y, 'or', 'LineWidth', 2); hold on; grid on; axis on; 
4. title(['Локальная интерполяция случайного набора точек', ... 
5. ' полиномом Лагранжа 3 степени, 2 алгоритма расчета.']); 
6. plot(x(1), y(1), 'b-', x(1), y(1), 'g:'); 
7. legend('точки', 'мой алгоритм', 'Фарроу'); 
8. for k = 1:(N-2) 
9.     % мой алгоритм: 1 умножение, 1 сдвиг, 6 сложений 
10.     c = d; e = y(k+1) - y(k); d = y(k+2) - y(k+1) - e; 
11.     b2 = c/2; b3 = (d - c)/6; b1 = e - b2 - b3; 
12.      
13.     if (k == 1) continue; end 
14.     t = 0:h:1; f = b3.*t.^3 + b2.*t.^2 + b1.*t + y(k); 
15.     plot(t+x(k), f, 'b-') 
16.      
17.     % Фарроу: 1 умножение, 2 сдвига, 8 сложений 
18.     a3 = (y(k+2) - y(k-1))/6 + (y(k) - y(k+1))/2; 
19.     a1 = (y(k+2) - y(k))/2 - a3; 
20.     a2 =  y(k+2) - y(k+1) - a1 - a3; 
21.      
22.     t = -1:h:0; f = a3.*t.^3 + a2.*t.^2 + a1.*t + y(k+1); 
23.     plot(t+1+x(k), f, 'g:') 
24. end 

[profrotter]

_Ivana в сообщении #785475 писал(а):

ему посвящается много статей, например вот,

Статью в топку, автора - на костёр! За ортогональные многочлены Лагранжа.

[_Ivana]

На другом форуме один участник выложил блок-схему своего варианта оптимизации, я перевел ее в свой скрипт. Его вариант оптимизации основан на алгоритме, представленном в ссылке первого поста, и уменьшает последний на 2 сложения, уступая теперь моему алгоритму только один сдвиг. Мой же алгоритм, в свою очередь, построен на другой математической модели, нежели стандартная постановка задачи построения полинома по 4 точкам, что и позволяет оптимизировать стандартный вариант.
Скрипт в теге оффтопа:

(Оффтоп)

1. clf reset 
2. N = 10; x = 1:N; y = rand(1, N); h = 0.01; d = 0; 
3. plot(x, y, 'or', 'LineWidth', 2); hold on; grid on; axis on; 
4. title(['Локальная интерполяция случайного набора точек', ... 
5. ' полиномом Лагранжа 3 степени, 2 алгоритма расчета.']); 
6. plot(x(1), y(1), 'b-', x(1), y(1), 'g:'); 
7. legend('точки', 'мой алгоритм', 'Фарроу по Anatoliy'); 
8. for k = 1:(N-2) 
9.     % мой алгоритм: 1 умножение, 1 сдвиг, 6 сложений 
10.     c = d; e = y(k+1) - y(k); d = y(k+2) - y(k+1) - e; 
11.     b2 = c/2; b3 = (d - c)/6; b1 = e - b2 - b3; 
12.      
13.     if (k == 1) continue; end 
14.     t = 0:h:1; f = b3.*t.^3 + b2.*t.^2 + b1.*t + y(k); 
15.     plot(t+x(k), f, 'b-') 
16.      
17.     % Фарроу по _Anatoliy: 1 умножение, 2 сдвига, 6 сложений        
18.     p = y(k) - y(k+1); q = (y(k+2) - y(k))/2; 
19.     a3 = ( (y(k+2) - y(k-1))/3 + p )/2; 
20.     a2 = p + q; a1 = q - a3; 
21.      
22.     t = -1:h:0; f = a3.*t.^3 + a2.*t.^2 + a1.*t + y(k+1); 
23.     plot(t+1+x(k), f, 'g:') 
24. end 

[_Ivana]

На досуге подумал, убрал еще одно сложение. Теперь у меня 1 умножение, 1 сдвиг и 5 сложений против лучшей представленной альтернативы 1 умножение, 2 сдвига и 6 сложений:

Используется синтаксис Matlab M

clf reset
N = 10; x = 1:N; y = rand(1, N);
h = 0.01; g = y(3) - y(2); d = g - y(2) + y(1);
plot(x, y, 'or', 'LineWidth', 2); hold on; grid on; axis on;
title(['Локальная интерполяция случайного набора точек', ...
' полиномом Лагранжа 3 степени, 2 алгоритма расчета.']);
plot(x(1), y(1), 'b-', x(1), y(1), 'g:');
legend('точки', 'мой алгоритм', 'Фарроу по Anatoliy');
for k = 2:(N-2)
    % мой алгоритм: 1 умножение, 1 сдвиг, 5 сложений
    c = d; e = g; g = y(k+2) - y(k+1); d = g - e;
    b2 = c/2; b3 = (d - c)/6; b1 = e - b2 - b3;
    t = 0:h:1; f = b3.*t.^3 + b2.*t.^2 + b1.*t + y(k);
    plot(t+x(k), f, 'b-')
    
    % Фарроу по _Anatoliy: 1 умножение, 2 сдвига, 6 сложений       
    p = y(k) - y(k+1); q = (y(k+2) - y(k))/2;
    a3 = ( (y(k+2) - y(k-1))/3 + p )/2;
    a2 = p + q; a1 = q - a3;
    
    t = -1:h:0; f = a3.*t.^3 + a2.*t.^2 + a1.*t + y(k+1);
    plot(t+1+x(k), f, 'g:')
end

[ewert]

_Ivana в сообщении #785475 писал(а):

локального сплайна, составленного именно из полиномов Лагранжа,

Так не бывает. Или сплайн, или Лагранж.

[_Ivana]

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

[ewert]

_Ivana в сообщении #864465 писал(а):

Я считаю, что функция удовлетворяет всем условиям, чтобы носить гордое имя сплайна.

Носить она может что угодно, однако чтобы носить с достоинством -- должна быть хотя бы один раз гладкой.

[_Ivana]

На это у меня есть два пункта:
1) Ломаная по заданным точкам ни разу не гладкая, однако носит это имя, не знаю насчет достоинства.
2) Самое забавное, что рассматриваемая в теме интерполяция может быть названа квазигладкой в определенном смысле (я старался сказать как можно более аккуратно, чтобы не вызвать праведного негодования ). И именно этот факт позволил мне оптимизировать процедуру ее численного расчета по сравнению с подходом "в лоб", лучший результат которого продемонстрирован в альтернативном расчете. И, кстати, этот факт я вскользь отмечал с той моей так и нигде не опубликованной статейке, которую вы некоторое время назад просмотрели и покритиковали.

[profrotter]

И там критиковали и тут укажем, что для любой локальной регулярной линейной интерполяции, при которой на интервале $t\in[t_k,t_{k+1}]$ интерполирующая функция $\psi(t)$ зависит от $M$ предшествующих и $M$ последующих отсчётов, фрагмент инерполирующей функции $\psi_k(t)$, соответствующий этому интервалу, описывается выражением:

$$\psi_k(t)=\sum\limits_{m=-M}^{M+1}s_{k+m}\varphi_0(t-t_{k+m}),$$

где $s_k$ - дискретные значения сигнала, $\varphi_0(t)$ - порождающая функция метода интерполяции, $t_k=kT$ - момент дискретизации, $T$ - период дискретизации, $k=0,\pm 1,\pm2,...$

Конкретно в рассматриваемом вами случае мы имеем дело с локальной полиномиальной интерполяцией порядка $M=1$, при этом:

$$\psi_k(t)=s_{k-1}\varphi_0(t-t_{k-1})+s_k\varphi_0(t-t_k)+s_{k+1}\varphi_0(t-t_{k+1})+s_{k+2}\varphi_0(t-t_{k+1}).$$

Выражение для порождающей функции приведено, например, в статье http://www.radiotec.ru/catalog.php?cat=jr8&art=12716

$$\varphi_0(t)=\begin{cases}\frac 1 2 \left|\frac t T\right|^3-\left(\frac t T\right)^2-\frac 1 2 \left|\frac t T\right|+1,|t|\leqslant T\\-\frac 1 6 \left|\frac t T\right|^3+\left(\frac t T\right)^2-\frac {11} {6} \left|\frac t T\right|+1,T\leqslant |t|\leqslant 2T\end{cases}$$

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

Сам алгоритм получения нового отсчёта при локальной регулярной интерполяции первого порядка, таким образом, включает четыре умножения и три сложения, независимо от метода интерполяции.

[_Ivana]

profrotter, если я правильно понимаю, вы записали выражение для интерполирующей функции в виде линейной комбинации некоторой "порождающей функции" с весами, равными значениям сигнала. Тогда конкретно в рассматриваем мною случае, эта "порождающая функция" будет представлять собой базисные полиномы Лагранжа. И с этим никто не спорит, более того, это было отмечено и в статье по ссылке из первого поста этого топика, за которую вы хотели отдать автора на костер - за "ортогональность" этих базисных полиномов Лагранжа.
Но когда доходит до практических расчетов, тогда именно что рассчитывают по нескольким точкам коэффициенты полинома и потом его значения в любом количестве внутренних точек интервала - так вычислительно эффективнее. А ваша фраза

profrotter в сообщении #864521 писал(а):

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

Сам алгоритм получения нового отсчёта при локальной регулярной интерполяции первого порядка, таким образом, включает четыре умножения и три сложения, независимо от метода интерполяции.

- это представление полиномиальной интерполяции в виде полифазного КИХ фильтра с заранее рассчитанными коэффициентами для каждой фазы - и с этим тоже никто не спорит. Вопрос только в том, что не всегда количество фаз заранее известно и конечно, а даже если и конечно - то может быть настолько велико, что хранение всех заранее рассчитанных коэффициентов потребует много памяти.

ЗЫ вот по этой ссылке тема, где обсуждается вопрос "что лучше - Лагранж или полифазный фильтр"
ЗЗЫ и в любом случае, независимо от личных предпочтений и взглядов на то, кто что использует и как вычисляет, это не делает полученный мной результат неверным.

[profrotter]

_Ivana в сообщении #864532 писал(а):

Тогда конкретно в рассматриваем мною случае, эта "порождающая функция" будет представлять собой базисные полиномы Лагранжа.

В конкретно рассматриваемом вами случае выражение для порождающей функции было приведено. Это не базисная функция Лагранжа. Базисные функции Лагранжа не ортогональны.

_Ivana в сообщении #864532 писал(а):

так вычислительно эффективнее.

Нисколько. Рассчитав коэффициенты полинома делается три умножения и три сложения: $\psi_k(t)=a_{k3}t^3+a_{k2}t^2+a_{k1}t+a_{k0}$. Прибавьте к этому умножения, которые необходимы для расчёта коэффициентов и получите четыре (а то и больше) умножения и три сложения, о которых я вам писал.

_Ivana в сообщении #864532 писал(а):

Вопрос только в том, что не всегда количество фаз заранее известно и конечно,

Не надо инсинуаций. Если Вы хотите оперировать этим, то не ссылайтесь на цифровую обработку сигналов. А если нет - то отпадает и необходимость во всей вашей оптимизации.

_Ivana в сообщении #864532 писал(а):

что лучше - Лагранж или полифазный фильтр

Бессмысленно. Полифазный фильтр может быть построен для любого метода линейной локальной регулярной интерполяции, в том числе и для локальной полиномиальной, также для локальной сплайновой. Как это сделать я написал выше.

_Ivana в сообщении #864532 писал(а):

полученный мной результат

Результата Вы не предъявили. Копаться в чужом программном коде просто никто не стал и всё. Опишите алгоритм - будет о чём разговаривать.

[_Ivana]

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

[profrotter]

_Ivana в сообщении #864543 писал(а):

Я вам привел прямую ссылку на тему со специализированного форума, в которой люди, не менее вас

Ссылки на авторитеты не катят, тем более, что эти авторитеты считают ортогональными базисные функции Лагранжа. Я вам ещё раз пишу, что никакой алгоритм Вы не описали. Без этого тема ни о чём - капаться в чужом программном коде, реализующем непонятно что бессмысленно.

[_Ivana]

Хотелось бы подвести черту под бесплодными спорами и закончить их. А также резюмировать и напомнить: тема в разделе "Олимпиадные задачи." В моем коде вообще копаться не надо, надо задачу решить
А она такова: есть некая процедура интерполяции, описанная по ссылке. Если надо, могу продублировать сюда ее описание. В той же ссылке приведен так называемый "оптимизированный алгоритм", реализующий эту процедуру, и содержащий $8$ операций сложения, $2$ сдвига и $1$ умножение для получения коэффициентов полинома каждого интервала. Собственно, задача - предложить и реализовать более оптимальный алгоритм.
ЗЫ а в коде (в котором не надо копаться, там готовые ответы) - результат участника другого форума, который оптимальнее "оптимизированного" по ссылке на $2$ сложения и мой, который оптимальнее на $3$ сложения и $1$ сдвиг.