ftp.nice.ch/pub/next/unix/audio/Cmix.N.s.tar.gz#/cmix/rooms/place/placelib.c

This is placelib.c in view mode; [Download] [Up]

#define MACH1 1080.0
#define HPI 1.5707963268
#define IWRAP(x,y) (((x)>=(y))?(x)-(y):(((x)<0)?(x)+(y):(x)))
#define FWRAP(x,y) (((x)>=(y))?(x)-(y):(((x)<0.0)?(x)+(y):(x)))
#define FABS(x) (((x)<0.0)?-(x):(x))
#include <stdio.h>
#include <math.h>
#include "ugens.h"
#include "sfheader.h"
#include "firdata.h"
#include "debug.h"
#define RADPT 162.99746617261	/* converts rads to 1024-element array ptr */
#define RADPT2 325.9493234522
#define SIN(x) SINARRAY[(int)(wrap(x)*RADPT+0.5)]
#define COS(x) COSARRAY[(int)(wrap(x)*RADPT+0.5)]
#define ATAN(x) ATANARRAY[(int)((x)*RADPT2)+512]

extern SFHEADER sfdesc[NFILES];
extern int tapsize;
extern double SINARRAY[1025], COSARRAY[1025], ATANARRAY[1025];

/*
   This file contains the routines directly associated with the sound
   processor 'place'. This processor converts a mono source file into
   a stereo output file, in which the source will appear spacially lo-
   cated at any angle and distance from the listener, as well as simu-
   lating the acoustic response of a room of any size and (rectangular)
   dimensions. 
								        */

/* wrap converts negative or large polar angles (in rads) into positive */

double wrap(x)
register double x;
{
    register double val;
    val = x;
    while (val < 0.0)
	val += PI2;
    while (val >= PI2)
	val -= PI2;
    return (val);
}

/* this one converts all rads into vals -PI/2 <= val < PI/2 */

double wrap2(x)
register double x;
{
    register double val;
    val = x;
    while (val < -HPI)
	val += HPI;
    while (val >= HPI)
	val -= HPI;
    return (val);
}

/* binaural computes the stereo distance/angle pairs for roomtrig */

double binaural (R, T, X, Y, H, rho, theta)
    double R, T, X, Y, H;           /* distance, angle, coordinates, */
    double *rho, *theta;            /* "ear" dist., output arrays    */
{
    double h, Ys, HcosT;
    double asin(), sqrt(); 

    /* calculate distance pair */

    h  = H / 2.;
    Ys = Y * Y;
    rho[0] = sqrt ( Ys + ((X + h) * (X + h)) );
    rho[1] = sqrt ( Ys + ((X - h) * (X - h)) );

    /* calculate angle pair if H and T are not 0 */

    if (H != 0. && T != 0.)
    {
	HcosT = h * Y/R;
	theta[0] = wrap (-(T + asin(HcosT / rho[0])) );
	theta[1] = wrap (T - asin(HcosT / rho[1]) );
    }
    else
	theta[0] = theta[1] = T;
}

/* roomtrig calculates all distance/angle vectors for images up through the
   2nd generation reflections. 13 vectors total:  1 source, 4 1st gen., 8
   2nd gen. These are split into stereo pairs by binaural.               */


int roomtrig (A, B, H, rho, theta, cart)
    double A, B, H;	/* A is 'rho' or 'x'; B is 'theta' or 'y' */
    double rho[2][13], theta[2][13];
    int cart;
{
    register int i;
    double x[13], y[13], r[13], t[13], d[4], Ra[2], Ta[2];
    double X, Y, R, T;
    double sin(), cos(), hypot(), atan(), wrap();  
    extern double Dimensions[5];
    static double z = 0.017453292;    /* Pi / 180 */

    /* convert T into rads; calc. X & Y if entered in polar form only */

    if (!cart)		/* polar coordinates */
    {
	R = A;
        T = wrap(B*z);      /* convert into rads 0 <= T <= 2PI */ 
	X = A * SIN(T);
	Y = A * COS(T);
    }
    else
    {
        R = hypot(A, B);
	T = atan(A/B);
	X = A;
	Y = B;
    }
    /* Check to see that source loc. is inside room bounds */

    if (X < Dimensions[3] || X > Dimensions[1] || Y > Dimensions[0] ||
        Y < Dimensions[2])
    {
        fprintf (stderr,"Source loc outside room bounds!!\n\n");
	return (1);
    } 
    /* multiply global dimension array by 2 to save calc. time */

    for (i = 0; i < 4; ++i)
	d[i] = Dimensions[i] * 2.0;
    
    /* image calculations */

    /* source vector */

    x[0] = X;
    y[0] = Y;
    t[0] = T;
    r[0] = R;

    /* front wall vector */

    x[1] = X;
    y[1] = d[0] - Y;
    t[1] = atan( x[1]/y[1] );
    r[1] = y[1]/COS(t[1]); 

    /* right wall vector */

    x[2] = d[1] - X;
    y[2] = Y;
    t[2] = PI/2. - atan( y[2]/x[2] );
    r[2] = x[2] / SIN(t[2]);

    /* back wall vector */

    x[3] = X;
    y[3] = d[2] - Y;
    t[3] = PI + atan( x[3]/y[3] );
    r[3] = y[3] / COS(t[3]);

    /* left wall vector */

    x[4] = d[3] - X;
    y[4] = Y;
    t[4] = 3.* PI/2. - atan( y[4]/x[4] );
    r[4] = x[4] / SIN(t[4]);

    /* 2nd gen. images: 4 opposing wall, 4 adjacent wall reflections */

    /* front wall vector */

    x[5] = X;
    y[5] = d[0] - d[2] + Y;
    t[5] = atan( x[5]/y[5] );
    r[5] = hypot(X, y[5]);

    /* right wall vector */

    x[6] = d[1] - d[3] + X;
    y[6] = Y;
    t[6] = PI/2. - atan( y[6]/x[6] );
    r[6] = hypot(x[6], Y);

    /* back wall vector */

    x[7] = X;
    y[7] = d[2] - d[0] + Y;
    t[7] = PI + atan( x[7]/y[7] );
    r[7] = hypot(X, y[7]);

    /* left wall vector */

    x[8] = d[3] - d[1] + X;
    y[8] = Y;
    t[8] = 3.*PI/2. - atan( y[8]/x[8] );
    r[8] = hypot(x[8], Y);

    /* fr. rt. vector - double image in rectangular room, as are next 3 */

    x[9] = x[2];
    y[9] = y[1];
    t[9] = atan( x[9]/y[9] );
    r[9] = hypot(x[9], y[9]);

    /* back rt. vector */

    x[10] = x[2];
    y[10] = y[3];
    t[10] = PI/2. - atan( y[10]/x[10] );
    r[10] = hypot(x[10], y[10]);

    /* back lft. vector */

    x[11] = x[4];
    y[11] = y[3];
    t[11] = PI + atan( x[11]/y[11] );
    r[11] = hypot(x[11], y[11]);

    /* front lft. vector */

    x[12] = x[4];
    y[12] = y[1];
    t[12] = 3.*PI/2. - atan( y[12]/x[12] );
    r[12] = hypot(x[12], y[12]);

    /* calculate stereo vector pairs for each of these */

    for (i = 0; i < 13; ++i)
    {
	binaural(r[i], t[i], x[i], y[i], H, Ra, Ta);
	rho[0][i] = Ra[0];
	rho[1][i] = Ra[1];
	theta[0][i] = Ta[0];
	theta[1][i] = Ta[1];
    }

    /* Check to see that source distance is not "zero" */

    if(rho[0][0] < 0.001 || rho[1][0] < 0.001) 
    {
	fprintf(stderr, "Zero source distance not allowed!\n\n");
	return(1);
    }

    /* print out data for analysis */
#ifdef debug
    printf ("Source angles: %.2f    %.2f\n", theta[0][0]/z,theta[1][0]/z); 
    printf ("All others\n");
    for (i = 1; i < 13; i++)
	  printf ("%.2f     %.2f\n", theta[0][i]/z, theta[1][i]/z);
    printf ("Direct delays:\n");
    printf ("%d: %.2f ms.       %.2f ms.\n", i, rho[0][0]/1.08,rho[1][0]/1.08);
    printf ("Room delays:\n");
    for (i = 1; i < 13; i++)
	  printf ("%.2f ms.       %.2f ms.\n ", rho[0][i]/1.08, rho[1][i]/1.08);
    printf ("\n");   
    printf ("Direct dists:\n");
    printf ("%d: %.2f ft.       %.2f ft.\n", i, rho[0][0],rho[1][0]);
    printf ("Room dists:\n");
    for (i = 1; i < 13; i++)
	  printf ("%.2f ft.       %.2f ft.\n ", rho[0][i], rho[1][i]); 
    printf ("\n");  
#endif
    return (0); 
}


/* reset zeroes out all delay and filter histories in move and reverb   
   each time move is called */

rvb_reset (tapdel)
    double *tapdel;
{
    extern double Walldata[2][13][3], *Rvb_del[2][6];
    extern double Allpass_del[2][502], Rvb_air[2][6][3];
    extern int rvbdelsize;
    register int i, j, k;
    double *point;

    /* reset wall filter hists */

    for (i = 0; i < 2; ++i)
    {
	for (j = 0; j < 13; ++j)
	    Walldata[i][j][2] = 0.0;
    }
    /* reset reverb filters and delays */


    for (i = 0; i < 2; ++i)
    {
	for (j = 0, k = 0; j < 6; ++j)
	{
	    Rvb_air[i][j][2] = 0.0;
	    point = Rvb_del[i][j];   
	    while(k++ < rvbdelsize) *point++ = 0.0;
        } 
	for (j = 2; j < 502; ++j)
	    Allpass_del[i][j] = 0.0;
    }
    /* reset tap delay for move */

    for (i = 0; i < tapsize; ++i)
	*(tapdel+i) = 0.0;
}


/* setair calculates and loads the gain factors (G1 & G2) for the tone
   filters used to simulate air absorption. The values for G1 are stored
   in AIRCOEFFS by space(). G2 takes the 1/r-squared attenuation factor into
   account.  Histories are reset to 0 if flag = 1.                    */


setair (rho, flag, coeffs)
    double rho, *coeffs;
    int flag;
{
    extern float AIRCOEFFS[512];
    double G1, G2;
    float fpoint, frac;
    register int gpoint;

    /* pointer into array. Max distance for rho is 300 ft. */

    rho = (rho > 300 ? 300 : rho);
    fpoint = rho * 511/300.0;
    gpoint = fpoint;
    frac = fpoint - (float)gpoint;
    G1 = AIRCOEFFS[gpoint] + frac * (AIRCOEFFS[gpoint+1]-AIRCOEFFS[gpoint]);
    G2 = (1.0/rho) * (1.0 - G1);

    /* load into output array */

    coeffs[0] = G2;
    coeffs[1] = G1;

    /* reset if flag = 1 */

    if (flag)
	coeffs[2] = 0.0;
}


/* airfil_set is called by move to set all airfilters (2 for each image
   pair) whenever source moves.                                           */

airfil_set (rho, flag, airdata)
    double rho[2][13], airdata[2][13][3];
    int flag;
{
    register int i, j;

    for (i = 0; i < 2; ++i)
    {
	for (j = 0; j < 13; ++j)
	{
            setair (rho[i][j], flag, &(airdata[i][j][0]));
        }
    }
}


/* fir is a f.i.r. filter, set up by setfir, which filters each image
   sample pair to simulate binaural spectral cues according to angle */

double fir (in, counter, nterms, coeffs, firtap)
    double *coeffs, *firtap, in;
    long counter;
    int nterms;
{
    double out;
    register int intap, outtap, length, i;

    /* tap delay must be one longer than # of coeffs */

    length = nterms + 1;
    intap = counter % length;    /* input location */
    *(firtap+intap) = in;        /* load input sample into delay */

    /* for each of the delayed output samps, mult. by coeff and sum */

    out = 0.0;
    for (i = 1; i < length; ++i)
    {
	outtap = IWRAP((intap-i), length);
	out += *(firtap+outtap) * coeffs[i-1];
    }
    return (out);
}


/* setfir looks up the coeffs for the particular binaural angle, rho,
   and loads them into the firfilter memory. Also resets hist if needed. */

setfir (theta, nterms, flag, coeffs, firtap)
    double theta, *coeffs, *firtap;
    int nterms, flag;
{
    extern double FIRDATA[NANGLES][MAXTERMS]; /* the stored coeffs for firs */
    static double radmax = PI2;               /* 2PI rads */
    double rad_inc, angle, frac;
    register int lower, upper, skip, i, j;

    /* reset filter histories if flag = 1 */

    if (flag)
    {
	for (i = 0; i <= nterms; ++i)
	    firtap[i] = 0.0;
    }
    /* calculations to produce interpolated data */

    rad_inc = radmax / NANGLES;      /* distance in rads betw. data pts.  */
    angle = theta / rad_inc;         /* current angle in rad_incs */
    lower = angle;                   /* truncate to lower integer */
    frac = angle - lower;            /* for interpolating */
    upper = (lower + 1) % NANGLES;

    /* since not all firs use max # of terms stored, here is skip pt. */

    skip = (MAXTERMS - nterms) / 2;

    /* interpolate and load coefficients */

    for (i = 0; i < nterms; ++i)
    {   
	j = i + skip;
	coeffs[i] = FIRDATA[lower][j] + (FIRDATA[upper][j] -
		    FIRDATA[lower][j]) * frac;
    }
}


/* earfil_set is called by move to load coeffs for the fir filter bank
   used in ear, the binaural image filter.                            */

earfil_set (Theta, flag, Fircoeffs, Firtaps)
    double Theta[2][13], Fircoeffs[2][13][MAXTERMS];
    double Firtaps[2][13][MAXTERMS+1];
    int flag;
{
    register int i, j;
    extern int Nterms[13]; 

    for (i = 0; i < 2; ++i)
    {
	for (j = 0; j < 13; ++j)
	{
	    setfir (Theta[i][j], Nterms[j], flag, &(Fircoeffs[i][j][0]),
	            &(Firtaps[i][j][0]));
        }
    }
}


/* tap_set determines the number of samps delay for each of the tap delay
   lines in move.  It compensates for the group phase delay in ear, and
   determines the max # of samps delay for loop purposes, later in move. */

long tap_set (Rho, outlocs, fir_flag)
    double Rho[2][13];
    long outlocs[2][13];
    int fir_flag;
{
    register int i, j;
    long int maxloc;
    double delay;
    extern int Group_delay[13];       /* from Firdata.h */
    maxloc = 0;
#ifdef debug
    printf("outlocs:\n"); 
#endif
    for (i = 0; i < 2; ++i)
    {
	for (j = 0; j < 13; ++j)
	{
            delay = Rho[i][j] / MACH1;    /* delay time */
	    if (fir_flag)
		outlocs[i][j] = delay * SR - Group_delay[j] + 0.5;
            else
		outlocs[i][j] = delay * SR + 0.5;
            if (outlocs[i][j] > (float)maxloc)
		maxloc = outlocs[i][j] + 0.5;    /* max # of samps delay */
#ifdef debug
		printf("%ld ", outlocs[i][j]);
#endif
        }
#ifdef debug
	printf("\n");
#endif

    }
#ifdef debug
    printf("\n");
#endif
    return (maxloc);
}


/* get_tap accesses the tap delay array and fills the del. signal array, Sig. 


get_tap (intap, outlocs, Sig, tapdel)
    double Sig[2][13], *tapdel;
    int intap;
    long outlocs[2][13];
{
    register int i, j, outtap;
    
    for (i = 0; i < 2; ++i)
    {
	for (j = 0; j < 13; ++j)
	{
            outtap = IWRAP((intap - outlocs[i][j]), tapsize);
	    Sig[i][j] = tapdel[outtap];
        }
    }
} 
*/
/* get_tap accesses the tap delay array and fills the del. signal array, Sig.
*/

get_tap (intap, outlocs, Sig, tapdel)
    double Sig[2][13], *tapdel;
    int intap;
    long outlocs[2][13];
{
    register int i, j, outtap, outtap1;
    register float floattap, frac;

    for (i = 0; i < 2; ++i)
    {
	for (j = 0; j < 13; ++j)
	{
            floattap = FWRAP(((float)intap - outlocs[i][j]), (float)tapsize);
            outtap = floattap;
            frac = floattap - (float)outtap;
            outtap1 = outtap + 1;
            outtap1 = (outtap1>=tapsize)?outtap1-tapsize:outtap1;
	    Sig[i][j] = tapdel[outtap] + 
                        frac * (tapdel[outtap1] - tapdel[outtap]);
        }
    }
}  
/* ear is the bank of fir filters that simulate binaural angle-dependent
   spectral cues, in move.                                            */

ear (Sig, N, Fircoeffs, Firtaps)
    long int N;
    double Sig[2][13], Fircoeffs[2][13][MAXTERMS];
    double Firtaps[2][13][MAXTERMS + 1];
{
    register int i, j;
    extern int Nterms[13]; 

    for (i = 0; i < 2; ++i)
    {
	for (j = 0; j < 13; ++j)
	{
	    Sig[i][j] = fir (Sig[i][j], N, Nterms[j], &(Fircoeffs[i][j][0]),
	                     &(Firtaps[i][j][0]));
        }
    }
}


/* air is the bank of tone filters that simulate air absorption */

air (Sig, airdata)
    double Sig[2][13], airdata[2][13][3];
{
    register int i, j;
    double tone();

    for (i = 0; i < 2; ++i)
    {
	for (j = 0; j < 13; ++j)
            Sig[i][j] = tone (Sig[i][j], &(airdata[i][j][0]));
    }
}

/* wall is the bank of tone filters that simulate absorption of high
   frequencies by the walls of the room. These are set up in the separate
   setup routine, space.                                              */

wall (Sig)
    double Sig[2][13];
{
    extern double Walldata[2][13][3];
    double tone(); 
    register int i, j;

    for (i = 0; i < 2; ++i)
    {
	for (j = 1; j < 13; ++j)         /* no filt of 0 loc (dir sig.) */
            Sig[i][j] = tone (Sig[i][j], &(Walldata[i][j][0]));  
    }
}

These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.