Быстрый эквивалент sin () для DSP, указанного в STK - PullRequest
3 голосов
/ 15 февраля 2012

Я использую кусочки Perry Cook's Synthesis Toolkit (STK) для генерации пильных и прямоугольных волн.STK включает в себя этот пилообразный генератор на базе BLIT:

inline STKFloat BlitSaw::tick( void ) {
  StkFloat tmp, denominator = sin( phase_ );
  if ( fabs(denominator) <= std::numeric_limits<StkFloat>::epsilon() )
      tmp = a_;
  else {
      tmp = sin( m_ * phase_ );
      tmp /= p_ * denominator;
  }

  tmp += state_ - C2_;
  state_ = tmp * 0.995;

  phase_ += rate_;
  if ( phase_ >= PI ) 
     phase_ -= PI;

  lastFrame_[0] = tmp;
     return lastFrame_[0];
}

Генератор прямоугольных импульсов очень похож.Вверху есть такой комментарий:

// A fully  optimized version of this code would replace the two sin 
// calls with a pair of fast sin oscillators, for which stable fast 
// two-multiply algorithms are well known.

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

#define AD_SIN(n) (n*(2.f- fabs(n))) 

Построение графика показывает, что это не совсем близкое приближение вне диапазона от -1 до 1, поэтому я недумаю, что я могу использовать его, когда phase_ находится в диапазоне от -pi до pi:

Plot of Sine vs. approximation

Здесь синус - синяя линия, а фиолетовая линия - приближение.

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

Спасибо

РЕДАКТИРОВАТЬ Спасибо за подробные и разнообразные ответы.Я изучу их и приму один на выходных.

РЕДАКТИРОВАТЬ 2 Не мог бы анонимный близкий избиратель, пожалуйста, разъяснить свой голос в комментариях?Спасибо.

Ответы [ 7 ]

6 голосов
/ 15 февраля 2012

По сути, синусоидальный генератор представляет собой одну (или более) переменную, которая изменяется с каждым шагом DSP, а не пересчитывается с нуля.

Простейшие основаны на следующих триггерных тождествах: (где d является константой, и, таким образом, также cos(d) и sin(d))

sin(x+d) = sin(x) cos(d) + cos(x) sin(d)
cos(x+d) = cos(x) cos(d) - sin(x) sin(d)

Однако для обновления требуется две переменные (одна для греха и одна для cos) и 4 умножения. Однако это все равно будет намного быстрее, чем вычисление полного синуса на каждом шаге.

Решение Оли Чарльзворта основано на решениях этого общего уравнения

A_{n+1} = a A_{n} + A_{n-1}

Где поиск решения вида A_n = k e^(i theta n) дает уравнение для theta.

e^(i theta (n+1) ) = a e^(i theta n ) + b e^(i theta (n-1) )

Что упрощает до

e^(i theta) - e^(-i theta ) = a
2 cos(theta) = a

Предоставление

A_{n+1} = 2 cos(theta) A_{n} + A_{n-1}

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

3 голосов
/ 15 февраля 2012

Насколько вам это нужно?

Эта функция, f (x) = 0.398x * (3.1076- | x |), выполняет достаточно хорошую работу для x между -pi и pi.

Редактировать

Еще лучшим приближением является f (x) = 0,38981969947653056 * (pi- | x |), в котором абсолютная ошибка сохраняется равной 0,038158444604 или меньше для x между -пи и пи.

Минимизация по методу наименьших квадратов даст немного другую функцию.

2 голосов
/ 15 февраля 2012

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

Например, учтите, что следующее разностное уравнение даст вам синусоиду:

y[n] = 2*cos(phi)*y[n-1] - y[n-2]

(где cos(phi) - это константа)

1 голос
/ 13 мая 2013

Вы можете посмотреть здесь:

http://devmaster.net/forums/topic/4648-fast-and-accurate-sinecosine/

Есть некоторый пример кода, который вычисляет очень хорошую аппроксимацию sin / cos, используя только умножения, сложения и abs () функция.Довольно быстро тоже.Комментарии также хорошо читаются.

Это сводится к следующему:

float sine(float x)
{
    const float B = 4/pi;
    const float C = -4/(pi*pi);
    const float P = 0.225;    

    float y = B * x + C * x * abs(x);

    return P * (y * abs(y) - y) + y;
}

и работает в диапазоне от -PI до PI

1 голос
/ 16 марта 2013

(От первоначального автора кода VST BLT).

На самом деле я переносил генераторы VST BLT на C #, поэтому я искал хорошие синус-генераторы.Вот что я придумал.Перевод на C ++ прост.См. Заметки в конце о накопленных ошибках округления.

public class FastOscillator
{

    private double b1;
    private double y1, y2;

    private double fScale;

    public void Initialize(int sampleRate)
    {
        fScale = AudioMath.TwoPi / sampleRate;
    }
     // frequency in Hz. phase in radians.
    public void Start(float frequency, double phase)
    {
        double w = frequency * fScale;
        b1 = 2.0 * Math.Cos(w);
        y1 = Math.Sin(phase - w);
        y2 = Math.Sin(phase - w * 2);
    }
    public double Tick()
    {
        double y0 = b1 * y1 - y2;
        y2 = y1;
        y1 = y0;
        return y0;
    }
}

Обратите внимание, что эта конкретная реализация генератора будет дрейфовать со временем, поэтому ее необходимо периодически переинициализировать.В этой конкретной реализации величина синусоиды уменьшается со временем.В первоначальных комментариях в коде STK предлагался генератор с двумя множителями.Фактически существуют два умножителя, которые достаточно стабильны во времени.Но в ретроспективе необходимость синхронизации синусоидальных (фазовых) и синусоидальных (m * фазовых) осцилляторов, вероятно, означает, что в любом случае их необходимо повторно синхронизировать.Ошибки округления между фазой и фазой m * означают, что, даже если осцилляторы были стабильными, они в конечном итоге дрейфовали, создавая значительный риск создания больших всплесков значений вблизи нулей функций BLT.Можно также использовать генератор с одним множителем.

Эти конкретные генераторы, вероятно, следует повторно инициализировать каждые 30–100 циклов (или около того).Моя реализация C # основана на фреймах (то есть вычисляет массив результатов float [] в методе void Tick (int count, float [] result). Осцилляторы повторно синхронизируются в конце каждого вызова Tick. Примерно так:

   void Tick(int count, float[] result)   
   {
        for (int i = 0; i < count; ++i)
        {
                ...
            result[i] = bltResult;
        }
        // re-initialize the oscillators to avoid accumulated drift.
        this.phase = (this.phase + this.dPhase*count) % AudioMath.TwoPi;
        this.sinOsc.Initialize(frequency,this.phase);
        this.mSinOsc.Initialize(frequency*m,this.phase*m);
   }

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

Существует своеобразный взлом, который предотвращает дрейф постоянного тока осцилляторов, даже когда происходит качание частоты генераторов.является то, что генераторы должны быть запущены с начальной настройкой фазы dPhase / 2. Это так же, как и для запуска генераторов с нулевым дрейфом постоянного тока, без необходимости выяснять правильное начальное состояние для различных интеграторов в каждом из генераторов BLT.

Странно, еслирегулировка переустанавливается всякий раз, когда изменяется частота генератора, тогда это также предотвращает дикий дрейф постоянного тока на выходе при качании частоты генератора.Всякий раз, когда частота изменяется, вычтите dPhase / 2 из предыдущего значения фазы, пересчитайте dPhase для новой частоты, а затем добавьте dPhase / 2. Я подозреваю, что это может быть формально доказано;но я так не смог.Все, что я знаю, это то, что это просто работает.

Для реализации блока осцилляторы должны фактически инициализироваться следующим образом, вместо переноса регулировки фазы в текущем значении this.phase.

        this.sinOsc.Initialize(frequency,phase+dPhase*0.5);
        this.mSinOsc.Initialize(frequency*m,(phase+dPhase*0.5)*m);
0 голосов
/ 16 февраля 2012

Общая идея получения периодически дискретизированных результатов с помощью функции синуса или косинуса заключается в использовании рекурсии триггера или инициализированного (едва ли) стабильного БИХ-фильтра (который может в конечном итоге быть почти такими же вычислениями).В литературе, посвященной DSP, есть такие группы, с различной точностью и стабильностью.Тщательно выбирайте.

0 голосов
/ 15 февраля 2012

Если вы можете, вы должны рассмотреть методы, основанные на запоминании.По существу, храните значения sin (x) и cos (x) для значений группы.Чтобы вычислить sin (y), найдите a и b, для которых существуют предварительно вычисленные значения, такие что a <= y <= b.Теперь используя sin (a), sin (b), cos (a), cos (b), ya и yb приблизительно вычислим sin (y). </p>

...