Question

I'm attempting to implement an FM synthesis operator with feedback using a phase accumulator in C. In Tomisawa's original patent, the phase accumulator going into the adder counts over both negative and positive indices, from -2^(n-1} at a sine wave phase of -pi to 2^(n-1) at a phase of pi. For simplicity's sake, I'd like to use a phase accumulator that counts only over positive values, using the top bytes of an unsinged 32 bit integer as an index into a sine table lookup.

I have experimented with this and unfortunately I can't seem to get the algorithm to produce the expected results when using feedback. Adding the sine wave output to the phase accumulator should produce a sawtooth like waveform, but I can't figure out how to properly add the output sine wave (which is a 16 bit signed int) to the unsigned phase accumulator to produce this. Any suggestions would be appreciated.

Edit:

Some clarification is probably in order. Here are some diagrams from Tomisawa's original patent:

algorithm

output

When both the phase accumulator and the sine wave output are signed, the algorithm is straightforward enough to implement. The phase accumulator starts at -1 and runs to 1, and the sine wave output is also between -1 and 1. In Python, the algorithm looks sort of like this to generate 1000 samples:

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

Which produces an ouput that looks as follows:

output

I am trying to adapt this algorithm to C. In C, for purposes of computational efficiency, I would like the phase accumulator to be a 32 bit unsigned integer rather than a signed integer. That way, I can use the first two bits of the high byte of the accumulator as a quadrant index, and the second high byte as an index into an array of 256 16 bit sine values, for a 1024 value sine table. Like:

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

My problem is that I am having difficulty adapting the FM algorithm as given to an unsigned phase accumulator. If the phase accumulator is an unsigned 32 bit int, and the sine wave table output is a (signed or unsigned) 16 bit integer, how can I adapt the algorithm as shown in the patent and the Python code above to work with this format, and produce the same output?

Was it helpful?

Solution

First of all we can try to write you pyton code on 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);
    }
}

Next step - use calculated values for 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);
    }
}    

Next step - use 16-bits values of sin and output. In this version value in "output" like XXXXXXQQ.IIIIIIII.XXXXXXXX.XXXXXXXX Also, we lost some accuracy.

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

But we lost some time for conversion int->double->int If we change ONE constant we well lose opportunity to fast getting quadrant, but get rid from conversion

#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);
    }
}    
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top