Помогите опознать алгоритм пожалуйста.

Для размерности 3 выглядит вот так:



a, b - исходные векторы размерности 3, a0 - результат.

В коде видны явные закономерности, похоже на какую-то стандартную операцию над векторами, но не могу понять какую именно.

Пусть у Вас есть матрица сдвига P

(0 1 0)
(0 0 1)
(1 0 0)

a и b - векторы, тогда матрица H=
( a, b ) ^T P ( a, b )
на диагонали будет содержать a0 и a2, а сумма внедиагональных членов будет равна a1.

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

Вот более полный код, для N=3 эквивалентен вышеприведеному, за исключением посленей части:
Это кусок мат. модели, векторы a и b вычисляются из физических параметров системы.

После этого из массивов а0 и a1 формируются коэфициенты полинома (комплексные):
t0 и t1 - тоже вычисляются из параметров, представляют собой выражения вида t0 = A * sin(W), t1 = A * cos(W).

После этого находятся корни полинома, и они используются дальше в модели.

a1 инициализирован только в первых 4 элементах. А дальше, в общем случае, тоже нули или модель просто не гоняли? Циклы в общем случае по i,j,k всегда до одного и того же N бегут? С виду очень походе на линейное предсказание с преобразованием исходных реальных векторов в комплексные, но и на сто процентов сказать не могу, если есть еще какая-то информация об алгоритме, говорите!

Да, нули. Но в жизни он будет использоваться только с 2 и 3 элементами.

Да, до одного (иногда до N, иногда до N+1).

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

Вот код инициализации a[] и b[]:

Тут arr[] - массив немного отличающихся частот, имитация 3 немного расстоенных струн в хоре одной клавиши, похоже. overtone - номер текущего обертона.

На выходе (после нахождения корней полинома) получается 3 комплексных числа, действительная часть которых деленная на 2*pi равна "настоящему" (тому что используется непостредственно при синтезе) смещению частоты от центральной (частоты ноты), а мнимая - напрямую скорости затухания струны (показателю экспоненты).

Массив arr[] получается путем сравнительно сложных табличных преобразований, но это не важно, можно просто считать, что он содержит числа примерно равные частоте струн данной ноты, с небольшой расстройкой. Расстройка растет с номером гармоники, на выходе получается похожая расстройка (тоже растет с номером гармоники), но большая по абсолютной величине.

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

PPS: В худшем случае я просто выкину все эти пробразования, заменив их на таблицы/аппроксимации, но интересно будет понять что-же тут происходит, может когда пригодится еще.

Немного становится понятнее. В цикле тут происходит перемножение полиномов. Начальный полином - это константа 1. после чего он домножается на (b[j] x + a[j]), где j != k. и результат складывается.

Теперь нужно понять какой в этом физический смысл.

Вот эквивалентный код на Matlab-e:


Интересное наблюдение - если в arr все числа одинаковые, то два корня полинома z равны нулю, а один не равен, именно он определяет частоту и скорость затухания. То-есть вся эта процедура стремится распределить скорости затухания по струнам так, чтобы имитировать двухстадийное затухание (резкий спад, потом плавное затухание), присущее настояшему музыкальному инструменту.

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

Добый день, Александр,

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

С уважением

ИИВ

LP - это линейное программирование?

Если что, вот статья про которую я говорю: https://ccrma.stanford.edu/~stilti/papers/rltalk.pdf

Вот немного более читаемый матлабовский код:

В результате что полностью согласуется с выводами в статье - получилось 3 "струны", одна затухает очень быстро (7.8407), а две остальные медленно (0.0279 и 0.0056) при этом они расстроены по частоте так, что образуются биения.

Я с корневыми годографами на "вы", так что осознать пока все это не могу, но похоже, что t1 - это то, что в статье обозвано коэффициентом связи (coupling coefficient). В параметрах модели он плавно меняется от 0 до -15 градусов (тут он переведен в радианы уже) при изменении номера клавиши от первой до последней.

Я думал, что линейное предсказание, но, похоже, ошибался.

Если представить $a_i$, $b_i$ как входные векторы, то результат в матричном виде и латехе можно записать примерно так:

$$\left( \product_i (a_i I + b_i P) \right) \sum_i (a_i I + b_i P)^{-1} e,$$

где $e = (1, 0, ..., 0)^T$, а P матрица сдвига, как я сказал выше.

Правда мне такое представление ничего не дало чтобы упростить или как-то серьезно перекаверкать этот алгоритм...

Похоже и не получится, похоже, что это результат не матричных операций. Плохо, что не ясен даже физический смысл a[], b[], t0 и t1.

В любом случае я продолжаю читать разные статьи и пробовать разные варианты, возможно получится решить задачу самостоятельно и другим методом.

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

Похоже, что решение состоит из комбинации экспонент (как минимум 4 штуки на 1 график).

PS: на изменение цвета линии внимание можно не обращать - это просто корни меняются местами.