c - 使用相位累加器的 FM 合成

Question

我正在尝试使用 C 中的相位累加器实现带反馈的 FM 合成运算符。在 Tomisawa 的 original patent 中，进入加法器的相位累加器对负索引和正索引进行计数，从 -pi 正弦波相位的 -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

首先我们可以尝试用 C 语言编写 pyton 代码

#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);
    }
}

下一步 - 使用 sin 的计算值

#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 值和输出。在此“输出”中的版本值如 XXXXXXQQ.IIIIIIIII.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);
    }
}