FM синтез с использованием фазового аккумулятора

Question

Я пытаюсь реализовать оператор синтеза FM с обратной связью, используя фазовый аккумулятор в C. В оригинальном патенте Tomisawa, фазовый аккумулятор, поступающий в сумматор, учитывает как отрицательные, так и положительные индексы, начиная с -2^(n-1} на фаза синусоидальной волны от -pi до 2^(n-1) в фазе пи. Для простоты я хотел бы использовать фазовый аккумулятор, который считает только положительные значения, используя верхние байты 32-разрядного целого числа без смещения в качестве индекса в поиске таблицы синусов.

Я экспериментировал с этим, и, к сожалению, я не могу получить алгоритм для получения ожидаемых результатов при использовании обратной связи. Добавление выхода синусоидальной волны к фазовому аккумулятору должно привести к образованию пилообразной формы волны, но я не могу понять, как правильно добавить выходную синусоидальную волну (которая представляет собой 16-битовое целое число со знаком) к аккумулятору фазы без знака, чтобы произвести это. Мы ценим любые предложения.

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

Некоторое разъяснение, вероятно, в порядке. Вот некоторые диаграммы из оригинального патента Tomisawa:

Когда и фазовый аккумулятор, и выход синусоидальной волны подписаны, алгоритм достаточно прост для реализации. Аккумулятор фазы начинается с -1 и продолжается до 1, а выход синусоидальной волны также находится между -1 и 1. В Python алгоритм выглядит примерно так, чтобы генерировать 1000 выборок:

table = []
feedback = 0.25
accumulator = -1

for i in xrange(1000):
    output = math.sin(math.pi*(accumulator + feedback*output)
    table[i] = output
    accumulator += 0.005
    if accumulator > 1:
        accumulator = -1

Который производит вывод, который выглядит следующим образом:

Я пытаюсь адаптировать этот алгоритм к C. В C, в целях вычислительной эффективности, я хотел бы, чтобы фазовым аккумулятором было 32-разрядное целое число без знака, а не целое число со знаком. Таким образом, я могу использовать первые два бита старшего байта аккумулятора в качестве индекса квадранта, а второй старший байт в качестве индекса в массиве 256 16-битных значений синусов для таблицы синусов 1024 значений. Подобно:

XXXXXXQQ.IIIIIIII.XXXXXXXX.XXXXXXXX
      ^^ ^^^^^^^^
 quadrant    index

Моя проблема в том, что у меня возникают трудности с адаптацией алгоритма FM, который задан для неподписанного фазового аккумулятора. Если фазовый накопитель представляет собой 32-разрядное целое число без знака, а выходная таблица синусоидальной волны представляет собой 16-разрядное целое число (со знаком или без знака), как я могу адаптировать алгоритм, как показано в патенте и коде Python выше, для работы с этим форматом, и производить тот же вывод?

Answer 1

Прежде всего мы можем попытаться написать вам код Python на C

#include <stdio.h>
#include <math.h>

void main() {
    double table[1000];
    double feedback = 0.25;
    double accumulator = -1;

    int i;
    for (i=0;i<1000;i++) {
        double output = sin(M_PI*(accumulator + feedback*output));
        table[i]=output;
        accumulator += 0.005;
        if (accumulator > 1)
            accumulator = -1;
        printf("%f\n",output);
    }
}

Следующий шаг - использовать рассчитанные значения для греха

#include <stdio.h>
#include <math.h>

void main() {
    double table[1000];
    double feedback = 0.25;
    double accumulator = 1;

    int i;

    double sinvalue[1024];
    for (i=0;i<1024;i++) {
        sinvalue[i]=sin(M_PI*i/512);
    }

    for (i=0;i<1000;i++) {
        double output = sinvalue[(int)(512*(accumulator + feedback*output))%1024];
        printf("%0.6f %0.6f %0.6f\t",accumulator,feedback,output);
        table[i]=output;
        accumulator += 0.005;
        if (accumulator > 2)
            accumulator = 0;
        printf("%f\n",output);
    }
}    

Следующий шаг - использование 16-битных значений sin и output. В этой версии значение в "выводе", например, XXXXXXQQ.IIIIIIII.XXXXXXXX.XXXXXXXX Также мы потеряли некоторую точность.

#include <stdio.h>
#include <math.h>

#define ONE ((int)(2*256*256*256/M_PI))

void main() {
    double table[1000];
    double feedback = 0.25;
    double accumulator = 1;
    double accumulatorDelta = 0.005;

    unsigned int feedback_i = ONE*feedback/32768;
    unsigned int accumulator_i = ONE*accumulator;
    unsigned int accumulatorDelta_i = ONE*accumulatorDelta;

    int i;

    double sinvalue[1025];
    short int sinvalue_i[1025];
    for (i=0;i<1025;i++) {
        sinvalue[i]=sin(M_PI*i/512);
        sinvalue_i[i]=32786*sinvalue[i];
        if (sinvalue[i]*32768>32768) sinvalue_i[i]=32768;
        if (sinvalue[i]*32768<-32767) sinvalue_i[i]=-32767;
    }

    for (i=0;i<1000;i++) {

        double output = sin(M_PI*(accumulator + feedback*output));
        short int output_i = sinvalue_i[ ((unsigned int) ((accumulator_i + feedback_i*output_i)*M_PI)>>16)%1024 ];
        table[i]=output;

        accumulator += 0.005;
        if (accumulator > 2)
            accumulator = 0;

        accumulator_i += accumulatorDelta_i;
        if (accumulator_i > 2*ONE)
            accumulator_i = 0;

        printf("%f %f %04X\n",output,(float)output_i/32768,(unsigned short int)output_i);
    }
}

Но мы потеряли некоторое время на преобразование int->double->int Если мы изменим ОДНУ константу, мы упустим возможность быстро получить квадрант, но избавимся от преобразования

#include <stdio.h>
#include <math.h>

#define ONE ((int)(2*256*256*256))

void main() {
    short int table[1000];

    unsigned int feedback_i = ONE*0.25/32768;
    unsigned int accumulator_i = ONE*1;
    unsigned int accumulatorDelta_i = ONE*0.005;

    int i;

    short int sinvalue_i[1025];
    for (i=0;i<1025;i++) {
        double sinvalue=sin(M_PI*i/512);
        sinvalue_i[i]=32786*sinvalue;
        if (sinvalue*32768>32768) sinvalue_i[i]=32768;
        if (sinvalue*32768<-32767) sinvalue_i[i]=-32767;
    }

    for (i=0;i<1000;i++) {
        short int output_i = sinvalue_i[ ( (accumulator_i + feedback_i*output_i)>>16)%1024 ];
        table[i]=output_i;

        accumulator_i += accumulatorDelta_i;
        if (accumulator_i > 2*ONE)
            accumulator_i = 0;

        printf("%f %04X\n",(float)output_i/32768,(unsigned short int)output_i);
    }
}