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

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

Но я не могу придумать/найти способ записать это математически в виде пригодном для моделирования.

Есть идеи как это можно сделать?

Или же помогите понять что тут происходит: алгоритм получен в результате реверс-инжениринга, так что есть неясности.

На входе есть массив arr = [261.4508114, 261.5063567, 261.6451613] - это частота основного тона ноты с небольшой расстройкой (в данном случае это нота C4).

и пара чисел t0 = 199.1735851 и t1 = -0.00705058376. Первое высчитывается сложным образом по таблицам, второе - плавно изменятется в зависимости от основной частоты в пределах от 0 до -0.08 радиан.

Алгоритм выглядит так:

avg = sum(arr) / length(arr);

t[i] = 1 - arr[i] / avg;
a[i] = (t[i]+1) / (avg * 2);
b[i] = t[i] * (partial) * pi; // partial - номер гармоники (1 для основного тона).

Далее составляются полиномы (правое число - это свободный член, дальльше идут по увеличению степени): P0 = {a[0], b[0]}, P1 = {a[1], b[1]}, P2 = {a[2], b[2]}, Px = {a[1] b[2]}.

Далее считаются полиномы A0 = P1*P2 + P0*P2 + P0*P1 и A1 = P0*P1*Px;

Из всего этого считается полином Z = t0*A1*exp(j*t1) - j*A0, где t0, t1 - константы (см.выше), A0, A1 - полиномы, j - мнимая единица.

Далее ищутся корни полинома Z, действительная часть которых имеет смысл частоты, а минмая - скорости затухания данной гармоники.

Есть ощущение, что это какой-то стандартный алгоритм.

В общем случае, для любой длинны массива A0 считается как сумма произведений всех полиномов, только из каждого из произведений один исходный полином выкинут, Px всегда состоит из элементов {a[N-2], b[N-1]}, A1 - это произведение последнего члена в A0 на Px.