我正在尝试使用C中的相位累加器来实现具有反馈的FM合成运算符.在Tomisawa的
original patent中,进入加法器的相位累加器计数负指针和正指数,从正弦的-2 ^(n-1}开始计数在pi的一个相位上-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位int,并且正弦波表输出是(有符号或无符号)16位整数,我如何调整上述专利和Python代码中所示的算法来处理这种格式,并产生相同的输出?
最佳答案 首先,我们可以尝试在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.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);
}
}