[ataradov 0]

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

 

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

 

      a0[0] = a[1] * a[2];
      a0[1] = b[1] * a[2] + a[1] * b[2];
      a0[2] = b[1] * b[2];

      a0[0] += a[0] * a[2];
      a0[1] += b[0] * a[2] + a[0] * b[2];
      a0[2] += b[0] * b[2];

      a0[0] += a[0] * a[1];
      a0[1] += b[0] * a[1] + a[0] * b[1];
      a0[2] += b[0] * b[1];

 

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

 

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

 

[iiv 29]

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

 

(0 1 0)

(0 0 1)

(1 0 0)

 

a и b - векторы, тогда матрица H=

( a, b ) ^T P ( a, b )

на диагонали будет содержать a0 и a2, а сумма внедиагональных членов будет равна a1.

 

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

[ataradov 0]

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

  for (int k = 0; k < N; k++)
  {
    a1[0] = 1.0; a1[1] = a1[2] = a1[3] = 0.0;

    for (int j = 0; j < N; j++)
    {
      if (j == k)
        continue;

      a1_old = 0.0;
      for (int i = 0; i < N+1; i++)
      {
        t = a1[i] * a[j] + a1_old * b[j];
        a1_old = a1[i];
        a1[i] = t;
      }
    }

    for (int i = 0; i < N+1; i++)
      a0[i] += a1[i];
  }

  a1_old = 0.0;
  for (int i = 0; i < N+1; i++)
  {
    t = a1[i] * a[N-2] + a1_old * b[N-1];
    a1_old = a1[i];
    a1[i] = t;
  }

Это кусок мат. модели, векторы a и b вычисляются из физических параметров системы.

 

После этого из массивов а0 и a1 формируются коэфициенты полинома (комплексные):

  for (int i = 0; i < N+1; i++)
  {
    c[i].real= t0 * a1[i];
    c[i].imag = t1 * a1[i];

    if (i > 0)
      c[i].imag -= a0[i-1];
  }

t0 и t1 - тоже вычисляются из параметров, представляют собой выражения вида t0 = A * sin(W), t1 = A * cos(W).

 

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

Изменено 10 октября, 2011 пользователем Taradov Alexander

[iiv 29]

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

[ataradov 0]

a1 инициализирован только в первых 4 элементах. А дальше, в общем случае, тоже нули или модель просто не гоняли?

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

 

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

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

 

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

 

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

  avg = 0.0;
  for (int i = 0; i < N; i++)
    avg += arr[i];
  avg = avg / N;

  a[N] = b[N] = 0.0;
  for (int i = 0; i < N; i++)
  {
    a[i] = (2.0 - arr[i] / avg) / (2.0 * avg);
    b[i] = (1.0 - arr[i] / avg) * overtone * M_PI;
  }

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

 

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

 

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

 

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

Но результат - напрямую частоты и экспоненты, говорит о том, что используются какие-то свойства этого полинома и его корней.

 

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

Изменено 10 октября, 2011 пользователем Taradov Alexander

[ataradov 0]

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

 

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

[ataradov 0]

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

arr = [261.4508114 261.5063567 261.6451613];
avg = sum(arr) / length(arr);

a = (2.0 - arr / avg) / (2.0 * avg);
b = (1.0 - arr / avg) * (1) * pi;

az1 = conv([1.0, 0.0], [a(2) b(2)]);
az1 = conv(az1, [a(3) b(3)]);

az2 = conv([1.0, 0.0], [a(1) b(1)]);
az2 = conv(az2, [a(3) b(3)]);

az3 = conv([1.0, 0.0], [a(1) b(1)]);
az3 = conv(az3, [a(2) b(2)]);

a0 = az1+az2+az3;

a1 = conv(az3, [a(2) b(3)]);

tt0 = 222.9164178;
tt1 = -1.571716919;

for i = 1:4,
  z(i) = tt0 * a1(i) + j*(tt1 * a1(i));
  if i > 1,
    z(i) = z(i) - j*a0(i-1);
  end;
end;

r = roots(z);

 

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

[ataradov 0]

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

[iiv 29]

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

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

 

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

 

С уважением

 

ИИВ

[ataradov 0]

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

 

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

 

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

clear all;

arr = [261.4508114 261.5063567 261.6451613];
avg = sum(arr) / length(arr);

t = 1 - arr / avg;
a = (t+1) / (avg * 2);
b = t * (1) * pi;

az1 = conv(conv([1.0, 0.0], [a(2) b(2)]), [a(3) b(3)]);
az2 = conv(conv([1.0, 0.0], [a(1) b(1)]), [a(3) b(3)]);
az3 = conv(conv([1.0, 0.0], [a(1) b(1)]), [a(2) b(2)]);

a0 = az1+az2+az3;
a0 = [0 a0];

a1 = conv(az3, [a(2) b(3)]);

t0 = 199.1735851;
t1 = -0.00705058376;

z = t0 * a1(1:4) * exp(j*t1) - j*a0(1:4);

r = roots(z);

В результате

r =

  -0.0568 + 7.8407i
   0.3645 + 0.0279i
  -0.3631 + 0.0056i

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

 

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

 

[iiv 29]

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

 

Если представить $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 матрица сдвига, как я сказал выше.

 

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

[ataradov 0]

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

 

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

[ataradov 0]

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

 

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

 

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