4

我正在尝试使用 C 中的相位累加器来实现一个带有反馈的 FM 合成算子。在 Tomisawa 的原始专利中,进入加法器的相位累加器对负数和正数都进行计数,从 -2^(n-1} 在-pi 的正弦波相位在 pi 的相位处变为 2^(n-1)。为简单起见,我想使用一个相位累加器,它只对正值进行计数,使用一个未分类的 32 位整数的顶部字节作为正弦表查找的索引。

我已经对此进行了实验,不幸的是,在使用反馈时,我似乎无法让算法产生预期的结果。将正弦波输出添加到相位累加器应该会产生类似锯齿的波形,但我不知道如何正确地将输出正弦波(这是一个 16 位有符号整数)添加到无符号相位累加器来产生这个。任何建议,将不胜感激。

编辑:

一些澄清可能是为了。以下是富泽的原始专利中的一些图表:

算法

输出

当相位累加器和正弦波输出都被签名时,算法很容易实现。相位累加器从 -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 代码中所示的算法以使用这种格式,并产生相同的输出?

4

1 回答 1

2

首先,我们可以尝试在 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);
    }
}    

下一步 - 使用 sin 和 output 的 16 位值。在这个版本中,“输出”中的值就像 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);
    }
}    
于 2013-06-04T11:27:01.230 回答