/*

  fract.c: Generator of fractional table io auodac.c

  This is a standard C program, not intended to run on VS1005.
  It calculates the values used in the fractional table needed for
  VS1005h accurate sample rate calculation.

 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

struct FractRate {
  int x, y;
  double val, fHighLimit;
  int highLimit;
  int bitPattern[8];
  int regVal;
} fractRate[256];

static int FractRateCompare(const struct FractRate *p1,
			    const struct FractRate *p2) {
  return (p1->val >= p2->val);
}

void FindBestFract(int inRate) {
  unsigned int f;
  int rate = inRate << 7;
  int over = 0;
  f = (unsigned int)((131072.0*8.0*256.0) * rate / 12288000 + 0.5);
  double backF1 = 0.0, backF2 = 0.0, backF3 = 0.0;
  unsigned int regVal = 0;
  static double totErr1 = 0.0, totErr2 = 0.0, totErr3 = 0.0;

  /* Rounding method */
  backF1 = ((f+128)&~0xFF) / 256.0 / 12288000.0 * 1125000.0;

  printf("%5d Hz, SRC %05x = %13.7f, ",
	 inRate, (f+128)>>8, backF1);

  backF2 = f / 256.0 / 12288000.0 * 1125000.0;

  {
    unsigned int ff = (unsigned int)f&0xFF;
    int i;

    if (ff >= 240) {
      /* Fractional part was so high that closest rounding is to next
	 integer value. */
      f += 256;
    } else {
      const struct FractRate *fr = fractRate;
      while (ff > fr->highLimit) {
	fr++;
      }
      regVal = fr->regVal;
    }
    for (i=0; i<=(regVal&7); i++) {
      if (regVal & (0x10<<i)) {
	over++;
      }
    }
  }
  
  backF3 = ((f>>8) + (regVal ? (double)over/((regVal&7)+1) : 0)) / 12288000.0 * 1125000.0;
  totErr1 += fabs(backF1-inRate);
  totErr2 += fabs(backF2-inRate);
  totErr3 += fabs(backF3-inRate);

  printf(": SRC %05x.%02x = %13.7f, regVal %04x, fract %13.7f, tE %7.5f %7.5f %7.5f\n",
    f>>8, f&0xFF, backF2, regVal, backF3,
    totErr1, totErr2, totErr3);
}

int main(int argc, char **argv) {
  int i, j;
  int n = 1;

  fractRate[0].val = 0.0;
  fractRate[0].y = 1;

  for (i=5; i<=8; i++) {
    int j;
    for (j=1; j<i; j++) {
      int k;
      fractRate[n].x = j;
      fractRate[n].y = i;
      fractRate[n].val = (double)j/i;
      fractRate[n].regVal = (i<<0)-1;
      for (k=0; k<j; k++) {
	int bit = k*i/j;
	fractRate[n].bitPattern[bit] = 1;
	fractRate[n].regVal |= (1<<4) << bit;
      }
      n++;
    }
  }
  fractRate[n].x = 8;
  fractRate[n].y = 8;
  fractRate[n].val = (double)8/8;
  n++;

  qsort(fractRate, n, sizeof(fractRate[0]), FractRateCompare);

  for (i=1; i<n; i++) {
    if (fractRate[i-1].val == fractRate[i].val) {
      printf("Purging first of:\n");
      printf("%d/%d = limit %8.6f, val %8.6f\n",
	     fractRate[i-1].x, fractRate[i-1].y, fractRate[i-1].fHighLimit, fractRate[i-1].val);
      printf("%d/%d = limit %8.6f, val %8.6f\n",
	     fractRate[i].x, fractRate[i].y, fractRate[i].fHighLimit, fractRate[i].val);
      memmove(fractRate+i-1, fractRate+i, (n-i)*sizeof(fractRate[0]));
      n--;
    }
  }
  printf("Purge complete\n\n");

  for (i=0; i<n; i++) {
    fractRate[i].fHighLimit = (fractRate[i+1].val + fractRate[i].val)*0.5;
    printf("  {%16.14f, 0x%04x}, /* %d/%d -> val %8.6f, dif %8.6f 0b",
	   fractRate[i].fHighLimit, fractRate[i].regVal,
	   fractRate[i].x, fractRate[i].y, fractRate[i].val,
	   (i) ? fractRate[i].val-fractRate[i-1].val : 0.0
	   );
    for (j=fractRate[i].y-2; j>=0; j--) {
      printf("%d", fractRate[i].bitPattern[j]);
    }
    printf(" */\n");
  }

  for (i=0; i<n-1; i++) {
    fractRate[i].fHighLimit = (fractRate[i+1].val + fractRate[i].val)*0.5;
    fractRate[i].highLimit = (int)(fractRate[i].fHighLimit*256+0.5);
    printf("  {%3d, 0x%04x}, /* %d/%d -> val %8.6f, limit %8.6f, 0b",
	   fractRate[i].highLimit, fractRate[i].regVal,
	   fractRate[i].x, fractRate[i].y, fractRate[i].val, fractRate[i].fHighLimit);
    for (j=fractRate[i].y-2; j>=0; j--) {
      printf("%d", fractRate[i].bitPattern[j]);
    }
    printf(" */\n");
  }

  printf("--- standard rates ---\n");
  FindBestFract(8000);
  FindBestFract(11025);
  FindBestFract(12000);
  FindBestFract(16000);
  FindBestFract(22050);
  FindBestFract(24000);
  FindBestFract(32000);
  FindBestFract(44100);
  FindBestFract(48000);
  printf("--- others ---\n");
  for (i=48000; i<48200; i++) {
    FindBestFract(i);
  }

  return EXIT_SUCCESS;
}
