Поиск аналитической формулы для набора данных

Метод наименьших квадратов.

если сумма экспонент, то надо сделать преобразование Фурье и посмотреть, что же это за экспоненты

Я его полностью не понимаю :) Можно общую схему действий в данной ситуации?

Попробуйте что-нибудь типа A*exp( - (x - x0)^2/s^2) + B*exp(x-x1), только «хвост» больше похож на логарим или корень.

Пусть есть множество точек {(x_i, y_i)}, i=1..n. Пусть есть функция f(x, c_1, ..., c_m) зависяшая от параметров c_1, ..., c_m, которой мы хотим найти (например f(x, c_1, c_2) = c_1*x + c_1).

Постоим функцию s(c_1, ..., c_2) = \sum_{i=1}^{n} ( y_i - f(x_i, c_1, ..., c_m)^2. Чем меньше s, тем ближе f проходит к исходным точкам. Так же очевидно, что s >= 0. Поэтому будем искать её минимум. Для этого нужно посчитать и приравнять к нулю частные производные по каждому из c_1, ..., c_m. Получим систему из m уравнений на m неизвестных c_1, ..., с_m.

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

Во. Нормальное краткое и ясное объяснение. Спасибо.

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

Да, всё так. Можно конечно попробовать в конечных разностях посчитать, но думаю точность не очень получиться. А от экспонент производные считать одно удовольствие :) Саму систему можно решить каким-нибудь численным методом.

Как уже заметили, тут хорош метод наименьших квадратов. Это умеют делать многие программы, например gnuplot.

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

Чем больше у вас параметров, тем хуже будет работать МНК. Да и зависимость у вас экспоненциальная - тут вполне реально, что «решение» у вас будет очень далеким от идеала.

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

Конец главы 2 (вкратце).

Мне нужно аналитическое решение. Это я пытаюсь понять алгоритм из этого топика http://www.linux.org.ru/forum/development/6871204. Эта функция - это график действительной части одного из корней при изменении параметра t1 и фиксированных всех остальных.

Я вижу глазами как изменение параметров вляет на внешний вид графиков, но нужно понять как он вляет на что-то численное, что определяет этот внешний вид.

Вы уверены, что там именно экспоненты? Для затухания звука - да, моделируется экспонентой, но все-таки, это гармонические колебания же.

// насчет некорректности МНК с экспонентами: вот пример аппроксимации гауссианы. Штука довольно «вредная» иногда параметры подбираются ну совершенно не так, как надо (особенно если данных немного, а шумов предостаточно)

Не уверен, видимо действительно лучше еще подольше глазами посмотерть.

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

А можно вопрос: чем вам не нравится таблично-волновой синтез?

Тут модель строится полностью на физических данных инструмента - длинна струны, сечение, модуль Юнга металла струн. То-есть она очень детально контролируема. Ну и вообще так интереснее.

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

Сложно ведь точно построить модель реального звука. Ну воспроизведете вы первые N гармоник струны, учитывая затухание и распространение волны в реальной среде конечных размеров. Но ведь сама по себе струна не звучит, ей же еще и резонатор нужен! А вот с резонатора и начнется самое веселье.

Хотя, конечно, современные вычислительные мощности вполне могут вам позволить смоделировать звучание гитары или скрипки. Только сомневаюсь, что хватит мощей на realtime.

Послушайте синтезатор с pianoteq.com, у меня знакомая пианистка его от реального инструмента не отличила.

Да, там не все на 100% делается на основе физики. На выходе даже первой ступени уже получается звук инструмента в целом, но те параметры, что можно расчитать - расчитываются.

Потом там хорошая модель изменения звука в зависимости от скорости нажатия (это не только изменение громкости, гармонический состав меняется очень сильно).

При этом хочется получить модель дающую нормальный звук, пусть даже и с «хаками» и постепенно вытеснять их, чем начинать с 0 делать чисто физическую модель. Я попробовал реализовать несколько моделей с научных статей, где все по правилам - получается лажа не описуемая.

Там регистрироваться надо - не хочу.

http://dl.dropbox.com/u/6121480/synth/ref/mozart-fantasy1.ogg (9.5 MB)

http://dl.dropbox.com/u/6121480/synth/ref/bach-air.ogg (7.4 MB)

http://dl.dropbox.com/u/6121480/synth/ref/brahms-hungarian-dance-no.2.ogg (6.0 MB)

http://dl.dropbox.com/u/6121480/synth/ref/chopin-fantaisie-impromptu.ogg (8.4 MB)

Качаю Брамса (т.к. самый короткий). Минут через 20, как скачается, послушаю.

Пример для одной экспоненты (для суммы нескольких – по аналогии) У тебя есть вот такая функция (подчеркивание перед икс – чтобы отделить его от параметров):

f(a1, t1, _x) = a1 * exp (t1 * _x)

Если бы не было коэффициента перед _x, то можно было бы МНК применить – он как раз служит для поиска линейных коэффициентов. Но, раз он есть – все плохо, и простой и надежный среднеквадрат нам не подходит.

Можно сделать такой фокус: рассмотреть это безобразие, как функцию коэффициентов a1, t1. И найти ее минимум - в смысле, минимум отклонения d = (f(x0) - y0). Пометка раз: чтобы уравнять в правах и отрицательные значения, будем искать минимум квадрата отклонения, пометка два – минимум ищем у суммы. res = argmin(sum(d^2)). Это, кстати говоря, уже стандартная задача — «Nonlinear least square optimization».

Уж минимумы искать завсегда умеем, численными методами-то.

Matlab-подобные штуковины даже имеют встроенную функцию для таких мероприятий. Вот, например, копипаста из доков по SciLab для простенькой экспоненты (http://help.scilab.org/docs/5.3.2/en_US/leastsq.html ).

function y=yth(t, x)
   y  = x(1)*exp(-x(2)*t) 
endfunction  

// Отсчеты
m = 10;
tm = [0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5]';
ym = [0.79, 0.59, 0.47, 0.36, 0.29, 0.23, 0.17, 0.15, 0.12, 0.08]';
// Единичные веса
wm = ones(m,1); 

// Минимизируем вот такую функцию
//  f(x) = sum_i  wm(i)^2 ( yth(tm(i),x) - ym(i) )^2   

x0 = [1.5 ; 0.8];

// Функция и ее якобиан
function e=myfun(x, tm, ym, wm)
   e = wm.*( yth(tm, x) - ym )
endfunction

function g=mydfun(x, tm, ym, wm)
   v = wm.*exp(-x(2)*tm)
   g = [v , -x(1)*tm.*v]
endfunction

// Простейший вызов
[f,xopt, gopt] = leastsq(list(myfun,tm,ym,wm),x0)

// С использованием якобиана
[f,xopt, gopt] = leastsq(list(myfun,tm,ym,wm),mydfun,x0)

// Картинки
tt = linspace(0,1.1*max(tm),100)';
yy = yth(tt, xopt);
scf();
plot(tm, ym, "kx")
plot(tt, yy, "b-")
legend(["Отсчеты", "Интерполяция"]);
xtitle("Пример нелинейной интерполяции")

Сейчас нарежу по 30 секунд, не нужно над собой издеваться :)

Да не, скорость подскочила аж до 15КБ/с. Через минутку уже докачается.

Кстати, я здесь соврал. Там никакой не МНК применяется (естественно, т.к. здесь никаким МНК не получится), а итеративный процесс подбора аппроксимирующей функции по критерию хи² (из GSL)

Те же файлы, только с постфиксом «_short.ogg», все в пределах 1 Мб.

Вот это да!!!

Никогда бы не поверил, если бы сам не услышал. Звучит чище и красивее таблично-волнового синтеза.

Оно в режиме реального времени работает, или все-таки дольше считает? И много ли моделей инструментов уже наделали?

Сильно быстрее реального времени, там все распараллелено через MMX и SSE + многопоточность.

Быстрее реального времени работает даже реверс, где все распараллеливание убрано :)

Очень хотелось бы на алгоритмы взглянуть. Я так понимаю, там все сильно огорожено? Публикаций на эту тему никаких нет?

Я из пока не хочу светить, во избежание. Сам синтез - это набор генераторов - простых затухающих по экспоненте синусоид. На каждую гармонику каждой струны таких синусоид 2 или 3 штуки (по числу струн в хоре). Фокус в расчете параметров этих генераторов. Часть этого расчета и есть рассматриваемый алгоритм.

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

Скорость достигается за счет того, что затухающий по экспоненте синус - это БИХ фильтр 2 порядка, с вычисленияим и калькулятор справится :)

Сдается мне, модель все-таки намного более сложная.

500 строчек на Си - это вычисление начальных данных, не зависящих от скорости нажатия + еще 400 - это преобразование этих данных в зависимости от скорости нажатия.

Все остальное кругом очень легко объясняется и логично, если подумать. А вот что делает этоот алогритм - не ясно :)

400 сократятся строк до 200, так как сейчас там еще мешанина из Си и ассемблера переписанного на Си.

Так вы где-то умудрились исходники достать?

Кстати, если там всего лишь 900 строчек - это же вообще халява!

Реверс инжинирингом примерно в течение 3 месяцев всего своего свободного времени и 2 недель отпуска :)

Упорство и труд все перетрут! :)

Желаю удачи!

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

Реверс есть только для синтезатора, в оригинале там еще реверберация и всякая пост обработка добавляется, но это я и так знаю как делать.

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

Вам нужно равномерное приближение.

// тред не читал

В качестве целевой взять корреляционную ф-ю, дальше тупо метод градиентного спуcка?

>Существуют-ли отличные от полного перебора методы поиска аналитических выражений в такого рода ситуациях?

В вариационных методах и получается аналитическое решение в виде разложения по ортонормированному базису полной системы подходящих функций (в том числе и на базисе экспонент).

Но в данной задаче можно попробовать пакет ALGLIB: аппроксимацию линейным или нелинейным МНК, например, интерполяция сплайнами.

> http://dl.dropbox.com/u/6121480/synth/ref/chopin-fantaisie-impromptu.ogg (8.4 MB)

Обалденно! Это синтетика?

Да, из midi файла.

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

Анон, для таких случаев ещё на первом курсе физике советуют прологарифмировать экспоненту,а потом применять МНК. Жаль, что это работает только для одной экспоненты...