FM synthesis using phase accumulator

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:

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:

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?

Answer 1

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