PDL-Audio/xlib.c

#include "xlib.h"

/*
 * the following functions were originally taken from sox-12.16/libst.c
 * license is unclear, but the header file contained this notice:
 */

/*
** Copyright (C) 1989 by Jef Poskanzer.
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/

/*
** This routine converts from linear to ulaw.
**
** Craig Reese: IDA/Supercomputing Research Center
** Joe Campbell: Department of Defense
** 29 September 1989
**
** References:
** 1) CCITT Recommendation G.711  (very difficult to follow)
** 2) "A New Digital Technique for Implementation of Any
**     Continuous PCM Companding Law," Villeret, Michel,
**     et al. 1973 IEEE Int. Conf. on Communications, Vol 1,
**     1973, pg. 11.12-11.17
** 3) MIL-STD-188-113,"Interoperability and Performance Standards
**     for Analog-to_Digital Conversion Techniques,"
**     17 February 1987
**
** Input: Signed 16 bit linear sample
** Output: 8 bit ulaw sample
*/

#undef ZEROTRAP      /* turn off the trap as per the MIL-STD */
#define uBIAS 0x84   /* define the add-in bias for 16 bit samples */
#define uCLIP 32635
#define ACLIP 31744

unsigned char
st_linear_to_ulaw(int sample)
    {
    static int exp_lut[256] = {0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,
                               4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
                               5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
                               5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
                               6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
                               6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
                               6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
                               6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
                               7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7};
    int sign, exponent, mantissa;
    unsigned char ulawbyte;

    /* Get the sample into sign-magnitude. */
    sign = (sample >> 8) & 0x80;                /* set aside the sign */
    if ( sign != 0 ) sample = -sample;          /* get magnitude */
    if ( sample > uCLIP ) sample = uCLIP;               /* clip the magnitude */

    /* Convert from 16 bit linear to ulaw. */
    sample = sample + uBIAS;
    exponent = exp_lut[( sample >> 7 ) & 0xFF];
    mantissa = ( sample >> ( exponent + 3 ) ) & 0x0F;
    ulawbyte = ~ ( sign | ( exponent << 4 ) | mantissa );
#ifdef ZEROTRAP
    if ( ulawbyte == 0 ) ulawbyte = 0x02;       /* optional CCITT trap */
#endif

    return ulawbyte;
    }

/*
** This routine converts from ulaw to 16 bit linear.
**
** Craig Reese: IDA/Supercomputing Research Center
** 29 September 1989
**
** References:
** 1) CCITT Recommendation G.711  (very difficult to follow)
** 2) MIL-STD-188-113,"Interoperability and Performance Standards
**     for Analog-to_Digital Conversion Techniques,"
**     17 February 1987
**
** Input: 8 bit ulaw sample
** Output: signed 16 bit linear sample
*/

int
st_ulaw_to_linear(unsigned char ulawbyte)
    {
    static int exp_lut[8] = { 0, 132, 396, 924, 1980, 4092, 8316, 16764 };
    int sign, exponent, mantissa, sample;

    ulawbyte = ~ ulawbyte;
    sign = ( ulawbyte & 0x80 );
    exponent = ( ulawbyte >> 4 ) & 0x07;
    mantissa = ulawbyte & 0x0F;
    sample = exp_lut[exponent] + ( mantissa << ( exponent + 3 ) );
    if ( sign != 0 ) sample = -sample;

    return sample;
    }

/*
 * A-law routines by Graeme W. Gill.
 * Date: 93/5/7
 *
 * References:
 * 1) CCITT Recommendation G.711
 *
 * These routines were used to create the fast
 * lookup tables.
 */

#define ACLIP 31744

unsigned char
st_linear_to_Alaw(int sample)
    {
    static int exp_lut[128] = {1,1,2,2,3,3,3,3,
                               4,4,4,4,4,4,4,4,
                               5,5,5,5,5,5,5,5,
                               5,5,5,5,5,5,5,5,
                               6,6,6,6,6,6,6,6,
                               6,6,6,6,6,6,6,6,
                               6,6,6,6,6,6,6,6,
                               6,6,6,6,6,6,6,6,
                               7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7,
                               7,7,7,7,7,7,7,7};

    int sign, exponent, mantissa;
    unsigned char Alawbyte;

    /* Get the sample into sign-magnitude. */
    sign = ((~sample) >> 8) & 0x80;             /* set aside the sign */
    if ( sign == 0 ) sample = -sample;          /* get magnitude */
    if ( sample > ACLIP ) sample = ACLIP;       /* clip the magnitude */

    /* Convert from 16 bit linear to ulaw. */
    if (sample >= 256)
        {
        exponent = exp_lut[( sample >> 8 ) & 0x7F];
        mantissa = ( sample >> ( exponent + 3 ) ) & 0x0F;
        Alawbyte = (( exponent << 4 ) | mantissa);
        }
    else
        Alawbyte = (sample >> 4);
    Alawbyte ^= (sign ^ 0x55);

    return Alawbyte;
    }

int
st_Alaw_to_linear(unsigned char Alawbyte)
    {
    static int exp_lut[8] = { 0, 264, 528, 1056, 2112, 4224, 8448, 16896 };
    int sign, exponent, mantissa, sample;

    Alawbyte ^= 0x55;
    sign = ( Alawbyte & 0x80 );
    Alawbyte &= 0x7f;                   /* get magnitude */
    if (Alawbyte >= 16)
        {
        exponent = (Alawbyte >> 4 ) & 0x07;
        mantissa = Alawbyte & 0x0F;
        sample = exp_lut[exponent] + ( mantissa << ( exponent + 3 ) );
        }
    else
        sample = (Alawbyte << 4) + 8;
    if ( sign == 0 ) sample = -sample;

    return sample;
    }

/* adapted from clm.c, by Bill Schottstaedt C<bil@ccrma.stanford.edu> */

#define SRC_SINC_DENSITY 1000
#define SRC_SINC_WIDTH 10

static Float **sinc_tables = 0;
static int *sinc_widths = 0;
static int sincs = 0;

static Float *
init_sinc_table (int width)
{
  int i, size, loc;
  Float sinc_freq, win_freq, sinc_phase, win_phase;
  for (i = 0; i < sincs; i++)
    if (sinc_widths[i] == width)
      return (sinc_tables[i]);

  if (sincs == 0)
    {
      sinc_tables = (Float **) calloc (8, sizeof (Float *));
      sinc_widths = (int *) calloc (8, sizeof (int));
      sincs = 8;
      loc = 0;
    }
  else
    {
      loc = -1;
      for (i = 0; i < sincs; i++)
        if (sinc_widths[i] == 0)
          {
            loc = i;
            break;
          }

      if (loc == -1)
        {
          sinc_tables = (Float **) realloc (sinc_tables, (sincs + 8) * sizeof (Float *));
          sinc_widths = (int *) realloc (sinc_widths, (sincs + 8) * sizeof (int));
          for (i = sincs; i < (sincs + 8); i++)
            {
              sinc_widths[i] = 0;
              sinc_tables[i] = NULL;
            }

          loc = sincs;
          sincs += 8;
        }
    }
  sinc_tables[loc] = (Float *) calloc (width * SRC_SINC_DENSITY + 1, sizeof (Float));
  sinc_widths[loc] = width;
  size = width * SRC_SINC_DENSITY;
  sinc_freq = M_PI / (Float) SRC_SINC_DENSITY;
  win_freq = M_PI / (Float) size;
  sinc_tables[loc][0] = 1.0;

  for (i = 1, sinc_phase = sinc_freq, win_phase = win_freq; i < size; i++, sinc_phase += sinc_freq, win_phase += win_freq)
    sinc_tables[loc][i] = sin (sinc_phase) * (0.5 + 0.5 * cos (win_phase)) / sinc_phase;

  return (sinc_tables[loc]);
}

void
mus_src (Float *input, int inpsize, Float *output, int outsize, Float srate, Float *sr_mod, int width)
{
  int i, lim, len, fsx, k, loc;

  Float x, xx, factor, *data, *sinc_table, sum, zf, srx, incr;

  Float *in0 = input;
  Float *in1 = input + inpsize;
  Float *out1 = output + outsize;

  if (width == 0)
    width = SRC_SINC_WIDTH;

  x = 0.0;
  lim = 2 * width;
  len = width * SRC_SINC_DENSITY;
  data = (Float *) calloc (lim + 1, sizeof (Float));
  sinc_table = init_sinc_table (width);

  for (i = width; i < lim; i++) data[i] = *input++;

  while (output < out1)
    {
      fsx = (int)x;
      if (fsx > 0)
        {
          /* realign data, reset x */
          for (i = fsx, loc = 0; i < lim; i++, loc++)
            data[loc] = data[i];

          for (i = loc; i < lim; i++)
            {
              if (srx < 0)
                input = (input > in0 ? input : in1) - 1;
              else
                input = input < in1 ? input+1 : in0;

              data[i] = *input;
            }

          x -= fsx;
        }

      srx = srate + (sr_mod ? *sr_mod++ : 0);
      srx = srx ? fabs (srx) : 0.001;
      factor = srx > 1 ? 1 / srx : 1;

      sum = 0.0;
      zf = factor * SRC_SINC_DENSITY;
      xx = zf * (1.0 - x - width);
      for (i = 0; i < lim; i++)
        {
          /* we're moving backwards in the data array, so xx has to mimic that (hence the '1.0 - x') */
          k = abs ((int)xx);

          if (k < len)
            sum += data[i] * sinc_table[k];

          xx += zf;
        }

      x += srx;
      *output++ = sum * factor;
    }

  free (data);
}

static unsigned long randx = 1;

static int 
irandom (int amp)
{
  int val;

  randx = randx * 1103515245 + 12345;
  val = (unsigned int) (randx >> 16) & 0x7fff;
  return ((int) (amp * (((Float) val / 32768))));
}

#define max(a,b) ((a)>(b) ? (a) : (b))
#define min(a,b) ((a)<(b) ? (a) : (b))

void
mus_granulate (Float *input, int insize,
               Float *output, int outsize,
               Float expansion, Float flength, Float scaler,
               Float hop, Float ramp, Float jitter, int max_size)
/* hop, jitter, length (*= smapling_rate) */
{
  int length = (int)ceil (flength);
  int rmp = (int) (ramp * length);
  int output_hop = (int)hop;
  int input_hop = (int)(output_hop / expansion);
  int s20 = (int) (jitter / 20);
  int s50 = (int) (jitter / 50);
  int outlen = max_size > 0 ? min ((int)(hop + flength), max_size)  : (int)(hop + flength);
  int in_data_len = outlen + s20 + 1;
  int in_data_start = in_data_len;
  Float *data = (Float *) calloc (outlen, sizeof (Float));
  Float *in_data = (Float *) calloc (in_data_len, sizeof (Float));

  Float *in1  = input  + insize;
  Float *out1 = output + outsize;

  int ctr = 0;
  Float cur_out = 0;

  int start, len, end, i, j, k;
  int  steady_end, curstart;
  Float incr, result, amp;

  if (s50 > output_hop)
    s50 = output_hop;

  for(;;)
    {
      while (ctr < cur_out)
        {
          *output++ = data[ctr++];
          if (output >= out1)
            goto ok;
        }

      start = cur_out;
      end = length - start;
      
      if (end <= 0)
        end = 0;
      else
        for (i = 0, j = start; i < end; i++, j++)
          data[i] = data[j];

      for (i = end; i < outlen; i++)
        data[i] = 0;

      start = in_data_start;
      len = in_data_len;

      if (start > len)
        {
          input += start - len;
          input = input < in1 ? input : in1;
          start = len;
        }
      else if (start < len)
        for (i = 0, k = start; k < len; i++, k++)
          in_data[i] = in_data[k];

      for (i = len - start; i < len; i++)
        {
          in_data[i] = *input;
          input = input < in1 ? input+1 : input;
        }

      in_data_start = input_hop;

      amp = 0.0;
      incr = scaler / (Float) rmp;
      steady_end = length - rmp;
      curstart = irandom (s20);

      for (i = 0, j = curstart; i < length; i++, j++)
        {
          data[i] += (amp * in_data[j]);

          if (i < rmp)
            amp += incr;
          else if (i > steady_end)
            amp -= incr;
        }

      ctr -= cur_out;
      cur_out = output_hop + irandom (s50) - (s50 >> 1);
    }

  ok:
  free (data);
  free (in_data);
}

static void
mus_shuffle (Float* rl, Float* im, int n)
{
  /* bit reversal */

  int i,m,j;
  Float tempr,tempi;
  j=0;
  for (i=0;i<n;i++)
    {
      if (j>i)
        {
          tempr = rl[j];
          tempi = im[j];
          rl[j] = rl[i];
          im[j] = im[i];
          rl[i] = tempr;
          im[i] = tempi;
        }
      m = n>>1;
      while ((m>=2) && (j>=m))
        {
          j -= m;
          m = m>>1;
        }
      j += m;
    }
}

static void
mus_fft (Float *rl, Float *im, int n, int isign)
{
  /* standard fft: real part in rl, imaginary in im,        */
  /* rl and im are zero-based.                              */
  int mmax,j,pow,prev,lg,i,ii,jj,ipow;
  Float wrs,wis,tempr,tempi;
  double wr,wi,theta,wtemp,wpr,wpi;
  ipow = (int)(ceil(log(n)/log(2.0)));
  mus_shuffle(rl,im,n);
  mmax = 2;
  prev = 1;
  pow = (int)(n*0.5);
  theta = (M_PI*isign);
  for (lg=0;lg<ipow;lg++)
    {
      wpr = cos(theta);
      wpi = sin(theta);
      wr = 1.0;
      wi = 0.0;
      for (ii=0;ii<prev;ii++)
        {
          wrs = (Float) wr;
          wis = (Float) wi;
#ifdef LINUX 
          if (isnan(wis)) wis=0.0;
#endif
          i = ii;
          j = ii + prev;
          for (jj=0;jj<pow;jj++)
            {
              tempr = wrs*rl[j] - wis*im[j];
              tempi = wrs*im[j] + wis*rl[j];
              rl[j] = rl[i] - tempr;
              im[j] = im[i] - tempi;
              rl[i] += tempr;
              im[i] += tempi;
              i += mmax;
              j += mmax;
            }
          wtemp = wr;
          wr = (wr*wpr) - (wi*wpi);
          wi = (wi*wpr) + (wtemp*wpi);
        }
      pow = (int)(pow*0.5);
      prev = mmax;
      theta = theta*0.5;
      mmax = mmax*2;
    }
}

static void
mus_convolution (Float* rl1, Float* rl2, int n)
{
  /* convolves two real arrays.                                           */
  /* rl1 and rl2 are assumed to be set up correctly for the convolution   */
  /* (that is, rl1 (the "signal") is zero-padded by length of             */
  /* (non-zero part of) rl2 and rl2 is stored in wrap-around order)       */
  /* We treat rl2 as the imaginary part of the first fft, then do         */
  /* the split, scaling, and (complex) spectral multiply in one step.     */
  /* result in rl1                                                        */

  int j,n2,nn2;
  Float rem,rep,aim,aip,invn;

  mus_fft(rl1,rl2,n,1);
  
  n2=(int)(n*0.5);
  invn = 0.25/n;
  rl1[0] = ((rl1[0]*rl2[0])/n);
  rl2[0] = 0.0;

  for (j=1;j<=n2;j++)
    {
      nn2 = n-j;
      rep = (rl1[j]+rl1[nn2]);
      rem = (rl1[j]-rl1[nn2]);
      aip = (rl2[j]+rl2[nn2]);
      aim = (rl2[j]-rl2[nn2]);

      rl1[j] = invn*(rep*aip + aim*rem);
      rl1[nn2] = rl1[j];
      rl2[j] = invn*(aim*aip - rep*rem);
      rl2[nn2] = -rl2[j];
    }
  
  mus_fft(rl1,rl2,n,-1);
}

void
mus_convolve (Float * input, Float * output, int size, Float * filter, int fftsize, int filtersize)
{
  int fftsize2 = fftsize >> 1;
  Float *rl1, *rl2, *buf;
  Float *in1 = input + size;
  int ctr = fftsize2;
  int i, j;

  rl1 = (Float *) calloc (fftsize, sizeof (Float));
  rl2 = (Float *) calloc (fftsize, sizeof (Float));
  buf = (Float *) calloc (fftsize, sizeof (Float));

  while (size > 0)
    {
      ctr++;
      if (ctr >= fftsize2)
        {
          for (i = 0, j = fftsize2; i < fftsize2; i++, j++)
            {
              buf[i] = buf[j];
              buf[j] = 0.0;
              rl1[i] = *input;
              rl1[j] = 0.0;
              rl2[i] = 0.0;
              rl2[j] = 0.0;

              input = input < in1 ? input+1 : input;
            }

          for (i = 0; i < filtersize; i++)
            rl2[i] = filter[i];

          mus_convolution (rl1, rl2, fftsize);

          for (i = 0, j = fftsize2; i < fftsize2; i++, j++)
            {
              buf[i] += rl1[i];
              buf[j] = rl1[j];
            }

          ctr = 0;
        }

      *output++ = buf[ctr];
      size--;
    }
}
Revision:	1.1
Committed:	Tue Dec 28 01:05:16 2004 UTC (19 years, 4 months ago) by root
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rel-1_1, rel-1_2, HEAD
Log Message:	* empty log message *
#	User	Rev	Content
1	root	1.1	#include "xlib.h"
2
3			/*
4			* the following functions were originally taken from sox-12.16/libst.c
5			* license is unclear, but the header file contained this notice:
6			*/
7
8			/*
9			** Copyright (C) 1989 by Jef Poskanzer.
10			**
11			** Permission to use, copy, modify, and distribute this software and its
12			** documentation for any purpose and without fee is hereby granted, provided
13			** that the above copyright notice appear in all copies and that both that
14			** copyright notice and this permission notice appear in supporting
15			** documentation. This software is provided "as is" without express or
16			** implied warranty.
17			*/
18
19			/*
20			** This routine converts from linear to ulaw.
21			**
22			** Craig Reese: IDA/Supercomputing Research Center
23			** Joe Campbell: Department of Defense
24			** 29 September 1989
25			**
26			** References:
27			** 1) CCITT Recommendation G.711 (very difficult to follow)
28			** 2) "A New Digital Technique for Implementation of Any
29			** Continuous PCM Companding Law," Villeret, Michel,
30			** et al. 1973 IEEE Int. Conf. on Communications, Vol 1,
31			** 1973, pg. 11.12-11.17
32			** 3) MIL-STD-188-113,"Interoperability and Performance Standards
33			** for Analog-to_Digital Conversion Techniques,"
34			** 17 February 1987
35			**
36			** Input: Signed 16 bit linear sample
37			** Output: 8 bit ulaw sample
38			*/
39
40			#undef ZEROTRAP /* turn off the trap as per the MIL-STD */
41			#define uBIAS 0x84 /* define the add-in bias for 16 bit samples */
42			#define uCLIP 32635
43			#define ACLIP 31744
44
45			unsigned char
46			st_linear_to_ulaw(int sample)
47			{
48			static int exp_lut[256] = {0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,
49			4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
50			5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
51			5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
52			6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
53			6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
54			6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
55			6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
56			7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
57			7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
58			7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
59			7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
60			7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
61			7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
62			7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
63			7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7};
64			int sign, exponent, mantissa;
65			unsigned char ulawbyte;
66
67			/* Get the sample into sign-magnitude. */
68			sign = (sample >> 8) & 0x80; /* set aside the sign */
69			if ( sign != 0 ) sample = -sample; /* get magnitude */
70			if ( sample > uCLIP ) sample = uCLIP; /* clip the magnitude */
71
72			/* Convert from 16 bit linear to ulaw. */
73			sample = sample + uBIAS;
74			exponent = exp_lut[( sample >> 7 ) & 0xFF];
75			mantissa = ( sample >> ( exponent + 3 ) ) & 0x0F;
76			ulawbyte = ~ ( sign \| ( exponent << 4 ) \| mantissa );
77			#ifdef ZEROTRAP
78			if ( ulawbyte == 0 ) ulawbyte = 0x02; /* optional CCITT trap */
79			#endif
80
81			return ulawbyte;
82			}
83
84			/*
85			** This routine converts from ulaw to 16 bit linear.
86			**
87			** Craig Reese: IDA/Supercomputing Research Center
88			** 29 September 1989
89			**
90			** References:
91			** 1) CCITT Recommendation G.711 (very difficult to follow)
92			** 2) MIL-STD-188-113,"Interoperability and Performance Standards
93			** for Analog-to_Digital Conversion Techniques,"
94			** 17 February 1987
95			**
96			** Input: 8 bit ulaw sample
97			** Output: signed 16 bit linear sample
98			*/
99
100			int
101			st_ulaw_to_linear(unsigned char ulawbyte)
102			{
103			static int exp_lut[8] = { 0, 132, 396, 924, 1980, 4092, 8316, 16764 };
104			int sign, exponent, mantissa, sample;
105
106			ulawbyte = ~ ulawbyte;
107			sign = ( ulawbyte & 0x80 );
108			exponent = ( ulawbyte >> 4 ) & 0x07;
109			mantissa = ulawbyte & 0x0F;
110			sample = exp_lut[exponent] + ( mantissa << ( exponent + 3 ) );
111			if ( sign != 0 ) sample = -sample;
112
113			return sample;
114			}
115
116			/*
117			* A-law routines by Graeme W. Gill.
118			* Date: 93/5/7
119			*
120			* References:
121			* 1) CCITT Recommendation G.711
122			*
123			* These routines were used to create the fast
124			* lookup tables.
125			*/
126
127			#define ACLIP 31744
128
129			unsigned char
130			st_linear_to_Alaw(int sample)
131			{
132			static int exp_lut[128] = {1,1,2,2,3,3,3,3,
133			4,4,4,4,4,4,4,4,
134			5,5,5,5,5,5,5,5,
135			5,5,5,5,5,5,5,5,
136			6,6,6,6,6,6,6,6,
137			6,6,6,6,6,6,6,6,
138			6,6,6,6,6,6,6,6,
139			6,6,6,6,6,6,6,6,
140			7,7,7,7,7,7,7,7,
141			7,7,7,7,7,7,7,7,
142			7,7,7,7,7,7,7,7,
143			7,7,7,7,7,7,7,7,
144			7,7,7,7,7,7,7,7,
145			7,7,7,7,7,7,7,7,
146			7,7,7,7,7,7,7,7,
147			7,7,7,7,7,7,7,7};
148
149			int sign, exponent, mantissa;
150			unsigned char Alawbyte;
151
152			/* Get the sample into sign-magnitude. */
153			sign = ((~sample) >> 8) & 0x80; /* set aside the sign */
154			if ( sign == 0 ) sample = -sample; /* get magnitude */
155			if ( sample > ACLIP ) sample = ACLIP; /* clip the magnitude */
156
157			/* Convert from 16 bit linear to ulaw. */
158			if (sample >= 256)
159			{
160			exponent = exp_lut[( sample >> 8 ) & 0x7F];
161			mantissa = ( sample >> ( exponent + 3 ) ) & 0x0F;
162			Alawbyte = (( exponent << 4 ) \| mantissa);
163			}
164			else
165			Alawbyte = (sample >> 4);
166			Alawbyte ^= (sign ^ 0x55);
167
168			return Alawbyte;
169			}
170
171			int
172			st_Alaw_to_linear(unsigned char Alawbyte)
173			{
174			static int exp_lut[8] = { 0, 264, 528, 1056, 2112, 4224, 8448, 16896 };
175			int sign, exponent, mantissa, sample;
176
177			Alawbyte ^= 0x55;
178			sign = ( Alawbyte & 0x80 );
179			Alawbyte &= 0x7f; /* get magnitude */
180			if (Alawbyte >= 16)
181			{
182			exponent = (Alawbyte >> 4 ) & 0x07;
183			mantissa = Alawbyte & 0x0F;
184			sample = exp_lut[exponent] + ( mantissa << ( exponent + 3 ) );
185			}
186			else
187			sample = (Alawbyte << 4) + 8;
188			if ( sign == 0 ) sample = -sample;
189
190			return sample;
191			}
192
193			/* adapted from clm.c, by Bill Schottstaedt C<bil@ccrma.stanford.edu> */
194
195			#define SRC_SINC_DENSITY 1000
196			#define SRC_SINC_WIDTH 10
197
198			static Float **sinc_tables = 0;
199			static int *sinc_widths = 0;
200			static int sincs = 0;
201
202			static Float *
203			init_sinc_table (int width)
204			{
205			int i, size, loc;
206			Float sinc_freq, win_freq, sinc_phase, win_phase;
207			for (i = 0; i < sincs; i++)
208			if (sinc_widths[i] == width)
209			return (sinc_tables[i]);
210
211			if (sincs == 0)
212			{
213			sinc_tables = (Float *) calloc (8, sizeof (Float ));
214			sinc_widths = (int *) calloc (8, sizeof (int));
215			sincs = 8;
216			loc = 0;
217			}
218			else
219			{
220			loc = -1;
221			for (i = 0; i < sincs; i++)
222			if (sinc_widths[i] == 0)
223			{
224			loc = i;
225			break;
226			}
227
228			if (loc == -1)
229			{
230			sinc_tables = (Float *) realloc (sinc_tables, (sincs + 8) sizeof (Float *));
231			sinc_widths = (int ) realloc (sinc_widths, (sincs + 8) sizeof (int));
232			for (i = sincs; i < (sincs + 8); i++)
233			{
234			sinc_widths[i] = 0;
235			sinc_tables[i] = NULL;
236			}
237
238			loc = sincs;
239			sincs += 8;
240			}
241			}
242			sinc_tables[loc] = (Float ) calloc (width SRC_SINC_DENSITY + 1, sizeof (Float));
243			sinc_widths[loc] = width;
244			size = width * SRC_SINC_DENSITY;
245			sinc_freq = M_PI / (Float) SRC_SINC_DENSITY;
246			win_freq = M_PI / (Float) size;
247			sinc_tables[loc][0] = 1.0;
248
249			for (i = 1, sinc_phase = sinc_freq, win_phase = win_freq; i < size; i++, sinc_phase += sinc_freq, win_phase += win_freq)
250			sinc_tables[loc][i] = sin (sinc_phase) * (0.5 + 0.5 * cos (win_phase)) / sinc_phase;
251
252			return (sinc_tables[loc]);
253			}
254
255			void
256			mus_src (Float input, int inpsize, Float output, int outsize, Float srate, Float *sr_mod, int width)
257			{
258			int i, lim, len, fsx, k, loc;
259
260			Float x, xx, factor, data, sinc_table, sum, zf, srx, incr;
261
262			Float *in0 = input;
263			Float *in1 = input + inpsize;
264			Float *out1 = output + outsize;
265
266			if (width == 0)
267			width = SRC_SINC_WIDTH;
268
269			x = 0.0;
270			lim = 2 * width;
271			len = width * SRC_SINC_DENSITY;
272			data = (Float *) calloc (lim + 1, sizeof (Float));
273			sinc_table = init_sinc_table (width);
274
275			for (i = width; i < lim; i++) data[i] = *input++;
276
277			while (output < out1)
278			{
279			fsx = (int)x;
280			if (fsx > 0)
281			{
282			/* realign data, reset x */
283			for (i = fsx, loc = 0; i < lim; i++, loc++)
284			data[loc] = data[i];
285
286			for (i = loc; i < lim; i++)
287			{
288			if (srx < 0)
289			input = (input > in0 ? input : in1) - 1;
290			else
291			input = input < in1 ? input+1 : in0;
292
293			data[i] = *input;
294			}
295
296			x -= fsx;
297			}
298
299			srx = srate + (sr_mod ? *sr_mod++ : 0);
300			srx = srx ? fabs (srx) : 0.001;
301			factor = srx > 1 ? 1 / srx : 1;
302
303			sum = 0.0;
304			zf = factor * SRC_SINC_DENSITY;
305			xx = zf * (1.0 - x - width);
306			for (i = 0; i < lim; i++)
307			{
308			/* we're moving backwards in the data array, so xx has to mimic that (hence the '1.0 - x') */
309			k = abs ((int)xx);
310
311			if (k < len)
312			sum += data[i] * sinc_table[k];
313
314			xx += zf;
315			}
316
317			x += srx;
318			output++ = sum factor;
319			}
320
321			free (data);
322			}
323
324			static unsigned long randx = 1;
325
326			static int
327			irandom (int amp)
328			{
329			int val;
330
331			randx = randx * 1103515245 + 12345;
332			val = (unsigned int) (randx >> 16) & 0x7fff;
333			return ((int) (amp * (((Float) val / 32768))));
334			}
335
336			#define max(a,b) ((a)>(b) ? (a) : (b))
337			#define min(a,b) ((a)<(b) ? (a) : (b))
338
339			void
340			mus_granulate (Float *input, int insize,
341			Float *output, int outsize,
342			Float expansion, Float flength, Float scaler,
343			Float hop, Float ramp, Float jitter, int max_size)
344			/* hop, jitter, length (= smapling_rate) /
345			{
346			int length = (int)ceil (flength);
347			int rmp = (int) (ramp * length);
348			int output_hop = (int)hop;
349			int input_hop = (int)(output_hop / expansion);
350			int s20 = (int) (jitter / 20);
351			int s50 = (int) (jitter / 50);
352			int outlen = max_size > 0 ? min ((int)(hop + flength), max_size) : (int)(hop + flength);
353			int in_data_len = outlen + s20 + 1;
354			int in_data_start = in_data_len;
355			Float data = (Float ) calloc (outlen, sizeof (Float));
356			Float in_data = (Float ) calloc (in_data_len, sizeof (Float));
357
358			Float *in1 = input + insize;
359			Float *out1 = output + outsize;
360
361			int ctr = 0;
362			Float cur_out = 0;
363
364			int start, len, end, i, j, k;
365			int steady_end, curstart;
366			Float incr, result, amp;
367
368			if (s50 > output_hop)
369			s50 = output_hop;
370
371			for(;;)
372			{
373			while (ctr < cur_out)
374			{
375			*output++ = data[ctr++];
376			if (output >= out1)
377			goto ok;
378			}
379
380			start = cur_out;
381			end = length - start;
382
383			if (end <= 0)
384			end = 0;
385			else
386			for (i = 0, j = start; i < end; i++, j++)
387			data[i] = data[j];
388
389			for (i = end; i < outlen; i++)
390			data[i] = 0;
391
392			start = in_data_start;
393			len = in_data_len;
394
395			if (start > len)
396			{
397			input += start - len;
398			input = input < in1 ? input : in1;
399			start = len;
400			}
401			else if (start < len)
402			for (i = 0, k = start; k < len; i++, k++)
403			in_data[i] = in_data[k];
404
405			for (i = len - start; i < len; i++)
406			{
407			in_data[i] = *input;
408			input = input < in1 ? input+1 : input;
409			}
410
411			in_data_start = input_hop;
412
413			amp = 0.0;
414			incr = scaler / (Float) rmp;
415			steady_end = length - rmp;
416			curstart = irandom (s20);
417
418			for (i = 0, j = curstart; i < length; i++, j++)
419			{
420			data[i] += (amp * in_data[j]);
421
422			if (i < rmp)
423			amp += incr;
424			else if (i > steady_end)
425			amp -= incr;
426			}
427
428			ctr -= cur_out;
429			cur_out = output_hop + irandom (s50) - (s50 >> 1);
430			}
431
432			ok:
433			free (data);
434			free (in_data);
435			}
436
437			static void
438			mus_shuffle (Float* rl, Float* im, int n)
439			{
440			/* bit reversal */
441
442			int i,m,j;
443			Float tempr,tempi;
444			j=0;
445			for (i=0;i<n;i++)
446			{
447			if (j>i)
448			{
449			tempr = rl[j];
450			tempi = im[j];
451			rl[j] = rl[i];
452			im[j] = im[i];
453			rl[i] = tempr;
454			im[i] = tempi;
455			}
456			m = n>>1;
457			while ((m>=2) && (j>=m))
458			{
459			j -= m;
460			m = m>>1;
461			}
462			j += m;
463			}
464			}
465
466			static void
467			mus_fft (Float rl, Float im, int n, int isign)
468			{
469			/* standard fft: real part in rl, imaginary in im, */
470			/* rl and im are zero-based. */
471			int mmax,j,pow,prev,lg,i,ii,jj,ipow;
472			Float wrs,wis,tempr,tempi;
473			double wr,wi,theta,wtemp,wpr,wpi;
474			ipow = (int)(ceil(log(n)/log(2.0)));
475			mus_shuffle(rl,im,n);
476			mmax = 2;
477			prev = 1;
478			pow = (int)(n*0.5);
479			theta = (M_PI*isign);
480			for (lg=0;lg<ipow;lg++)
481			{
482			wpr = cos(theta);
483			wpi = sin(theta);
484			wr = 1.0;
485			wi = 0.0;
486			for (ii=0;ii<prev;ii++)
487			{
488			wrs = (Float) wr;
489			wis = (Float) wi;
490			#ifdef LINUX
491			if (isnan(wis)) wis=0.0;
492			#endif
493			i = ii;
494			j = ii + prev;
495			for (jj=0;jj<pow;jj++)
496			{
497			tempr = wrsrl[j] - wisim[j];
498			tempi = wrsim[j] + wisrl[j];
499			rl[j] = rl[i] - tempr;
500			im[j] = im[i] - tempi;
501			rl[i] += tempr;
502			im[i] += tempi;
503			i += mmax;
504			j += mmax;
505			}
506			wtemp = wr;
507			wr = (wrwpr) - (wiwpi);
508			wi = (wiwpr) + (wtempwpi);
509			}
510			pow = (int)(pow*0.5);
511			prev = mmax;
512			theta = theta*0.5;
513			mmax = mmax*2;
514			}
515			}
516
517			static void
518			mus_convolution (Float* rl1, Float* rl2, int n)
519			{
520			/* convolves two real arrays. */
521			/* rl1 and rl2 are assumed to be set up correctly for the convolution */
522			/* (that is, rl1 (the "signal") is zero-padded by length of */
523			/* (non-zero part of) rl2 and rl2 is stored in wrap-around order) */
524			/* We treat rl2 as the imaginary part of the first fft, then do */
525			/* the split, scaling, and (complex) spectral multiply in one step. */
526			/* result in rl1 */
527
528			int j,n2,nn2;
529			Float rem,rep,aim,aip,invn;
530
531			mus_fft(rl1,rl2,n,1);
532
533			n2=(int)(n*0.5);
534			invn = 0.25/n;
535			rl1[0] = ((rl1[0]*rl2[0])/n);
536			rl2[0] = 0.0;
537
538			for (j=1;j<=n2;j++)
539			{
540			nn2 = n-j;
541			rep = (rl1[j]+rl1[nn2]);
542			rem = (rl1[j]-rl1[nn2]);
543			aip = (rl2[j]+rl2[nn2]);
544			aim = (rl2[j]-rl2[nn2]);
545
546			rl1[j] = invn(repaip + aim*rem);
547			rl1[nn2] = rl1[j];
548			rl2[j] = invn(aimaip - rep*rem);
549			rl2[nn2] = -rl2[j];
550			}
551
552			mus_fft(rl1,rl2,n,-1);
553			}
554
555			void
556			mus_convolve (Float * input, Float * output, int size, Float * filter, int fftsize, int filtersize)
557			{
558			int fftsize2 = fftsize >> 1;
559			Float rl1, rl2, *buf;
560			Float *in1 = input + size;
561			int ctr = fftsize2;
562			int i, j;
563
564			rl1 = (Float *) calloc (fftsize, sizeof (Float));
565			rl2 = (Float *) calloc (fftsize, sizeof (Float));
566			buf = (Float *) calloc (fftsize, sizeof (Float));
567
568			while (size > 0)
569			{
570			ctr++;
571			if (ctr >= fftsize2)
572			{
573			for (i = 0, j = fftsize2; i < fftsize2; i++, j++)
574			{
575			buf[i] = buf[j];
576			buf[j] = 0.0;
577			rl1[i] = *input;
578			rl1[j] = 0.0;
579			rl2[i] = 0.0;
580			rl2[j] = 0.0;
581
582			input = input < in1 ? input+1 : input;
583			}
584
585			for (i = 0; i < filtersize; i++)
586			rl2[i] = filter[i];
587
588			mus_convolution (rl1, rl2, fftsize);
589
590			for (i = 0, j = fftsize2; i < fftsize2; i++, j++)
591			{
592			buf[i] += rl1[i];
593			buf[j] = rl1[j];
594			}
595
596			ctr = 0;
597			}
598
599			*output++ = buf[ctr];
600			size--;
601			}
602			}