/*************************************************************************\
  resample.c:
    Changes the size (and sampling rate) of an IFF 8SVX sample by an
    integer ratio (preferably 1:2), while avoiding aliasing effects by 
    calculating and applying a low pass filter first.
  arguments and their default values:
    <infile>                   ---
    <outfile>                  <infile>.out
    -l<filter length>          21
    -i<interpolate>            1
    -d<decimate>               2
    -g<gain>                   1.0
    -w<weight of stop band>    5.0
    -p<pass band>              0.7
    -v[erbose]
\*************************************************************************/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <signal.h>
#include <math.h>
#include "remez.h"

#define ID_FORM 0x464F524D
#define ID_8SVX 0x38535658
#define ID_VHDR 0x56484452
#define ID_BODY 0x424F4459

typedef unsigned long ULONG;
typedef unsigned short UWORD;
typedef unsigned char UBYTE;

struct {
  ULONG oneShotHiSamples,
        repeatHiSamples,
        samplesPerHiCycle;
  UWORD samplesPerSec;
  UBYTE ctOctave,
        sCompression;
  ULONG volume;  /* ULONG is not quite the truth, but close enough */
} vhdr;

void info();
void create_filter();
ULONG getlong();
long filter_body(long len);
void handle_iff();

FILE *infile, *outfile;
float fstop, stopwt = 5;  /* filter parameters */
float fpass = 0.7;
short a[NFMAX];  /* filter coefficients, fixed point representation */
int shift;       /* denotes fixed point unit */
int n = 0;
int intpl = 0, decmt = 0;
float gain = 1;
int verbose = 0;

main (int argc, char* argv[])
{
  ULONG l;
  int i,rawmode;
  long rawsize;
  char c, *s;
  char name1[80], name2[80];
  
  name1[0] = name2[0] = '\0';
  while (--argc) {
    s = *++argv;
    if (*s != '-') {
      if (name1[0] == '\0')
        strcpy(name1, s);
      else
        strcpy(name2, s);
    } else {
      switch (*++s) {
        case 'l':
          n = atoi(++s); break;
        case 'd': 
          decmt = atoi(++s); break;
        case 'i': 
          intpl = atoi(++s); break;
        case 'g':
          gain = atof(++s); break;
        case 'p':
          fpass = atof(++s); break;
        case 'w':
          stopwt = atof(++s); break;
        case 'v':
          verbose = 1; break;
        default:
          printf("Illegal parameter '%s', please use -ldiwgpv!\n", argv[0]);
      }
    }
  }
  /* add some default values: */
  if (intpl<1)  intpl = 1;
  if (decmt<1) {
    if (intpl == 1)  decmt = 2;
    else  decmt = 1;
  }
  if (n<1)
    n = (decmt > intpl) ? 10*decmt+1 : 10*intpl+1;
  /* make sure the arrays don't overflow */
  if (n>NFMAX)  n = NFMAX;
  if (name1[0] == '\0') {
    printf("Please specify an input file!\n");
    return 5;
  }
  if (name2[0] == '\0') {
    strcpy(name2, name1);
    strcat(name2, ".out");
  }
  if (!(infile = fopen(name1,"r"))) {
    printf("Can't open %s!\n", name1);
    return 10;
  }
  if (rawmode = (getlong() != ID_FORM)) {
    fseek(infile, 0L, SEEK_END); rawsize = ftell(infile);
  } else {
    l = getlong();  /* Dateigröße überlesen */
    if (getlong() != ID_8SVX) {
      printf("Not an IFF 8SVX file: %s!\n", name1);
      fclose(infile);
      return 5;
    }
  }
  if (!(outfile = fopen(name2,"w"))) {
    printf("Can't open %s for output!\n", name2);
    fclose(infile);
    return 10;
  }
  /* tell the user what his parameters will do */
  printf("Resampling \"%s\" as \"%s\",\n", name1, name2);
  printf("  using a %d-point FIR filter (-l), gain %.2f (-g),\n", n, gain);
  printf("  interpolation/decimation ratio is %d:%d (-i, -d),\n", intpl, decmt);
  printf("  pass band ends at %.2f below the stop band (-p),\n", fpass);
  printf("  relative weight of stop band errors is %.1f (-w).\n", stopwt);
  create_filter();
  rewind(infile);
  if (rawmode) {
    printf("Input file is not IFF, processing as raw 8bit-samples.\n");
    filter_body(rawsize);
  } else
    handle_iff();
  fclose(infile);
  fclose(outfile);

  return 0;
}

void create_filter()
/* Creates filter coefficients a[0..n-1], normalized, so that sum(a[i])==1,
 * and converted to a fixed point representation.
 */
{
  float sum=0, max=0;
  long isum=0, unit;
  int k;

  /* filter parameters: */
  fstop = (decmt > intpl) ? 0.5/decmt : 0.5/intpl;
  fpass *= fstop;   /* e. g. fpass = 0.7 * fstop */
  printf("Calculating filter parameters "); fflush(stdout);
  remez_dot = 1;    /* will create "....." output */
  if (makefilter(n) > 0) {   /* routine from "remez.o" */
    printf(" convergence failure, sorry.\n");
    if (n%2 == 0)
      printf("Better try '-l%d'.\n", n+1);
    exit(10);
  }
  /* convert the impulse response from float h[] to short a[], */
  /* but make sure it's normalized: */
  for (k=0; k<n; k++) {
    sum += h[k];
    if (fabs(h[k]) > max)
      max = fabs(h[k]);
  }
  /* find the maximum possible fixed point unit */
  shift = 0;
  do {
    unit = 1<<++shift;
  } while(unit * max*gain/sum < 0x8000);
  unit = 1<<--shift;
  for (k=0; k<n; k++) {
    a[k] = unit * h[k] * gain/sum;
    isum += a[k];
  }
  printf(" done, unit: 0x%x, checksum: 0x%x\n", unit, isum);
  if (verbose)  for (k=0; k<n; k++)
    printf("h[%2d] = %8.5f  (%6d)\n", k, h[k], a[k]);
  printf("Pass band [0:%.2f], ripple ± %.1f %%.\n", fpass,
    100*fabs(extreme(0)-1));
  printf("Stop band [%.2f:0.5], attenuation %.0f dB.\n", fstop,
    20*log10(fabs(extreme(0.5))));
}

float desire(float freq)
{
  if (freq > fstop)
    return 0;
  else
    return 1;
}

float weight(float freq)
{
  if (freq > fstop)
    return stopwt;
  else if (freq > fpass)
    return 0;
  else
    return 1;
}

int ctrl_c = 0;

void oops()
/* break-handler during filter operation */
{
  ctrl_c = 1;
}

long filter_body(long len)
/* Return value is the length of the filtered sample, may be more or less
 * than the original value, depending on what interpolation/decimation
 * ratio is applied.
 * Modulo operations are avoided, because they imply integer division, for
 * better runtime performance we use counters instead.
 */
{
  long max, loop, outlen = 0, warnings = 0;
  short b[NFMAX];  /* input signal buffer for the filter */
  int i, j, k, quiet = 0, n_half = n/2;
  int modi = 1, modd = 0, modth = len % 1000 + 1;   /* modulo counters */
  long sum, peak = 0;
  signed char c;  /* for the sample bytes */

  signal(SIGINT, *oops);  /* intercept user break (Ctrl-C) */
  for (i = 0; i<n; i++)
    b[i] = 0;
  i = 0;
  max = len*intpl + n_half;
  for (loop = 0; loop < max; loop++) {
    c = 0;
    if (--modi == 0) {
      modi = intpl;
      if (len-- > 0) {
        c = fgetc(infile);
        if (--modth == 0) {
          modth = 1000; 
          printf("\r%d input bytes to go. ",len+1); fflush(stdout);
          if (warnings>20 || (quiet && warnings>0))
            printf("%d warnings.\n", warnings);
          warnings = 0;
        }
      }
    }
    b[i] = c; for (j = 1; j<intpl; j++)  b[i] += c;    /* b[i] = intpl*c */
    i++; if (i == n)  i = 0;   /* points to the oldest value */
    if (loop == n_half)
      modd = 1;  /* enable output */
    if (--modd == 0) {
      modd = decmt;
      /* filter operation */
      sum = 0; k = i;
      for (j = 0; j<n; j++) {
        sum += a[j] * b[k];
        k++; if (k == n)  k = 0;
      }
      c = sum >>= shift;
      if (sum>127) c = 127;
      if (sum<-128) c = -128;
      if (c != sum) {
        if (!quiet) {
          printf("Warning: filter in saturation (%d)\n", sum);
          if (sum>peak)  peak = sum;
          if (warnings>50) {
            printf("Too many warnings. Maybe try '-g%.1f' or less.\n", 127*gain/peak);
            quiet = 1;
          }
        }
        warnings++;
      }
      fputc(c, outfile); outlen++;
    }
    if (ctrl_c) {
      printf("*** BREAK - output sample is incomplete!\n");
      break;
    }
  }
  printf("\r");
  signal(SIGINT, SIG_DFL);  /* restore normal handling */
  
  return outlen;
}

ULONG getlong()
/* returns 0 at EOF, not so brilliant, but sufficient for our purpose */
{
  ULONG l=0;
  int i, c;
  for (i=0; i<4; i++)
    l = (l<<8) | (c = fgetc(infile));
  if (c != EOF)
    return l;
  else
    return 0L;
}

void putlong(ULONG l)
{
  int i, byte[4];
  for (i=0; i<4; i++) {
    byte[3-i] = l & 0xff; l >>= 8;
  }
  for (i=0; i<4; i++)
    fputc(byte[i], outfile);
}

UWORD getword()
{
  UWORD w;
  w = fgetc(infile) << 8;
  return w | fgetc(infile);
}

void putword(UWORD w)
{
  fputc(w >> 8, outfile);
  fputc(w & 0xff, outfile);
}

void copy_chunk()
/* copy an IFF chunk to the output file (its ID has already been copied) */
{
  ULONG l;
  putlong(l = getlong());
  while (l--)
    fputc(fgetc(infile), outfile);
}

void handle_iff()
/* The input file has been recognized as valid IFF 8SVX, so copy it */
{
  ULONG l, ifflen, vhdlen, bodylen;
  fpos_t pos1, pos2, pos3;  /* markers for values calculated later */
  int i, done = 0;
  long delta;
  char s1[5];

  putlong(getlong());  /* "FORM" */
  fgetpos(outfile, &pos1);  /* pos1: for ifflen */
  putlong(ifflen = getlong());
  putlong(getlong());  /* "8SVX" */
  while (!done && (l = getlong())) {
    putlong(l);
    strncpy(s1, (char *)&l, 4); s1[4] = '\0'; printf("Chunk: %s\n",s1);
    if (l == ID_VHDR) {
      putlong(vhdlen = getlong());
      if (vhdlen < sizeof vhdr) {
        printf("Fatal error: incomplete VHDR chunk!\n");
        exit(10);  /* also closes the files! */
      }
      fgetpos(outfile, &pos2);  /* pos2: for the VHDR */
      vhdr.oneShotHiSamples = getlong();
      vhdr.repeatHiSamples = getlong();
      vhdr.samplesPerHiCycle = getlong();
      vhdr.samplesPerSec = getword();
      vhdr.ctOctave = fgetc(infile);
      vhdr.sCompression = fgetc(infile);
      vhdr.volume = getlong();
      for (i=0; i<vhdlen; i++) {
        fputc('\0', outfile);  /* will be filled in later */
        if (i >= sizeof vhdr)
          fgetc(infile);  /* shouldn't occur, but who knows */
      }
      if (vhdr.sCompression) {
        printf("Sorry, can't handle compressed samples ...\n");
        exit(10);
      }
    } else if (l == ID_BODY) {
      fgetpos(outfile, &pos3);  /* pos3: for bodylen */
      putlong(bodylen = getlong());
      delta = filter_body(bodylen) - bodylen;  /* -> i. a. (delta<0) */
      if ((bodylen + delta) & 1L) {  /* odd chunk size, fix it */
        fputc(0, outfile); delta++;
      }
      done = 1;   /* process no chunks after BODY */
    } else {
      copy_chunk();
    }
  }
  /* now adjust the sizes: */
  fsetpos(outfile, &pos1); putlong(ifflen + delta);
  fsetpos(outfile, &pos3); putlong(bodylen + delta);
  /* a bit more complicated for the VHDR: */
  fsetpos(outfile, &pos2);
  /* split the body in "one shot" and "repeat" part */
  l = intpl * vhdr.repeatHiSamples / decmt;
  putlong(bodylen + delta - l);  /* oneShotHiSamples */
  putlong(l);                    /* repeatHiSamples */
  putlong(intpl * vhdr.samplesPerHiCycle / decmt);
  l = intpl * vhdr.samplesPerSec / decmt;
  printf("Playback rate is now %d samples/sec\n", l);
  putword(l);
  fputc(vhdr.ctOctave, outfile);
  fputc(vhdr.sCompression, outfile);
  putlong(vhdr.volume);
}