/*************************************************************************\
  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.
\*************************************************************************/

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

char *version = "$VER: resample 1.1 (04.12.99)";


enum { FALSE, TRUE };

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

void help( char *s )
    {
    if( s )
        printf( "Illegal parameter '%s'!\n", s );
    printf( "Usage: resample <infile> [<outfile>] [options]\n" );
    printf( "Where <outfile> defaults to <infile>.out and options are:\n" );
    printf( "  -g<gain> (i.e. amplification)\n" );
    printf( "  -i<interpolation ratio>\n" );
    printf( "  -d<decimation ratio>\n" );
    printf( "  -l<filter length>\n" );
    printf( "  -p<pass band width>\n" );
    printf( "  -w<weight of stop band>\n" );
    printf( "  -v[erbose]\n" );
    exit( s ? 10 : 0 );
    }

int main( int argc, char* argv[] )
    {
    ULONG l;
    int rawmode = FALSE;
    long rawsize = 0;
    char *s;
    char name1[80], name2[80];

    if( argc < 2 )
        help( NULL );
    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:
                help( argv[ 0 ] );
            }
        }
    /* add some default values: */
    if( intpl < 1 )
        intpl = 1;
    if( decmt < 1 )
        decmt = (intpl == 1) ? 2 : 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( getlong() != ID_FORM )
        {
        fseek( infile, 0L, SEEK_END );
        rawmode = TRUE;
        rawsize = ftell( infile );
        }
    else
        {
        l = getlong();  /* skip file size */
        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%lx, checksum: 0x%lx\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;
    short b[ NFMAX ];           /* input signal buffer for the filter */
    int i, j, k, n_half = n/2;
    int modi = 1, modd = 0;     /* modulo counters */
    char modbyte = 0;
    time_t t0;
    long warnct = 0, warnbase = 0;
    long sum, peak = 0;
    signed char c;              /* for the sample bytes */

    signal(SIGINT, *oops);      /* intercept user break (Ctrl-C) */
    t0 = time( NULL );
    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 );
            }
        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;
                }
            sum >>= shift;
            c = sum;
            if( sum > 127 )
                c = 127;
            if( sum < -128 )
                c = -128;
            fputc( c, outfile );
            outlen++;
            /* overshoot statistics */
            warnbase++;
            if( c != sum )
                warnct++;
            if( peak < sum )
                peak = sum;
            if( peak < -sum )
                peak = -sum;
            }
        /* progress indicator, once every second */
        if( len == 0 || (modbyte++ == 0 && t0 != time( NULL )) )
            {
            float ratio, amount;

            t0 = time( NULL );
            printf( "\r%ld input samples to go, ", len );
            ratio = warnbase ? warnct * 100.0 / warnbase : 0;
            printf( "overshoot ratio %.1f %%", ratio );
            amount = peak * 100.0 / 128;
            printf( " (peak at %.0f %%)", amount );
            printf( "\e[K" );   /* that's ClrEoL */
            fflush( stdout );
            if( ratio > 2 )
                printf( "\n" );
            warnct = 0;
            warnbase = 0;
            peak = 0;
            }
        if( ctrl_c )
            {
            printf("\n*** BREAK - output sample is incomplete!");
            break;
            }
        }
    printf("\n");
    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));
    return ( c == EOF ) ? 0 : l;
    }

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 = 0;
    fpos_t pos1, pos2, pos3;  /* markers for values calculated later */
    int i, done = 0;
    long delta = 0;
    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);       /* implicitly 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;  /* -> usually (delta<0) */
            if( (bodylen + delta) & 1 )
                {               /* 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 %ld samples/sec\n", l );
    putword( l );
    fputc( vhdr.ctOctave, outfile );
    fputc( vhdr.sCompression, outfile );
    putlong( vhdr.volume );
    }