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 (21 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 *
#	Content
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	}