This is movelib.c in view mode; [Download] [Up]
#define SIZE 512
#define MACH1 1080.0
#define HPI 1.57079632679490
#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 "nobug.h"
#define IPTR(a, ptr) (a[(int)ptr]+(ptr-(int)ptr)*(a[(int)ptr+1]-a[(int)ptr]))
#define RADPT 162.97466172610083 /* converts rads to 1024-element array ptr */
#define RADPT2 325.9493234522
#define ATAN(x) ATANARRAY[(int)((x)*RADPT2)+512]
extern int tapsize;
extern double ATANARRAY[1025];
/*
This file contains the routines directly associated with the sound
processor 'move'. 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.
*/
double
SIN(val)
double val;
{
extern double SINARRAY[1025];
double wrap();
double point = wrap(val) * RADPT;
point = point > 1023.0 ? point - 1023.0 : point;
return(IPTR(SINARRAY, point));
}
double
COS(val)
double val;
{
extern double COSARRAY[1025];
double wrap();
double point = wrap(val) * RADPT;
point = point > 1023.0 ? point - 1023.0 : point;
return(IPTR(COSARRAY, point));
}
/* 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, xmh, xph;
double asin(), sqrt();
extern int earflag;
/* calculate distance pair */
h = H / 2.;
Ys = Y * Y;
xph = X + h;
xmh = X - h;
rho[0] = sqrt ( Ys + (xph * xph) );
rho[1] = sqrt ( Ys + (xmh * xmh) );
/* calculate angle pair if earflag set */
if (earflag)
{
if(Y >= 0.0) {
theta[0] = PI2 - asin(xph / rho[0]);
theta[1] = asin(xmh / rho[1]);
}
else {
theta[0] = PI + asin(xph / rho[0]);
theta[1] = PI - asin(xmh / rho[1]);
}
theta[0] = wrap(theta[0]);
theta[1] = wrap(theta[1]);
}
else
theta[0] = theta[1] = 0.0;
}
/* 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)
double A, B, H; /* A is 'rho' or 'x'; B is 'theta' or 'y' */
double rho[2][13], theta[2][13];
{
register int i;
double x[13], y[13], r[13], t[13], d[4], Ra[2], Ta[2];
double X, Y, R, T, temp;
static double cnv = 360.0/PI2;
double SIN(), COS(), hypot(), atan(), wrap();
extern double Dimensions[5];
extern int earflag, cartflag;
/* convert T into rads; calc. X & Y if entered in polar form only */
if (!cartflag) /* polar coordinates */
{
R = (A==0.0) ? 0.0001 : A;
T = wrap(B); /* 0 <= T <= 2PI */
X = A * SIN(T);
Y = A * COS(T);
}
else
{
A = (A==0.0) ? 0.0001 : A;
B = (B==0.0) ? 0.0001 : B;
R = hypot(A, B);
T = atan(A/B);
X = A;
Y = B;
}
/*
printf("raw source loc: x = %f y = %f rho = %f theta = %f\n",X,Y,R,T*cnv);
*/
/* 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: distances (& angles for those where always needed) */
/* source vector */
x[0] = X;
y[0] = Y;
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] = HPI - 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.* HPI - 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;
r[5] = hypot(X, y[5]);
/* right wall vector */
x[6] = d[1] - d[3] + X;
y[6] = Y;
r[6] = hypot(x[6], Y);
/* back wall vector */
x[7] = X;
y[7] = d[2] - d[0] + Y;
r[7] = hypot(X, y[7]);
/* left wall vector */
x[8] = d[3] - d[1] + X;
y[8] = Y;
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];
r[9] = hypot(x[9], y[9]);
/* back rt. vector */
x[10] = x[2];
y[10] = y[3];
r[10] = hypot(x[10], y[10]);
/* back lft. vector */
x[11] = x[4];
y[11] = y[3];
r[11] = hypot(x[11], y[11]);
/* front lft. vector */
x[12] = x[4];
y[12] = y[1];
r[12] = hypot(x[12], y[12]);
/* image calculations: angles for remaining (if earflag set) */
if(earflag) {
/* source vector */
t[0] = T;
r[0] = R;
/* 2nd gen. images: 4 opposing wall, 4 adjacent wall reflections */
t[5] = atan( x[5]/y[5] ); /* front wall vector */
t[6] = HPI - atan( y[6]/x[6] ); /* right wall vector */
t[7] = PI + atan( x[7]/y[7] ); /* back wall vector */
t[8] = 3.* HPI - atan( y[8]/x[8] ); /* left wall vector */
/* fr. rt. vector - dbl. image in rectangular room, as are next 3 */
t[9] = atan( x[9]/y[9] );
t[10] = HPI - atan( y[10]/x[10] ); /* back rt. vector */
t[11] = PI + atan( x[11]/y[11] ); /* back lft. vector */
t[12] = 3.* HPI - atan( y[12]/x[12] ); /* front lft. vector */
}
/* 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 */
/*
printf ("Source angles: %.2f %.2f\n", theta[0][0]*cnv,theta[1][0]*cnv);
printf ("All others\n");
for (i = 1; i < 13; i++)
printf ("%.2f %.2f\n", theta[0][i]*cnv, theta[1][i]*cnv);
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");
*/
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]);
}
}
}
/* 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],
Firtaps[i][j]);
}
}
}
/* 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];
double 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;
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];
else
outlocs[i][j] = delay * SR;
if (outlocs[i][j] > (float)maxloc)
maxloc = outlocs[i][j] + 0.5; /* max # of samps delay */
}
}
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;
double 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;
/* printf("floattap: %f outtap: %d outtap1: %d frac: %f\n", floattap, outtap, outtap1, frac); */
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],
Firtaps[i][j]);
}
}
}
/* 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]);
}
}
/* 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]);
}
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.