Избегание арифметики с плавающей запятой

Question

Я написал небольшой программный синтезатор для iPhone.
Чтобы дополнительно настроить производительность, я измерил свое приложение с помощью Shark и обнаружил, что теряю много времени при преобразованиях float / SInt16.
Поэтому я переписал некоторые части, чтобы обойти преобразования, предварительно рассчитав таблицы поиска, которые возвращают "готовые к использованию" образцы SInt16.Пока это работает нормально.
В настоящее время я пытаюсь переписать некоторые фильтры и мою реализацию ADSR envelope, чтобы использовать только целочисленную арифметику, но я мог бы использовать некоторые советы, как выполнять умножения / деления без чисел с плавающей точкой.
Я нацелен на Канонический формат iPhone:

• LPCM
• 16-битные целочисленные выборки

Каковы хорошие подходы для применения амплитуды к моему окончательному образцу без использования float?

Редактировать:
Единственное, что я выяснил на данный момент, это то, что я могу разделить на степени 2, сдвинув вправо мой текущий образец.

inBuffer[frame] = wavetable[i % cycleLengthInSamples] >> 4;

Но я не могу придумать ни одного элегантного способа создать с помощью этого гладкую оболочку ADSR.

Редактировать 2: Спасибо за все ваши замечательные ответы!
Мой нынешний подход:

• приведите все мои значения огибающей ADSR к положительному диапазону SInt16
• умножьте на текущее значение из wavetable (сохраните промежуточные значения как SInt32).
• сдвиньте результат на 16 вправо

кажется, это работает :)

Answer 1

Фиксированная точка — это хорошо, поскольку в этом случае вы используете 16 бит.Самый простой способ — умножить на 10 в зависимости от необходимой вам точности.Если вы можете использовать 32-битные целые числа в качестве промежуточного звена, вы сможете получить приличную точность.В конце вы можете преобразовать обратно в 16-битное целое число, округляя или усекая его по своему усмотрению.

Редактировать:Вам нужен сдвиг влево, чтобы увеличить значения.Сохраните результат сдвига в виде с большей точностью (32 или 64 бита в зависимости от того, что вам нужно). простой сдвиг не будет работать, если вы используете подписанные типы

Будьте осторожны, если вы умножаете или делите два числа с фиксированной точкой.Умножение получается (a*n) * (bn), и вы получитебn^2 вместо aбн.Дивизия – это (ап) / (бn), что представляет собой (a/b) вместо ((aп)/б).Вот почему я предложил использовать степени 10, чтобы было легче найти ошибки, если вы не знакомы с фиксированной точкой.

Когда вы закончите свои вычисления, вы сдвинетесь вправо, чтобы вернуться к 16-битному int.Если вы хотите пофантазировать, вы также можете выполнить округление перед сдвигом.

Я предлагаю вам почитать, если вы действительно заинтересованы в реализации эффективной фиксированной точки. http://www.digitalsignallabs.com/fp.pdf

Answer 2

Ответы на это ТАК вопрос довольно обширны с точки зрения реализации.Вот немного больше объяснений, чем я видел там:

Один из подходов — поместить все ваши числа в диапазон, скажем, [-1.0,1.0).Затем вы сопоставляете эти числа с диапазоном [-2^15,(2^15)-1].Например,

Half = round(0.5*32768); //16384
Third = round((1.0/3.0)*32768); //10923

Умножив эти два числа, вы получите

Temp = Half*Third; //178962432
Result = Temp/32768; //5461 = round(1.0/6.0)*32768

Делить на 32768 в последней строке и есть точка Патрос Сделал умножение, требующее дополнительного шага масштабирования.Это имеет больше смысла, если вы явно запишете масштабирование 2^N:

x1 = x1Float*(2^15);
x2 = x2Float*(2^15);
Temp = x1Float*x2Float*(2^15)*(2^15);
Result = Temp/(2^15); //get back to 2^N scaling

Итак, это арифметика.Обратите внимание, что для реализации умножение двух 16-битных целых чисел требует 32-битного результата, поэтому Temp должен быть 32-битным.Кроме того, число 32768 невозможно представить в виде 16-битной переменной, поэтому имейте в виду, что компилятор сразу создаст 32-битные переменные.И, как вы уже заметили, вы можете перейти к умножению/делению на степени 2, чтобы вы могли написать

N = 15;
SInt16 x1 = round(x1Float * (1 << N));
SInt16 x2 = round(x2Float * (1 << N));
SInt32 Temp = x1*x2;
Result = (SInt16)(Temp >> N);
FloatResult = ((double)Result)/(1 << N);

Но предположим, что [-1,1) не является правильным диапазоном?Если вы предпочитаете ограничить свои числа, скажем, [-4.0,4.0), вы можете использовать N = 13.Тогда у вас есть 1 знаковый бит, два бита до двоичной точки и 13 после.Они называются дробными типами с фиксированной запятой 1,15 и 3,13 соответственно.Вы обмениваете точность дроби на запас.

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

Quotient = (x1/x2) << N;

или, для сохранения точности

Quotient = (SInt16)(((SInt32)x1 << N)/x2); //x1 << N needs wide storage

Умножение и деление на целые числа работает нормально.Например, чтобы разделить на 6, вы можете просто написать

Quotient = x1/6; //equivalent to x1Float*(2^15)/6, stays scaled

А в случае деления на степень 2,

Quotient = x1 >> 3; //divides by 8, can't do x1 << -3 as Patros pointed out

Однако сложение и вычитание целых чисел не работает наивно.Сначала вам нужно посмотреть, соответствует ли целое число вашему типу xy, создать эквивалентный дробный тип и продолжить.

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

Answer 3

Взгляните на эту страницу, где описаны алгоритмы быстрого умножения.

http://www.newton.dep.anl.gov/askasci/math99/math99199.htm

Answer 4

В общем, допустим, вы будете использовать знаковое представление с фиксированной точкой 16.16.Таким образом, 32-битное целое число будет иметь 16-битную целочисленную часть со знаком и 16-битную дробную часть.Тогда я не знаю, какой язык используется при разработке iPhone (возможно, Objective-C?), но этот пример написан на C:

#include <stdint.h>

typedef fixed16q16_t int32_t ;
#define FIXED16Q16_SCALE 1 << 16 ;

fixed16q16_t mult16q16( fixed16q16_t a, fixed16q16_t b )
{
    return (a * b) / FIXED16Q16_SCALE ;
}

fixed16q16_t div16q16( fixed16q16_t a, fixed16q16_t b )
{
    return (a * FIXED16Q16_SCALE) / b ;
}

Обратите внимание, что приведенное выше является упрощенной реализацией и не обеспечивает защиты от арифметического переполнения.Например, в div16q16() я умножаю перед разделением, чтобы сохранить точность, но в зависимости от операндов операция может переполниться.Вы могли бы использовать 64-разрядный промежуточный сервер, чтобы преодолеть это.Кроме того, деление всегда округляется в меньшую сторону, потому что оно использует целочисленное деление.Это обеспечивает наилучшую производительность, но может повлиять на точность итерационных вычислений.Исправления просты, но увеличивают накладные расходы.

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

На языке OO fixed16q16_t, естественно, был бы кандидатом для класса с перегрузкой операторов, поэтому вы можете использовать его как обычный арифметический тип.

Возможно, вам будет полезно конвертировать между типами:

double fixed16q16_to_double( fixed16q16_t fix )
{
    return (double)fix / FIXED16Q16_SCALE ;
}

int fixed16q16_to_int( fixed16q16_t fix )
{
    // Note this rounds to nearest rather than truncates
    return ((fix + FIXED16Q16_SCALE/2)) / FIXED16Q16_SCALE ;
}

fixed16q16_t int_to_fixed16q16( int i )
{
    return i * FIXED16Q16_SCALE ;
}

fixed16q16_t double_to_fixed16q16( double d )
{
    return (int)(d * FIXED16Q16_SCALE) ;
}

Это основы, можно сделать более сложными и добавить тригонометрические и другие математические функции.

Исправлено сложение и вычитание, которые работают со встроенными операторами + и - и их вариантами.