// Fractional-N Frequency Synthesizer
// MASH Simulation
// www.holmea.demon.co.uk

#include <stdio.h>
#include <string.h>
#include <math.h>

const double PI = 3.14159265358979323846;

#define MODULO (1<<16)      // Modulus of overflowing accumulators
#define SEQLEN (4*MODULO)   // Length of sequence (2* for F even; 4* for F odd)

// Time represented as multiples of the carrier (VCO) period
double  tREF[SEQLEN];       // Position of REF edge
int     tVCO[SEQLEN];       // Position of divided VCO edge (always an integer)

// Sigma-delta sequence
char Delta[SEQLEN];

// Parameters
const int N=17, F=1;


class MASH
{
    int q[4];       // Latch outputs
    int c[4][4];    // Carry flip flops

public:
    MASH() // Constructor
    { 
        // Initial values
        memset(c, 0, sizeof c);
        q[0] = 1;
        q[1] = 0;
        q[2] = 0;
        q[3] = 0;

        // Fill flip-flops
        for (int i=0; i<4; i++) Clock();
    }

    int Clock()
    {
        int d[4]; // Adder output / latch data input

        // Carry flip-flops
        c[3][3] = c[3][2]; c[3][2] = c[3][1]; c[3][1] = c[3][0];
                           c[2][2] = c[2][1]; c[2][1] = c[2][0];
                                              c[1][1] = c[1][0];

        // Adders
        d[0] = (q[0]+   F) % MODULO;
        d[1] = (q[1]+d[0]) % MODULO;
        d[2] = (q[2]+d[1]) % MODULO;
        d[3] = (q[3]+d[2]) % MODULO;

        // Carries
        c[0][0] = (q[0]+   F) / MODULO;
        c[1][0] = (q[1]+d[0]) / MODULO;
        c[2][0] = (q[2]+d[1]) / MODULO;
        c[3][0] = (q[3]+d[2]) / MODULO;

        // Latch
        for (int i=0; i<4; i++) q[i] = d[i];

        return  c[0][0]
            +   c[1][0] -   c[1][1]
            +   c[2][0] - 2*c[2][1] +   c[2][2]
            +   c[3][0] - 3*c[3][1] + 3*c[3][2] - c[3][3]
            ;
    }

    bool operator != (MASH& m) // Debug
    {
        return memcmp(this, &m, sizeof(MASH)) != 0;
    }
};


void main()
{
    MASH m;

    for (int i=0; i<SEQLEN; i++) Delta[i]=m.Clock();

    if (m!=MASH()) return; // Debug
    for (i=0; i<SEQLEN; i++) if (Delta[i]!=m.Clock()) return; // Debug
    if (m!=MASH()) return; // Debug

    int t=0;

    for (i=0; i<SEQLEN; i++)
    {
        t += Delta[i];

        tVCO[i] = (i+1)*N + t;
        tREF[i] = (i+1)*N + (i+1)*double(F)/double(MODULO);

//      printf("%d  %d\n", i, Delta[i]);
    }
    
    // 2L = Sequence length (in VCO periods)
    const double L = SEQLEN/2 * (N + double(F)/double(MODULO));

    if (2*L != tVCO[SEQLEN-1]) return; // Debug

    for (double h=1; h<1e6; h=h<1000?h+1:h*1.001)
    {
        double a=0, b=0, n=floor(h);

        // Superposition of pulses
        for (i=0; i<SEQLEN; i++)
        {
            // Fourier series integrals over period 0....2L for nth harmonic
            a += 1.0/n/PI * (sin(n*PI/L*tVCO[i]) - sin(n*PI/L*(tREF[i])));
            b -= 1.0/n/PI * (cos(n*PI/L*tVCO[i]) - cos(n*PI/L*(tREF[i])));
        }

        printf("%g,%g\n", log10(h), 10*log10(a*a + b*b));
    }
}