Базовый алгоритм

Основной алгоритм обычно строится из предположения, что время перехода волны через границы секций происходит в моменты времени (T - t) и (T + t), исходя из чего период t выбирается в 2 раза меньше, чем необходимо для пробега звуковой волны. Это, в свою очередь, приводит к работе алгоритма на излишне высоких частотах.

Частота дискретизации рассчитывается по формуле:

Fd = C / dX

где:

Fd - частота дискретизации, Гц

C - скорость звука, м/с

dX - длина секции, м

При обычно используемой величине dX около 0.4 см частота дискретизации требуется около 80 кГц.

Алгоритм расчета для волн давления следующий:

1.Прочитать площади секций в массив A размерности N.

2.Рассчитать коэффициенты отражения Ri для i = {0,1,2, . . . , N-2}
R(i) = ( A(i) - A(i + 1) ) / ( A(i) + A(i + 1) )

3.Установить коэффициенты отражения со стороны голосовой щели и губ: Rg = +0.99, Rl = -0.99.

4.Создать массивы F(N), B(N) и F1(N), B1(N) и инициализировать их нулями

5.Установить текущее число сэмплов n = 0

6.Пока n <= числа сэмплов звука, повторять шаги 7 - 12.

7.F(0) установить равным сумме отражённой волны и текущего значения давления голосового источника, u(n) - объёмная скорость
F(0) = u(n) * ρ * Cx / A(0) + Rg * B(0)

8.Для всех i = {0, 1, 2, ..., N-2}, выполнить расчет переходов:
d = Ri * [F1(i) - B1(i + 1)]
F(i + 1) = F1(i) + d
B(i) = B1(i + 1) + d

9.Отразить прямую компоненту волны от губ назад:
b(N - 1) = Rl * F(N - 1)

10.Рассчитать и запомнить давление на губах и запомнить значение в массив Pout.
Pout(n) = F(N - 1) + B(N - 1)

11.Переписать значения из массивов F(N), B(N) в F1(N), B1(N)

12.Увеличить номер сэмпла n на 1. Перейти к пункту 6

Пример синтеза можно найти, например, здесь: http://www.falstad.com/vowel/.

Ниже приводится модифицированный алгоритм, позаимствованный из работы Siddharth Mathur, VARIABLE-LENGTH VOCAL TRACT MODELING FOR SPEECH SYNTHESIS.

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

Частота дискретизации рассчитывается по формуле:

Fd = C / 2 / dX

где:

Fd - частота дискретизации, Гц

C - скорость звука, м/с

dX - длина секции, м

При обычно используемой величине dX около 0.4 см, частота дискретизации требуется около 40 кГц.

Алгоритм расчета для волн давления следующий:

1.Прочитать площади секций в массив A размерности N.

2.Установить номер индекса массива для которого будет считаться излучение Radiation = N - 1

3.Если число элементов массива нечётное, добавить еще один элемент и присвоить ему значение A(N - 1) = A(N - 2). Размерность массива теперь увеличилась на 1. N = N + 1

4.Рассчитать коэффициенты отражения Ri для i = {0,1,2, . . . , N-2}
R(i) = ( A(i) - A(i + 1) ) / ( A(i) + A(i + 1) )

5.Установить коэффициенты отражения со стороны голосовой щели и губ: Rg = +0.99, Rl = -0.99.

6.Создать массивы F(N) и B(N) и инициализировать их нулями

7.Установить текущее число сэмплов n = 0

8.Пока n <= числа сэмплов звука, повторять шаги 9 - 14.

9.F(0) установить равным сумме отражённой волны и текущего значения давления голосового источника, u(n) - объёмная скорость
F(0) = u(n) * ρ * Cx / A(0) + Rg * B(0)

10.Для всех i = {1, 3, 5, ..., N-3}, выполнить расчет переходов:
d = Ri * [F(i) - B(i + 1)]
F(i + 1) = F(i) + d
B(i) = B(i + 1) + d

11.Для всех i = {0, 2, 4, ..., N-2}, выполнить расчет переходов:
d = Ri * [F(i) - B(i + 1)]
F(i + 1) = F(i) + d
B(i) = B(i + 1) + d

12.Отразить прямую компоненту волны от губ назад:
b(Radiation) = Rl * F(Radiation)

13.Рассчитать и запомнить давление на губах и запомнить значение в массив Pout.
Pout(n) = F(Radiation) + B(Radiation)

14.Увеличить номер сэмпла n на 1. Перейти к пункту 8

Код

 public static final double C0 = 330;//скорость звука

 public static final double Ro = 0.00114;//плотность воздуха

 public static double calcDX(int sampleRate) {

         return (double)C0 / 2.0 / (double)sampleRate;

 }

 public static int calcSampleRate(double DX) {

         return (int)(C0 / 2.0 / DX);

 }

 private double calcPresure(double U, double A) {

         if( A > 0 )//устранить неопределённость деления на 0

                 return U*Ro*C0/A;

         else

                 return 10000;

 }

 public void calcRWbase(double[] pArea, double[] snd, int pos, int sampleCount)

 {

         int radiation = pArea.length - 1;

         int N = pArea.length + (pArea.length & 1);

         //

         double[] pA = new double[N];

         System.arraycopy( pArea, 0, pA, 0, pArea.length);

         //создание массивов и инициализация их нулями

         double[] pF = new double[N];

         double[] pB = new double[N];

         //

         N--;

         pA[N] = pA[N - (pArea.length & 1)];

         //calc R

         double[] pR = new double[N];

         for( int i = 0; i < N; i++ ) {

                 pR[i] = (pA[i] - pA[i + 1]) / (pA[i + 1] + pA[i]);

         }

         //

         for( int n = 0; n < sampleCount; n++ ) {

                 double U = getUpulse();

                 pF[0] = calcPresure( U, pA[0] ) + 0.99 * pB[0];//glottal

                 for( int i = 1; i < N; i += 2 ) {

                         double d = pR[i]*(pF[i] - pB[i + 1]);

                         pF[i + 1] = pF[i] + d;

                         pB[i] = pB[i + 1] + d;

                 }

                 for( int i = 0; i < N; i += 2 ) {

                         double d = pR[i]*(pF[i] - pB[i + 1]);

                         pF[i + 1] = pF[i] + d;

                         pB[i] = pB[i + 1] + d;

                 }

                 pB[radiation] = -0.99 * pF[radiation];//lips

                 double P = pB[radiation] + pF[radiation];

                 snd[pos++] = P;

         }

 }