home *** CD-ROM | disk | FTP | other *** search
/ Enigma Amiga Life 109 / EnigmaAmiga109CD.iso / software / musica / resample / src / resample.c < prev    next >
C/C++ Source or Header  |  1999-12-04  |  14KB  |  505 lines
  1. /*************************************************************************\
  2.   resample.c:
  3.     Changes the size (and sampling rate) of an IFF 8SVX sample by an
  4.     integer ratio (preferably 1:2), while avoiding aliasing effects by
  5.     calculating and applying a low pass filter first.
  6. \*************************************************************************/
  7.  
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include <stdlib.h>
  11. #include <signal.h>
  12. #include <time.h>
  13. #include <math.h>
  14. #include "remez.h"
  15.  
  16. char *version = "$VER: resample 1.1 (04.12.99)";
  17.  
  18.  
  19. enum { FALSE, TRUE };
  20.  
  21. #define ID_FORM 0x464F524D
  22. #define ID_8SVX 0x38535658
  23. #define ID_VHDR 0x56484452
  24. #define ID_BODY 0x424F4459
  25.  
  26. typedef unsigned long ULONG;
  27. typedef unsigned short UWORD;
  28. typedef unsigned char UBYTE;
  29.  
  30. struct
  31.     {
  32.     ULONG oneShotHiSamples,
  33.           repeatHiSamples,
  34.           samplesPerHiCycle;
  35.     UWORD samplesPerSec;
  36.     UBYTE ctOctave,
  37.           sCompression;
  38.     ULONG volume;  /* ULONG is not quite the truth, but close enough */
  39.     } vhdr;
  40.  
  41. void info();
  42. void create_filter();
  43. ULONG getlong();
  44. long filter_body(long len);
  45. void handle_iff();
  46.  
  47. FILE *infile, *outfile;
  48. float fstop, stopwt = 5;  /* filter parameters */
  49. float fpass = 0.7;
  50. short a[NFMAX];  /* filter coefficients, fixed point representation */
  51. int shift;       /* denotes fixed point unit */
  52. int n = 0;
  53. int intpl = 0, decmt = 0;
  54. float gain = 1;
  55. int verbose = 0;
  56.  
  57. void help( char *s )
  58.     {
  59.     if( s )
  60.         printf( "Illegal parameter '%s'!\n", s );
  61.     printf( "Usage: resample <infile> [<outfile>] [options]\n" );
  62.     printf( "Where <outfile> defaults to <infile>.out and options are:\n" );
  63.     printf( "  -g<gain> (i.e. amplification)\n" );
  64.     printf( "  -i<interpolation ratio>\n" );
  65.     printf( "  -d<decimation ratio>\n" );
  66.     printf( "  -l<filter length>\n" );
  67.     printf( "  -p<pass band width>\n" );
  68.     printf( "  -w<weight of stop band>\n" );
  69.     printf( "  -v[erbose]\n" );
  70.     exit( s ? 10 : 0 );
  71.     }
  72.  
  73. int main( int argc, char* argv[] )
  74.     {
  75.     ULONG l;
  76.     int rawmode = FALSE;
  77.     long rawsize = 0;
  78.     char *s;
  79.     char name1[80], name2[80];
  80.  
  81.     if( argc < 2 )
  82.         help( NULL );
  83.     name1[0] = name2[0] = '\0';
  84.     while( --argc )
  85.         {
  86.         s = *++argv;
  87.         if (*s != '-')
  88.             {
  89.             if( name1[0] == '\0' )
  90.                 strcpy(name1, s);
  91.             else
  92.                 strcpy(name2, s);
  93.             }
  94.         else switch (*++s)
  95.             {
  96.             case 'l':
  97.                 n = atoi(++s);
  98.                 break;
  99.             case 'd':
  100.                 decmt = atoi(++s);
  101.                 break;
  102.             case 'i':
  103.                 intpl = atoi(++s);
  104.                 break;
  105.             case 'g':
  106.                 gain = atof(++s);
  107.                 break;
  108.             case 'p':
  109.                 fpass = atof(++s);
  110.                 break;
  111.             case 'w':
  112.                 stopwt = atof(++s);
  113.                 break;
  114.             case 'v':
  115.                 verbose = 1;
  116.                 break;
  117.             default:
  118.                 help( argv[ 0 ] );
  119.             }
  120.         }
  121.     /* add some default values: */
  122.     if( intpl < 1 )
  123.         intpl = 1;
  124.     if( decmt < 1 )
  125.         decmt = (intpl == 1) ? 2 : 1;
  126.     if( n < 1 )
  127.         n = (decmt > intpl) ? 10*decmt+1 : 10*intpl+1;
  128.     /* make sure the arrays don't overflow */
  129.     if( n > NFMAX )
  130.         n = NFMAX;
  131.     if( name1[0] == '\0' )
  132.         {
  133.         printf("Please specify an input file!\n");
  134.         return 5;
  135.         }
  136.     if (name2[0] == '\0')
  137.         {
  138.         strcpy(name2, name1);
  139.         strcat(name2, ".out");
  140.         }
  141.     if (!(infile = fopen(name1,"r")))
  142.         {
  143.         printf("Can't open %s!\n", name1);
  144.         return 10;
  145.         }
  146.     if( getlong() != ID_FORM )
  147.         {
  148.         fseek( infile, 0L, SEEK_END );
  149.         rawmode = TRUE;
  150.         rawsize = ftell( infile );
  151.         }
  152.     else
  153.         {
  154.         l = getlong();  /* skip file size */
  155.         if (getlong() != ID_8SVX)
  156.             {
  157.             printf("Not an IFF 8SVX file: %s!\n", name1);
  158.             fclose(infile);
  159.             return 5;
  160.             }
  161.         }
  162.     if (!(outfile = fopen(name2,"w")))
  163.         {
  164.         printf("Can't open %s for output!\n", name2);
  165.         fclose(infile);
  166.         return 10;
  167.         }
  168.     /* tell the user what his parameters will do */
  169.     printf("Resampling \"%s\" as \"%s\",\n", name1, name2);
  170.     printf("  using a %d-point FIR filter (-l), gain %.2f (-g),\n", n, gain);
  171.     printf("  interpolation/decimation ratio is %d:%d (-i, -d),\n", intpl, decmt);
  172.     printf("  pass band ends at %.2f below the stop band (-p),\n", fpass);
  173.     printf("  relative weight of stop band errors is %.1f (-w).\n", stopwt);
  174.     create_filter();
  175.     rewind( infile );
  176.     if( rawmode )
  177.         {
  178.         printf( "Input file is not IFF, processing as raw 8bit-samples.\n" );
  179.         filter_body( rawsize );
  180.         }
  181.     else
  182.         handle_iff();
  183.     fclose( infile );
  184.     fclose( outfile );
  185.  
  186.     return 0;
  187.     }
  188.  
  189.  
  190.  
  191. void create_filter()
  192. /* Creates filter coefficients a[0..n-1], normalized, so that sum(a[i])==1,
  193.  * and converted to a fixed point representation.
  194.  */
  195.     {
  196.     float sum=0, max=0;
  197.     long isum=0, unit;
  198.     int k;
  199.  
  200.     /* filter parameters: */
  201.     fstop = (decmt > intpl) ? 0.5/decmt : 0.5/intpl;
  202.     fpass *= fstop;   /* e. g. fpass = 0.7 * fstop */
  203.     printf("Calculating filter parameters "); fflush(stdout);
  204.     remez_dot = 1;    /* will create "....." output */
  205.     if( makefilter( n ) > 0 )   /* routine from "remez.o" */
  206.         {
  207.         printf(" convergence failure, sorry.\n");
  208.         if( n%2 == 0 )
  209.             printf("Better try '-l%d'.\n", n+1);
  210.         exit( 10 );
  211.         }
  212.     /* convert the impulse response from float h[] to short a[], */
  213.     /* but make sure it's normalized: */
  214.     for( k=0; k<n; k++ )
  215.         {
  216.         sum += h[k];
  217.         if (fabs(h[k]) > max)
  218.             max = fabs(h[k]);
  219.         }
  220.     /* find the maximum possible fixed point unit */
  221.     shift = 0;
  222.     do  {
  223.       unit = 1 << ++shift;
  224.     } while( unit * max*gain/sum < 0x8000 );
  225.     unit = 1 << --shift;
  226.     for( k=0; k<n; k++ )
  227.         {
  228.         a[k] = unit * h[k] * gain/sum;
  229.         isum += a[k];
  230.         }
  231.     printf( " done, unit: 0x%lx, checksum: 0x%lx\n", unit, isum );
  232.     if( verbose )
  233.         for( k=0; k<n; k++ )
  234.             printf( "h[%2d] = %8.5f  (%6d)\n", k, h[k], a[k] );
  235.     printf( "Pass band [0:%.2f], ripple ± %.1f %%.\n", fpass,
  236.         100*fabs(extreme(0)-1));
  237.     printf( "Stop band [%.2f:0.5], attenuation %.0f dB.\n", fstop,
  238.         20*log10(fabs(extreme(0.5))) );
  239.     }
  240.  
  241.  
  242. float desire( float freq )
  243.     {
  244.     if( freq > fstop )
  245.         return 0;
  246.     else
  247.         return 1;
  248.     }
  249.  
  250. float weight(float freq)
  251.     {
  252.     if (freq > fstop)
  253.         return stopwt;
  254.     else if (freq > fpass)
  255.         return 0;
  256.     else
  257.         return 1;
  258.     }
  259.  
  260.  
  261.  
  262. int ctrl_c = 0;
  263.  
  264. void oops()
  265. /* break-handler during filter operation */
  266.     {
  267.     ctrl_c = 1;
  268.     }
  269.  
  270.  
  271.  
  272. long filter_body(long len)
  273. /* Return value is the length of the filtered sample, may be more or less
  274.  * than the original value, depending on what interpolation/decimation
  275.  * ratio is applied.
  276.  * Modulo operations are avoided, because they imply integer division, for
  277.  * better runtime performance we use counters instead.
  278.  */
  279.     {
  280.     long max, loop, outlen = 0;
  281.     short b[ NFMAX ];           /* input signal buffer for the filter */
  282.     int i, j, k, n_half = n/2;
  283.     int modi = 1, modd = 0;     /* modulo counters */
  284.     char modbyte = 0;
  285.     time_t t0;
  286.     long warnct = 0, warnbase = 0;
  287.     long sum, peak = 0;
  288.     signed char c;              /* for the sample bytes */
  289.  
  290.     signal(SIGINT, *oops);      /* intercept user break (Ctrl-C) */
  291.     t0 = time( NULL );
  292.     for (i = 0; i<n; i++)
  293.         b[i] = 0;
  294.     i = 0;
  295.     max = len*intpl + n_half;
  296.     for( loop = 0; loop < max; loop++ )
  297.         {
  298.         c = 0;
  299.         if( --modi == 0 )
  300.             {
  301.             modi = intpl;
  302.             if( len-- > 0 )
  303.                 c = fgetc( infile );
  304.             }
  305.         b[i] = c;
  306.         for( j = 1; j<intpl; j++ )
  307.             b[i] += c;          /* b[i] = intpl*c */
  308.         i++;
  309.         if( i == n )
  310.             i = 0;              /* points to the oldest value */
  311.         if( loop == n_half )
  312.           modd = 1;             /* enable output */
  313.         if( --modd == 0 )
  314.             {
  315.             modd = decmt;
  316.             /* filter operation */
  317.             sum = 0;
  318.             k = i;
  319.             for (j = 0; j<n; j++)
  320.                 {
  321.                 sum += a[j] * b[k];
  322.                 k++;
  323.                 if( k == n )
  324.                     k = 0;
  325.                 }
  326.             sum >>= shift;
  327.             c = sum;
  328.             if( sum > 127 )
  329.                 c = 127;
  330.             if( sum < -128 )
  331.                 c = -128;
  332.             fputc( c, outfile );
  333.             outlen++;
  334.             /* overshoot statistics */
  335.             warnbase++;
  336.             if( c != sum )
  337.                 warnct++;
  338.             if( peak < sum )
  339.                 peak = sum;
  340.             if( peak < -sum )
  341.                 peak = -sum;
  342.             }
  343.         /* progress indicator, once every second */
  344.         if( len == 0 || (modbyte++ == 0 && t0 != time( NULL )) )
  345.             {
  346.             float ratio, amount;
  347.  
  348.             t0 = time( NULL );
  349.             printf( "\r%ld input samples to go, ", len );
  350.             ratio = warnbase ? warnct * 100.0 / warnbase : 0;
  351.             printf( "overshoot ratio %.1f %%", ratio );
  352.             amount = peak * 100.0 / 128;
  353.             printf( " (peak at %.0f %%)", amount );
  354.             printf( "\e[K" );   /* that's ClrEoL */
  355.             fflush( stdout );
  356.             if( ratio > 2 )
  357.                 printf( "\n" );
  358.             warnct = 0;
  359.             warnbase = 0;
  360.             peak = 0;
  361.             }
  362.         if( ctrl_c )
  363.             {
  364.             printf("\n*** BREAK - output sample is incomplete!");
  365.             break;
  366.             }
  367.         }
  368.     printf("\n");
  369.     signal(SIGINT, SIG_DFL);  /* restore normal handling */
  370.  
  371.     return outlen;
  372.     }
  373.  
  374.  
  375. ULONG getlong()
  376. /* returns 0 at EOF, not so brilliant, but sufficient for our purpose */
  377.     {
  378.     ULONG l=0;
  379.     int i, c;
  380.   
  381.     for( i=0; i<4; i++ )
  382.         l = (l<<8) | (c = fgetc(infile));
  383.     return ( c == EOF ) ? 0 : l;
  384.     }
  385.  
  386. void putlong(ULONG l)
  387.     {
  388.     int i, byte[4];
  389.  
  390.     for( i=0; i<4; i++ )
  391.         {
  392.         byte[3-i] = l & 0xff;
  393.         l >>= 8;
  394.         }
  395.     for( i=0; i<4; i++ )
  396.         fputc(byte[i], outfile);
  397.     }
  398.  
  399. UWORD getword()
  400.     {
  401.     UWORD w;
  402.     w = fgetc( infile ) << 8;
  403.     return w | fgetc(infile);
  404.     }
  405.  
  406. void putword( UWORD w )
  407.     {
  408.     fputc( w >> 8, outfile );
  409.     fputc( w & 0xff, outfile );
  410.     }
  411.  
  412. void copy_chunk()
  413. /* copy an IFF chunk to the output file (its ID has already been copied) */
  414.     {
  415.     ULONG l;
  416.     putlong( l = getlong() );
  417.     while (l--)
  418.         fputc( fgetc(infile), outfile );
  419.     }
  420.  
  421. void handle_iff()
  422. /* The input file has been recognized as valid IFF 8SVX, so copy it */
  423.     {
  424.     ULONG l, ifflen, vhdlen, bodylen = 0;
  425.     fpos_t pos1, pos2, pos3;  /* markers for values calculated later */
  426.     int i, done = 0;
  427.     long delta = 0;
  428.     char s1[5];
  429.  
  430.     putlong( getlong() );       /* "FORM" */
  431.     fgetpos( outfile, &pos1 );  /* pos1: for ifflen */
  432.     putlong( ifflen = getlong() );
  433.     putlong( getlong() );       /* "8SVX" */
  434.     while( !done && (l = getlong()) )
  435.         {
  436.         putlong( l );
  437.         strncpy( s1, (char *)&l, 4 );
  438.         s1[ 4 ] = '\0';
  439.         printf( "Chunk: %s\n", s1 );
  440.         if (l == ID_VHDR)
  441.             {
  442.             putlong( vhdlen = getlong() );
  443.             if( vhdlen < sizeof( vhdr ) )
  444.                 {
  445.                 printf("Fatal error: incomplete VHDR chunk!\n");
  446.                 exit(10);       /* implicitly closes the files! */
  447.                 }
  448.             fgetpos( outfile, &pos2 );  /* pos2: for the VHDR */
  449.             vhdr.oneShotHiSamples = getlong();
  450.             vhdr.repeatHiSamples = getlong();
  451.             vhdr.samplesPerHiCycle = getlong();
  452.             vhdr.samplesPerSec = getword();
  453.             vhdr.ctOctave = fgetc( infile );
  454.             vhdr.sCompression = fgetc( infile );
  455.             vhdr.volume = getlong();
  456.             for( i=0; i<vhdlen; i++ )
  457.                 {
  458.                 fputc('\0', outfile);   /* will be filled in later */
  459.                 if( i >= sizeof( vhdr ) )
  460.                     fgetc( infile );    /* shouldn't occur, but who knows */
  461.                 }
  462.             if( vhdr.sCompression )
  463.                 {
  464.                 printf("Sorry, can't handle compressed samples...\n");
  465.                 exit(10);
  466.                 }
  467.             }
  468.         else if (l == ID_BODY)
  469.             {
  470.             fgetpos( outfile, &pos3 );  /* pos3: for bodylen */
  471.             putlong( bodylen = getlong() );
  472.             delta = filter_body(bodylen) - bodylen;  /* -> usually (delta<0) */
  473.             if( (bodylen + delta) & 1 )
  474.                 {               /* odd chunk size, fix it */
  475.                 fputc( 0, outfile );
  476.                 delta++;
  477.                 }
  478.             done = 1;           /* process no chunks after BODY */
  479.             }
  480.         else
  481.             copy_chunk();
  482.         }
  483.     /* now adjust the sizes: */
  484.     fsetpos( outfile, &pos1 );
  485.     putlong( ifflen + delta );
  486.     fsetpos( outfile, &pos3 );
  487.     putlong( bodylen + delta );
  488.  
  489.     /* a bit more complicated for the VHDR: */
  490.     fsetpos( outfile, &pos2 );
  491.     /* split the body in "one shot" and "repeat" part */
  492.     l = intpl * vhdr.repeatHiSamples / decmt;
  493.     putlong( bodylen + delta - l );     /* oneShotHiSamples */
  494.     putlong( l );                       /* repeatHiSamples */
  495.     putlong( intpl * vhdr.samplesPerHiCycle / decmt );
  496.     l = intpl * vhdr.samplesPerSec / decmt;
  497.     printf( "Playback rate is now %ld samples/sec\n", l );
  498.     putword( l );
  499.     fputc( vhdr.ctOctave, outfile );
  500.     fputc( vhdr.sCompression, outfile );
  501.     putlong( vhdr.volume );
  502.     }
  503.  
  504.  
  505.