This is cmus.c in view mode; [Download] [Up]
/* adjunct to cmus.lisp
*
* the simple cases, both arithmetic and unit-generator, are expanded in-line.
*
*/
#include <stdio.h>
#include <math.h>
#include "cmus.h"
extern int clm_read(int fd, char *buf, int n);
extern int clm_write(int fd, char *buf, int n);
extern short clm_short(short n);
extern int clm_int(int n);
extern int clm_create(char *name);
extern long clm_seek(int tfd, long offset, int origin);
/* io file structure reflected in datai */
#define io_fd 0
#define io_chans 1
#define io_size 2
#define io_beg 3
#define io_end 4
#define io_bufsiz 5
#define io_data_start 6
#define io_data_end 7
#define io_open_index 8
#define io_dat_A 9
#define io_dat_B 10
#define io_dat_C 11
#define io_dat_D 12
#define io_dir 13
#define io_loc 14
#define io_hdr_end 15
#define c_aref_block 0
#define c_aref_type 1
#define c_aref_size 2
#define c_aref_dims 3
#define c_aref_dim_list_adr 4
#define io_in_f 0
#define io_out_f 1
/* from io.c */
extern cwritez(int fd, int num);
extern crdmono(int fd, int beg, int end, int *buf);
extern cwrtmono(int fd, int beg, int end, int *unpackedbuf);
extern crdstereo(int fd, int beg, int end, int *bufa, int *bufb);
extern cwrtstereo(int fd, int beg, int end, int *unpackedbufa, int *unpackedbufb);
extern crdquad(int fd, int beg, int end, int *bufa, int *bufb, int *bufc, int *bufd);
extern cwrtquad(int fd, int beg, int end, int *unpackedbufa, int *unpackedbufb, int *unpackedbufc, int *unpackedbufd);
/* io support functions translated from io.lisp */
void c_check_stream_indices (int *io)
{
if ((io[io_data_start] < 0) || (io[io_data_end] < 0))
{
printf("file buffer index is negative! -- will try to fix it..."); fflush(stdout);
if (io[io_data_start] < 0) io[io_data_start] = 0; else io[io_data_end] = 0;
}
}
void c_put_data (int *io)
{
if (io[io_fd] != -1) /* -1 -> closed file */
{
c_check_stream_indices(io);
switch (io[io_chans])
{
case 0: printf("aw good grief -- no IO buffers!"); fflush(stdout); break;
case 1: cwrtmono(io[io_fd],io[io_data_start],io[io_data_end],(int *)io[io_dat_A]); break;
case 2: cwrtstereo(io[io_fd],io[io_data_start],io[io_data_end],(int *)io[io_dat_A],(int *)io[io_dat_B]); break;
case 4: cwrtquad(io[io_fd],io[io_data_start],io[io_data_end],(int *)io[io_dat_A],
(int *)io[io_dat_B],(int *)io[io_dat_C],(int *)io[io_dat_D]); break;
}
}
}
void c_get_data (int *io, int size)
{
if (size > io[io_bufsiz])
{
printf("input request is too big: %d", size); fflush(stdout);
}
else
{
if (size == 0)
{
printf("get-data asked to read 0 bytes!?"); fflush(stdout);
}
else
{
c_check_stream_indices(io);
switch (io[io_chans])
{
case 0: printf("aw ferchrissake -- no IO channels!"); fflush(stdout); break;
case 1: crdmono(io[io_fd],0,size-1,(int *)io[io_dat_A]); break;
case 2: crdstereo(io[io_fd],0,size-1,(int *)io[io_dat_A],(int *)io[io_dat_B]); break;
case 4: crdquad(io[io_fd],0,size-1,(int *)io[io_dat_A],
(int *)io[io_dat_B],(int *)io[io_dat_C],(int *)io[io_dat_D]); break;
}
}
}
}
int c_sample_position(int *io, int loc)
{
clm_seek(io[io_fd],io[io_hdr_end]+(2*io[io_chans]*loc),0);
return(loc);
}
void c_io_bufclr (int *io, int beg)
{
if (io[io_dat_A] == 0) {printf("attempt to clear deallocated IO buffer"); fflush(stdout);}
else
{
int i,end;
int *j;
end=io[io_bufsiz];
j=((int *)io[io_dat_A]);
for (i=beg; i<end;i++) j[i]=0;
if (io[io_chans] > 1)
{
j = (int *)(io[io_dat_B]);
for (i=beg; i<end;i++) j[i]=0;
if (io[io_chans] > 2)
{
j = (int *)(io[io_dat_C]);
for (i=beg; i<end;i++) j[i]=0;
j = (int *)(io[io_dat_D]);
for (i=beg; i<end;i++) j[i]=0;
}
}
}
}
void c_flush_snd(int *io)
{
if ((io[io_fd] != -1) &&
(io[io_dir] != io_in_f) &&
(io[io_data_end] != 0) &&
(io[io_data_start] != io[io_data_end]))
{
if (io[io_data_end] > io[io_bufsiz])
{
printf("data end indication in IO buffer is too big: %d > %d",io[io_data_end],io[io_bufsiz]); fflush(stdout);
}
else
{
c_sample_position(io,io[io_beg]+io[io_data_start]);
c_put_data(io);
io[io_data_start] = io[io_data_end];
if (io[io_size] < (io[io_beg] + io[io_data_end])) io[io_size] = (io[io_beg] + 1 + io[io_data_end]);
}
}
}
void clm_file_reset(int loc0, int *io)
{
/* called when loc is outside the current in-core frame for the file pointed to by io */
/* equivalent to the io.lisp function file-check + read-in */
int file_end,bytes,loc;
loc = loc0;
if (io[io_dir] != io_in_f) c_flush_snd(io);
if ((loc < io[io_beg]) && ((loc + (int)(.9*io[io_bufsiz])) > io[io_beg]))
{
if ((loc + 10) > io[io_beg]) loc -= (int)(.75*io[io_bufsiz]);
if (loc < 0) loc = 0;
if (io[io_chans] == 1) loc = (2 * (int)(loc / 2));
}
file_end = io[io_size];
bytes = file_end - loc;
if (bytes > io[io_bufsiz]) bytes=io[io_bufsiz];
if (bytes < 0) /* tried to access beyond current end of file */
{
if (io[io_dir] == io_in_f) {printf("attempt to read past EOF"); fflush(stdout);}
else
{
c_io_bufclr(io,0);
bytes = io[io_bufsiz];
io[io_data_start] = 0;
io[io_data_end] = 0;
clm_seek(io[io_fd],0,2); /* go to end of file */
switch (io[io_chans])
{
case 4: cwritez(io[io_fd],2*(loc-file_end)); io[io_beg]=loc; break;
case 2: cwritez(io[io_fd],loc-file_end); io[io_beg]=loc; break;
case 1: cwritez(io[io_fd],(int)((loc-file_end)/2)); if ((loc%2)==0) io[io_beg]=loc; else io[io_beg]=loc-1; break;
}
}
}
else /* bytes is positive or 0 */
{
io[io_beg] = c_sample_position(io,loc);
if (bytes > 0) c_get_data(io,bytes);
if (bytes < io[io_bufsiz]) c_io_bufclr(io,bytes);
io[io_data_start] = 0;
io[io_data_end] = 0;
if (io[io_dir] == io_in_f)
{
int bufend,filbytes;
bufend = io[io_bufsiz]-1;
filbytes = file_end-io[io_beg]-1;
if (filbytes < bufend) bufend = filbytes;
if (bufend > 0) io[io_data_end] = bufend;
}
}
io[io_end] = io[io_beg]+io[io_bufsiz]-1;
/* io[io_data_start] = io[io_data_end]; */
io[io_loc] = (loc0-io[io_beg]);
}
void clm_full_file_reset(int loc, int *io)
{
/* update data_end before calling clm_file_reset -- we are about to write at loc, so loc-1 is last output */
/* dangerous assumption!? */
int lim;
lim = ((loc - io[io_beg]) - 1);
if (lim >= io[io_bufsiz]) lim = (io[io_bufsiz]-1);
if (lim > io[io_data_end]) io[io_data_end] = lim;
clm_file_reset(loc,io);
}
void c_newline(void) {printf("\n"); fflush(stdout);}
/* sine wave via table lookup -- normally faster than the built-in sin call */
/* on the Mac, however, the goddam table won't fit */
#ifndef MAC
static float c_sine_table[514];
static float c_sine_table_scaler = (float)(512.0/two_pi);
static int c_sine_table_is_inited = 0;
#endif
c_init_sine (void)
{
#ifndef MAC
int i;
float phase,incr;
if (c_sine_table_is_inited == 0)
{
incr = (float)1.0/c_sine_table_scaler;
for (i=0,phase=0.0;i<514;i++,phase+=incr)
{
c_sine_table[i]=(float)sin((double)phase);
}
c_sine_table_is_inited = 1;
}
#endif
}
float c_sin_lookup(float phase)
{
#ifndef MAC
/* in-coming phase is in radians, can be anything */
/* use of sin here was noticeably slower on the 68040 */
float mag_phase,frac;
int index;
float *sint;
if (phase < 0.0)
{
index=(int)(phase*inverse_pi);
phase+=(float)(((1-index)*(two_pi)));
}
mag_phase = phase*c_sine_table_scaler;
index = (int)mag_phase;
frac = mag_phase - index;
index = index%512; /* % is mod */
sint = c_sine_table+index;
return ((*sint + (frac*(*(sint+1)-*sint))));
#else
return(sin(phase));
#endif
}
#ifndef MAC
/* from Cmix -- appears to be more reliable than the 20 bit form */
static long randx = 1;
#define inverse_rscl .000061035156
float c_frandom(float amp)
{
int i = ((randx = randx*1103515245 + 12345)>>16) & 077777;
return(amp * (((float)i)*inverse_rscl - 1.0));
}
#else
float c_frandom(float amp)
{
int i;
i = rand();
return(amp * ((float)(i/16384.0) - 1.0));
}
#endif
void c_fft_window(float *rl, float *win, int size) {int i; for (i=0;i<size;i++) rl[i]=(rl[i]*win[i]);}
void c_clear_block(float *rl, int size) {int i; for (i=0;i<size;i++) rl[i]=0.0;}
void c_blti(int *arr0, int beg, int *arr1, int len) {int i,k; for (i=0,k=beg;i<len;i++,k++) arr0[k]=arr1[i];}
void c_bltf(float *arr0, int beg, float *arr1, int len) {int i,k; for (i=0,k=beg;i<len;i++,k++) arr0[k]=arr1[i];}
double c_gcd(double a, double b)
{
/* everything double here to mimic main c math libraries */
double rem;
if (a == b) return(a);
else
{
if (a>b)
{
rem = fmod(a,b);
if (rem == 0.0) return(b);
else return(c_gcd(b,rem));
}
else return(c_gcd(b,a));
}
}
double c_lcm(double a, double b)
{
return((a*b)/c_gcd(a,b));
}
float c_table_interp (float *fn, float xx, int len)
{
int int_part;
float frac_part,val0,val1;
if (xx < 0.0)
{
int_part = (int)(xx/len);
xx += ((1-int_part)*len);
}
int_part = (int)(xx);
frac_part = xx-(float)int_part;
if (int_part >= len) int_part = (int_part%len);
if (frac_part == 0.0) return(fn[int_part]);
val0 = fn[int_part];
if (int_part == (len-1)) val1 = fn[0];
else val1 = fn[int_part+1];
return (val0 + (frac_part*(val1-val0)));
}
float c_table_lookup (int *tbl, float *datar, float fm)
{
float *wave;
float val,freq,phase;
int size;
freq = datar[tbl[0]];
phase = datar[tbl[0]+1];
wave = (float *)(datar+tbl[1]);
size = tbl[1+c_aref_size];
val = c_table_interp(wave,phase,size);
phase += (freq + fm);
while (phase >= size) phase -= size;
while (phase < 0.0) phase += size;
datar[tbl[0]+1] = phase;
return(val);
}
float c_direct_flt (int m, float input, float *flt_a, float *flt_b, float *flt_d)
{
float xout;
int j;
xout = 0.0;
flt_d[0] = input;
for (j=m;j>=1;j--)
{
xout += flt_d[j]*flt_b[j];
flt_d[0] -= flt_a[j]*flt_d[j];
flt_d[j] = flt_d[j-1];
}
return(xout+(flt_d[0]*flt_b[0]));
}
float c_lattice_flt (int m, float input, float *flt_a, float *flt_b, float *flt_d)
{
float xout,inp;
int i,j;
xout = 0.0;
inp = input;
for (i=m-1,j=m;i>=0;i--,j--)
{
inp -= flt_a[i]*flt_d[i];
flt_d[j] = flt_d[i] + (flt_a[i]*inp);
xout += flt_d[j]*flt_b[j];
}
flt_d[0] = inp;
return(xout + (inp*flt_b[0]));
}
float c_ladder_flt (int m, float input, float so, float *flt_a, float *flt_b, float *flt_c, float *flt_d)
{
float xout,inp;
int i;
xout = 0.0;
inp = input;
for (i=m-1;i>=0;i--)
{
flt_d[i+1] = (inp*flt_a[i]) + (flt_d[i]*flt_c[i]);
inp = (inp*flt_c[i]) - (flt_d[i]*flt_a[i]);
xout += flt_d[i+1]*flt_b[i+1];
}
flt_d[0] = inp;
return(so*(xout+(flt_d[0]*flt_b[0])));
}
int c_y_or_n_p (void)
{
int c,i;
char s[256];
/* question string has already gone out */
LOOP:
i=0;
while ((i<256) && ((c = getchar()) != '\n') && (c != EOF)) {s[i]=c; i++;}
if ((s[0] == 'y') || (s[0] == 'Y')) return(1);
if ((s[0] == 'n') || (s[0] == 'N')) return(0);
printf("Type \"y\" for yes or \"n\" for no. "); fflush(stdout);
goto LOOP;
}
int c_yes_or_no_p (void)
{
int c,i;
char s[256];
/* question string has already gone out */
LOOP:
i=0;
while ((i<256) && ((c = getchar()) != '\n') && (c != EOF)) {s[i]=c; i++;}
if (((s[0] == 'y') || (s[0] == 'Y')) && ((s[1] == 'e') || (s[1] == 'E')) && ((s[2] == 's') || (s[2] == 'S')) && (i == 3)) return(1);
if (((s[0] == 'n') || (s[0] == 'N')) && ((s[1] == 'o') || (s[1] == 'O')) && (i == 2)) return(0);
printf("Type \"yes\" or \"no\". "); fflush(stdout);
goto LOOP;
}
/* fft and convolution of real data in zero-based arrays */
/* plus code to tie into clm's convolve generator */
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;
}
}
c_fft (float* rl, float* im, int n, int isign, int ipow)
{
/* standard fft: real part in rl, imaginary in im, */
/* ipow = ceiling (log n / log 2), isign=1 fft, =-1 ifft. */
/* rl and im are zero-based. */
/* */
/* oddly enough, the integer version (using integer ops */
/* and "block floating point" scaling) was much slower, */
/* and splitting out the no-ops (twiddle factors==0 etc) */
/* made no difference at all. */
int mmax,j,pow,prev,lg,i,ii,jj;
float wrs,wis,tempr,tempi;
double wr,wi,theta,wtemp,wpr,wpi;
shuffle(rl,im,n);
mmax = 2;
prev = 1;
pow = n*0.5;
theta = (one_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;
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 = pow*0.5;
prev = mmax;
theta = theta*0.5;
mmax = mmax*2;
}
}
convolve (float* rl1, float* rl2, int n, int ipow)
{
/* 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;
c_fft(rl1,rl2,n,1,ipow);
n2=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];
}
c_fft(rl1,rl2,n,-1,ipow);
}
/* now functions to read/write linear 16-bit 1/2 channel sound files converting to/from floats */
floadarray (float* rl0, int rlbeg, int beg, int end, int step, char* data)
{
char *k;
float *x;
float rscl;
int i,sample;
rscl = 1.0/32678.0;
for (i=beg,k=data+beg,x=rl0+rlbeg;i<end;i++,x++,k+=step)
{
sample = ( (((*k)<<8)&0xff00) | ((*(k+1))&0xff));
if ((*k) & 0x8000) sample |= 0xffff0000;
*x = (float) sample*rscl;
}
}
fwritearray(int chan, int wrds, float* rl0, int rlbeg, char* buffer)
{
char *k;
float *x;
int i,sample;
for (i=0,k=buffer,x=rl0+rlbeg;i<wrds;i++,x++,k+=2)
{
sample = (int) 32768*(*x);
*(k+1)=(sample&0xff);
*k=((sample>>8)&0xff);
}
clm_write(chan,buffer,(wrds*2));
}
/* convolve either all at once (ipow==allpow) or by overlap-add */
/* filter is assumed to be set up correctly (i.e. in wrap-around order if that is what is wanted) */
/* but in any case, we have to re-set it up in the new (assured-to-be-power-of-2) buffer */
void c_convolve (char *fname,
int filec, int filehdr, int filesize,
int filterc, int filterhdr, int filtersize,
int start, int fftsize, int ipow, int allpow, int samps, int fullsize, int chans, int chan)
{
#ifndef AKCL
extern void write_next_header (int chan, int srate, int chans, int loc, int siz, int format, char *comment, int len);
#endif
#ifdef MAC
#pragma unused(filesize,chan)
#endif
float *rl0,*rl1,*temprl,*tempout;
char *buffer;
int tempfile,i,readsize,startloc,j,step,fsize,filterwrds,sampctr,fft2,fft4;
fft2 = (int) fftsize*.5;
fft4 = (int) fft2*.5;
step = 2;
if (chans == 2) step = 4;
startloc = filehdr + start*step;
readsize = fftsize*step;
if (allpow != ipow) readsize = (int) readsize*.5;
fsize = 2*fullsize;
/* get to start point in the two sound files */
lseek(filec,startloc,0);
lseek(filterc,filterhdr,0);
rl0 = (float *)calloc(fftsize,sizeof(float));
rl1 = (float *)calloc(fftsize,sizeof(float));
buffer = (char *)calloc(readsize,sizeof(char));
/* read in the filter data */
clm_read(filterc,buffer,filtersize*2);
floadarray(rl1,0,0,filtersize,2,buffer);
/* make sure filter (i.e. impulse response) is in wrap around order in the new array */
if (filtersize != fftsize)
{
filterwrds = (int) filtersize*.5;
for (i=(fftsize-1),j=(filtersize-1);j>filterwrds;i--,j--)
{
rl1[i] = rl1[j];
rl1[j] = 0.0;
}
}
if (allpow != ipow)
{
/* save filter samples -- rl1 will be clobbered by the convolution */
temprl = (float *)calloc(fftsize,sizeof(float));
for (i=0;i<fftsize;i++) temprl[i] = rl1[i];
/* if overlap-add (since ipow != allpow) create buffer to hold intermediate results */
tempout = (float *)calloc(fftsize,sizeof(float));
}
/* need file to hold convolution output */
tempfile=clm_create(fname);
write_next_header(tempfile,22050,1,28,fsize,1,"",0);
/* split out the easy case (where only one convolution is performed) */
if (allpow == ipow)
{
clm_read(filec,buffer,readsize);
floadarray(rl0,0,2*(chans-1),fftsize,step,buffer);
convolve(rl0,rl1,fftsize,ipow);
fwritearray(tempfile,fftsize,rl0,0,buffer);
}
else
{
/* the messy case -- here we loop through the signal doing overlap-adds */
for (sampctr=0;sampctr<samps;sampctr+=fft2)
{
clm_read(filec,buffer,readsize);
for (i=0;i<fft4;i++) rl0[i] = 0.0;
floadarray(rl0,fft4,2*(chans-1),fft2,step,buffer);
for (i=fft2+fft4;i<fftsize;i++) rl0[i] = 0.0;
if (sampctr != 0)
{
for (i=0;i<fftsize;i++) rl1[i] = temprl[i];
}
convolve(rl0,rl1,fftsize,ipow);
for (i=0;i<fftsize;i++) tempout[i] += rl0[i];
if (sampctr != 0)
{
fwritearray(tempfile,fft2,tempout,0,buffer);
}
else
{
fwritearray(tempfile,fft4,tempout,fft4,buffer);
}
for (i=0,j=fft2;i<fft2;i++,j++)
{
tempout[i] = tempout[j];
tempout[j] = 0.0;
}
}
fwritearray(tempfile,fft4,tempout,fft2,buffer);
free(temprl);
free(tempout);
}
close(tempfile);
free(rl0);
free(rl1);
free(buffer);
}
#define env_current_value 0
#define env_rate 1
#define env_base 2
#define env_offset 3
#define env_scaler 4
#define env_power 5
#define env_init_y 6
#define env_init_power 7
#define env_type 0
#define env_index 1
#define env_top 2
#define env_pass 3
#define env_end 4
#define env_restart 5
void c_env_restart (float *ef, int *ei, float *data)
{
int diff,i;
ef[env_current_value] = ef[env_init_y];
ef[env_power] = ef[env_init_power];
ef[env_rate] = 0.0;
ei[env_index] = 0;
diff = (ei[env_pass] - ei[env_restart]);
ei[env_restart] = ei[env_pass];
ei[env_end] += diff;
for (i=0;i<ei[env_top];i+=2) data[i] += diff;
}
#ifdef MAC
float c_env (float *ef, int *ei, float *data)
{
int loc,passpt;
float val;
loc = ei[env_index];
val = ef[env_current_value];
passpt = data[loc]; /* believe it or not, (int) here, which goddam well ought to be a no-op causes the */
/* the result to be 4 times the true value as an integer! There are lots more bugs */
/* along these lines in MPW C -- surely the worst implementation of C on the planet */
if ((loc < ei[env_top]) && (ei[env_pass] >= passpt))
{
ef[env_rate] = data[loc+1];
ei[env_index] += 2;
}
ei[env_pass] += 1;
if (ei[env_pass]<=ei[env_end])
{
if (ei[env_type] == 0)
{ /* line-segment or step envelope */
if (ef[env_base] != 0.0)
ef[env_current_value] += ef[env_rate];
else ef[env_current_value] = ef[env_rate];
}
else /* exponential envelope */
{
if (ef[env_rate] != 0.0)
{
ef[env_power] += ef[env_rate];
ef[env_current_value] = (ef[env_offset] + (ef[env_scaler] * (exp(ef[env_base]*ef[env_power]) - 1.0)));
}
}
}
return(val);
}
#endif
#ifndef __LITTLE_ENDIAN__
#ifndef MAC
/* exp from Jos Vermaseren */
/* on the NeXT this code is about 4 times as fast as the equivalent call on the libm version of exp */
/* on little endian machines this code dies because it assumes big endian IEEE format */
double jv_exp(double x)
{
/*
routine for the evaluation of the exp (IEEE)
Method:
1: y = x/log(2) -> exp(x) = 2^y;
2: z = n+y, -0.5 <= y < 0.5
3: exp(x) = exp(z*log(2))*2^n
For the last series we need only terms till 13!
Routine made by J.A.M.Vermaseren 20-apr-1992
*/
double y,z,scr;
long n;
int m,sflag;
y = x*(1.0/0.6931471805599453094);
if ( y < 0 ) { sflag = 1; y = -y; }
else sflag = 0;
scr = y + 0.5;
n = *((long *)(&scr));
m = (n >> 20) - 0x3FF;
if ( m > 7 ) {
fprintf(stderr,"Overflow in exponent of %e\n",x);
exit(-1);
}
n &= 0xFFFFF;
n |= 0x100000;
n >>= (20-m);
z = (y - n)*0.6931471805599453094;
if ( sflag ) { n = -n; z = -z; }
scr = 0.99999999999999999693+z*(
0.99999999999999999850+z*(
0.50000000000000183922+z*(
0.16666666666666701804+z*(
0.04166666666648804594+z*(
0.00833333333330984930+z*(
0.00138888889523274996+z*(
0.00019841269908532560+z*(
0.00002480148546912669+z*(
0.00000275572255877329+z*(
0.00000027632644397330+z*(
0.00000002511466084085)))))))))));
n <<= 20;
*((long *)(&scr)) += n;
return(scr);
}
#else
double jv_exp(double x) {return(exp(x));}
#endif
#else
double jv_exp(double x) {return(exp(x));}
#endif
/* RUN-BLOCK */
#define blk_loc 0
#define blk_siz 1
#define blk_ctr 0
float c_run_block (int *blk, float *fblk, float *buf)
{
int loc,i,j;
float val;
loc = blk[blk_loc];
if (loc < blk[blk_siz]) val = buf[loc];
else val = 0.0;
loc++;
if (loc >= fblk[blk_ctr])
{
if (loc < blk[blk_siz])
{
for (i=0;i<loc;i++) buf[i]=0.0;
for (i=0,j=loc;j<blk[blk_siz];i++,j++)
{
buf[i]=buf[j];
buf[j]=0.0;
}
}
else
{
for (i=0;i<blk[blk_siz];i++) buf[i]=0.0;
}
fblk[blk_ctr] -= loc;
blk[blk_loc] = 0;
}
else blk[blk_loc] = loc;
return(val);
}
void c_base_run_block (int *blk, float *fblk, float *buf)
{
int loc,i,j;
loc = blk[blk_loc];
if (loc < blk[blk_siz])
{
for (i=0;i<loc;i++) buf[i]=0.0;
for (i=0,j=loc;j<blk[blk_siz];i++,j++)
{
buf[i]=buf[j];
buf[j]=0.0;
}
}
else
{
for (i=0;i<blk[blk_siz];i++) buf[i]=0.0;
}
fblk[blk_ctr] -= loc;
blk[blk_loc] = 0;
}
/* READIN */
float c_readin (int *rd, int *io)
{
float val;
if ((rd[0] >= 0) && (rd[0] < io[io_size]))
{
if ((rd[0] < io[io_beg]) || (rd[0] > io[io_end])) clm_file_reset(rd[0],io);
val = (float)(_sndflt_ * (((int *)(io[io_dat_A+rd[2]]))[rd[0] - io[io_beg]]));
}
else val = 0.0;
rd[0] += rd[1];
return(val);
}
/* SRC */
float sinc(float x)
{
float pi_x;
if (x == 0.0)
return 1.0;
else
{
pi_x = x * one_pi;
return(c_sin_lookup(pi_x) / pi_x);
}
}
float smoothed_table (float *data, float x, int size)
{
int i;
float sum;
if (size == 1)
return x;
else
{
size -= 1;
for (i=0;i<size;i++) data[i]=data[i+1];
data[size]=x;
sum = 0.5*(data[0]+x);
for (i=1;i<size;i++) sum += data[i];
return (sum/(float)size);
}
}
#define src_x 0
#define src_incr 1
#define src_offset 2
#define src_left 0
#define src_right 1
#define src_width 2
#define src_use_filter 3
#define src_filter_size 4
float c_src(int *rd, int *rio, float *srf, int *sri, float *data, float *filt, float fm)
{
float sum,start_x;
int loc,lim,fsx,i,k,right;
right=sri[src_right];
lim=2*sri[src_width];
start_x=srf[src_x]+srf[src_offset]-sri[src_width];
fsx=0;
loc=0;
if (start_x > 1.0) fsx=(int)start_x;
if ((sri[src_left] < fsx) || ((fsx+lim) > right))
{
if (fsx <= right)
{
for (i=fsx,k=(lim-(right-fsx));i<=right;i++,k++)
{
data[loc]=data[k];
loc++;
}
}
else
if (right != -1) rd[0] += (fsx-right); /* just flush these samples */
if (sri[src_use_filter] != 0)
{
for (i=loc;i<=lim;i++) data[i]=smoothed_table(filt,c_readin(rd,rio),sri[src_filter_size]);
}
else
for (i=loc;i<=lim;i++) data[i]=c_readin(rd,rio);
}
sri[src_left]=fsx;
sri[src_right]=fsx+lim;
start_x=(sri[src_left]-srf[src_x])-srf[src_offset];
sum = 0.0;
for (i=0;i<=lim;i++)
{
sum += data[i]*sinc(start_x);
start_x+=1.0; /* bizarre -- this is implicit in the mus.lisp loop */
/* the 56k version windows the sinc with a hamming window */
}
srf[src_x] += (fm+srf[src_incr]);
if (srf[src_x] > 1024.0) {srf[src_offset] += 1024.0; srf[src_x] -= 1024.0;}
return(sum);
}
/* EXPAND */
#define spd_len 0
#define spd_rmp 1
#define spd_in_spd 2
#define spd_out_spd 3
#define spd_cur_in 4
#define spd_cur_out 5
#define spd_s20 6
#define spd_s50 7
#define spd_ctr 8
#define spd_amp 0
void shift_block_back (int *spidat, float *b)
{
int start,end,i,j;
start=spidat[spd_cur_out];
end=spidat[spd_len]-start;
for (i=0,j=start;i<end;i++,j++) b[i]=b[j];
for (i=end;i<spidat[spd_len];i++) b[i]=0.0;
}
void add_one_segment (int *rd, int *spidat, float *spfdat, int *rio, float *b)
{
float curamp,incr;
int i,k,steady_time,ramp_time;
curamp = 0.0;
ramp_time = spidat[spd_rmp];
incr = spfdat[spd_amp] / (float)ramp_time;
steady_time = spidat[spd_len] - (2*ramp_time);
k=0;
for (i=0;i<ramp_time;i++,k++)
{
b[k] += (curamp * c_readin(rd,rio));
curamp += incr;
}
for (i=0;i<steady_time;i++,k++)
{
b[k] += (curamp * c_readin(rd,rio));
}
for (i=0;i<ramp_time;i++,k++)
{
b[k] += (curamp * c_readin(rd,rio));
curamp -= incr;
}
}
void set_expansion_triggers (int *rddat, int *spidat)
{
spidat[spd_ctr] -= spidat[spd_cur_out];
spidat[spd_cur_out] = (spidat[spd_out_spd] + (int)c_frandom((float)spidat[spd_s20]));
rddat[0] = spidat[spd_cur_in] + (int)c_frandom((float)spidat[spd_s50]);
spidat[spd_cur_in] += spidat[spd_in_spd];
}
float c_spd (int *rddat, int *spidat, float *spfdat, int *io, float *b)
{
float val;
val=b[spidat[spd_ctr]];
spidat[spd_ctr]++;
if (spidat[spd_ctr] >= spidat[spd_cur_out])
{
shift_block_back(spidat,b);
add_one_segment(rddat,spidat,spfdat,io,b);
set_expansion_triggers(rddat,spidat);
}
return(val);
}
/* FFT-FILTER */
#define fftf_siz 0
#define fftf_hop 1
#define fftf_half_siz 2
#define fftf_pow 3
void readin_data (int *rddat, int *fftdat, int *rio, float *datar, float *datai)
{
int start,amount,i,k;
float scaler;
start=(int)(fftdat[fftf_siz]/4);
amount=fftdat[fftf_half_siz];
scaler = 1.0/(float)fftdat[fftf_siz];
for (i=0,k=start+amount;i<start;i++,k++)
{
datar[i]=0.0;
datai[i]=0.0;
datar[k]=0.0;
datai[k]=0.0;
}
if ((rio != NULL) && (rio[io_fd] != -1))
{
for (k=0,i=start;k<amount;k++,i++)
{
datar[i]=(scaler * c_readin(rddat,rio));
datai[i]=0.0;
}
}
else
{
for (k=0,i=start;k<amount;k++,i++)
{
datar[i]=c_frandom(scaler);
datai[i]=0.0;
}
}
}
void apply_envelope_and_shift(int *fftdat, float *datar, float *datai, float *fenv)
{
int i,k,lim,mid;
float val;
lim=fftdat[fftf_siz]-1;
mid=fftdat[fftf_half_siz];
for (i=0,k=lim;i<mid;i++,k--)
{
val=fenv[i];
datar[i]=(datar[i]*val);
datai[i]=(datai[i]*val);
datar[k]=(datar[k]*val);
datai[k]=(datai[k]*val);
}
}
void readout_first_data(int *fftdat, float *runbuf, float *datar)
{
int start,amount,i,k;
start=(int)(fftdat[fftf_siz]/4);
amount=start*3;
for (i=0,k=start;i<amount;i++,k++)
{
runbuf[i]=datar[k];
}
}
void readout_data(int *fftdat, float *runbuf, float *datar)
{
int stop,i,k;
stop=fftdat[fftf_hop];
for (i=0;i<stop;i++) runbuf[i] += datar[i];
for (i=0,k=stop;i<stop;i++,k++) runbuf[k] = datar[k];
}
float c_fft_filter (int *rddat, int *io, int *runblk, float *fblk, int *fftdat,
float *runbuf, float *datar, float *datai, float *fenv)
{
if (runblk[blk_loc] == 0)
{
readin_data(rddat,fftdat,io,datar,datai);
c_fft(datar,datai,fftdat[fftf_siz],1,fftdat[fftf_pow]);
apply_envelope_and_shift(fftdat,datar,datai,fenv);
c_fft(datar,datai,fftdat[fftf_siz],-1,fftdat[fftf_pow]);
if (fftdat[fftf_hop] < 0)
{
fftdat[fftf_hop] = fftdat[fftf_half_siz];
fblk[blk_ctr] = (float)(fftdat[fftf_siz]/4);
readout_first_data(fftdat,runbuf,datar);
}
else
{
readout_data(fftdat,runbuf,datar);
fblk[blk_ctr] = (float)fftdat[fftf_hop];
}
}
return(c_run_block(runblk,fblk,runbuf));
}
load_one_sine (float partial, float amp, float *table, int table_size, float phase)
{
float angle,freq;
int i;
if (amp != 0.0)
{
freq = (partial*two_pi)/((float)table_size);
for (i=0,angle=phase;i<table_size;i++,angle+=freq)
table[i]+=amp*sin(angle);
}
}
float c_convolve_or_readin (int *rddat, int *io, int *runblk, float *fblk, int *fftdat,
float *runbuf, float *datar, float *datai, float *filtr, float *filti)
{
/* can be just a readin call if two very large files were convolved into a temp */
if (filtr == NULL)
{
return(c_readin(rddat,io));
}
else
{
if (runblk[blk_loc] == 0)
{
float scale,temp;
int i;
readin_data(rddat,fftdat,io,datar,datai);
c_fft(datar,datai,fftdat[fftf_siz],1,fftdat[fftf_pow]);
scale = (float)(1.0/fftdat[fftf_siz]);
for (i=0;i<fftdat[fftf_siz];i++)
{
temp=datar[i];
datar[i]=scale*((datar[i]*filtr[i]) - (datai[i]*filti[i]));
datai[i]=scale*((temp*filti[i]) + (datai[i]*filtr[i]));
}
c_fft(datar,datai,fftdat[fftf_siz],-1,fftdat[fftf_pow]);
if (fftdat[fftf_hop] < 0)
{
fftdat[fftf_hop] = fftdat[fftf_half_siz];
fblk[blk_ctr] = (float)(fftdat[fftf_siz]/4);
readout_first_data(fftdat,runbuf,datar);
}
else
{
readout_data(fftdat,runbuf,datar);
fblk[blk_ctr] = (float)fftdat[fftf_hop];
}
}
return(c_run_block(runblk,fblk,runbuf));
}
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.