AuroraMiddleware/AuroraOpenALSoft
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
603e3e23ac
AuroraOpenALSoft/utils/makehrtf.c
Chris Robinson 57e516720e Avoid a potential calloc of 0
2017-12-17 14:30:51 -08:00

3475 lines
104 KiB
C
Raw Blame History

/*
* HRTF utility for producing and demonstrating the process of creating an
* OpenAL Soft compatible HRIR data set.
*
* Copyright (C) 2011-2017 Christopher Fitzgerald
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along
* with this program; if not, write to the Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
* Or visit: http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
*
* --------------------------------------------------------------------------
*
* A big thanks goes out to all those whose work done in the field of
* binaural sound synthesis using measured HRTFs makes this utility and the
* OpenAL Soft implementation possible.
*
* The algorithm for diffuse-field equalization was adapted from the work
* done by Rio Emmanuel and Larcher Veronique of IRCAM and Bill Gardner of
* MIT Media Laboratory. It operates as follows:
*
* 1. Take the FFT of each HRIR and only keep the magnitude responses.
* 2. Calculate the diffuse-field power-average of all HRIRs weighted by
* their contribution to the total surface area covered by their
* measurement.
* 3. Take the diffuse-field average and limit its magnitude range.
* 4. Equalize the responses by using the inverse of the diffuse-field
* average.
* 5. Reconstruct the minimum-phase responses.
* 5. Zero the DC component.
* 6. IFFT the result and truncate to the desired-length minimum-phase FIR.
*
* The spherical head algorithm for calculating propagation delay was adapted
* from the paper:
*
* Modeling Interaural Time Difference Assuming a Spherical Head
* Joel David Miller
* Music 150, Musical Acoustics, Stanford University
* December 2, 2001
*
* The formulae for calculating the Kaiser window metrics are from the
* the textbook:
*
* Discrete-Time Signal Processing
* Alan V. Oppenheim and Ronald W. Schafer
* Prentice-Hall Signal Processing Series
* 1999
*/
#include "config.h"
#define _UNICODE
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <stddef.h>
#include <string.h>
#include <limits.h>
#include <ctype.h>
#include <math.h>
#ifdef HAVE_STRINGS_H
#include <strings.h>
#endif
#ifdef HAVE_GETOPT
#include <unistd.h>
#else
#include "getopt.h"
#endif
#include "win_main_utf8.h"
/* Define int64_t and uint64_t types */
#if defined(__STDC_VERSION__) && __STDC_VERSION__ >= 199901L
#include <inttypes.h>
#elif defined(_WIN32) && defined(__GNUC__)
#include <stdint.h>
#elif defined(_WIN32)
typedef __int64 int64_t;
typedef unsigned __int64 uint64_t;
#else
/* Fallback if nothing above works */
#include <inttypes.h>
#endif
#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif
#ifndef HUGE_VAL
#define HUGE_VAL (1.0 / 0.0)
#endif
// The epsilon used to maintain signal stability.
#define EPSILON (1e-9)
// Constants for accessing the token reader's ring buffer.
#define TR_RING_BITS (16)
#define TR_RING_SIZE (1 << TR_RING_BITS)
#define TR_RING_MASK (TR_RING_SIZE - 1)
// The token reader's load interval in bytes.
#define TR_LOAD_SIZE (TR_RING_SIZE >> 2)
// The maximum identifier length used when processing the data set
// definition.
#define MAX_IDENT_LEN (16)
// The maximum path length used when processing filenames.
#define MAX_PATH_LEN (256)
// The limits for the sample 'rate' metric in the data set definition and for
// resampling.
#define MIN_RATE (32000)
#define MAX_RATE (96000)
// The limits for the HRIR 'points' metric in the data set definition.
#define MIN_POINTS (16)
#define MAX_POINTS (8192)
// The limit to the number of 'distances' listed in the data set definition.
#define MAX_FD_COUNT (16)
// The limits to the number of 'azimuths' listed in the data set definition.
#define MIN_EV_COUNT (5)
#define MAX_EV_COUNT (128)
// The limits for each of the 'azimuths' listed in the data set definition.
#define MIN_AZ_COUNT (1)
#define MAX_AZ_COUNT (128)
// The limits for the listener's head 'radius' in the data set definition.
#define MIN_RADIUS (0.05)
#define MAX_RADIUS (0.15)
// The limits for the 'distance' from source to listener for each field in
// the definition file.
#define MIN_DISTANCE (0.05)
#define MAX_DISTANCE (2.50)
// The maximum number of channels that can be addressed for a WAVE file
// source listed in the data set definition.
#define MAX_WAVE_CHANNELS (65535)
// The limits to the byte size for a binary source listed in the definition
// file.
#define MIN_BIN_SIZE (2)
#define MAX_BIN_SIZE (4)
// The minimum number of significant bits for binary sources listed in the
// data set definition. The maximum is calculated from the byte size.
#define MIN_BIN_BITS (16)
// The limits to the number of significant bits for an ASCII source listed in
// the data set definition.
#define MIN_ASCII_BITS (16)
#define MAX_ASCII_BITS (32)
// The limits to the FFT window size override on the command line.
#define MIN_FFTSIZE (65536)
#define MAX_FFTSIZE (131072)
// The limits to the equalization range limit on the command line.
#define MIN_LIMIT (2.0)
#define MAX_LIMIT (120.0)
// The limits to the truncation window size on the command line.
#define MIN_TRUNCSIZE (16)
#define MAX_TRUNCSIZE (512)
// The limits to the custom head radius on the command line.
#define MIN_CUSTOM_RADIUS (0.05)
#define MAX_CUSTOM_RADIUS (0.15)
// The truncation window size must be a multiple of the below value to allow
// for vectorized convolution.
#define MOD_TRUNCSIZE (8)
// The defaults for the command line options.
#define DEFAULT_FFTSIZE (65536)
#define DEFAULT_EQUALIZE (1)
#define DEFAULT_SURFACE (1)
#define DEFAULT_LIMIT (24.0)
#define DEFAULT_TRUNCSIZE (32)
#define DEFAULT_HEAD_MODEL (HM_DATASET)
#define DEFAULT_CUSTOM_RADIUS (0.0)
// The four-character-codes for RIFF/RIFX WAVE file chunks.
#define FOURCC_RIFF (0x46464952) // 'RIFF'
#define FOURCC_RIFX (0x58464952) // 'RIFX'
#define FOURCC_WAVE (0x45564157) // 'WAVE'
#define FOURCC_FMT (0x20746D66) // 'fmt '
#define FOURCC_DATA (0x61746164) // 'data'
#define FOURCC_LIST (0x5453494C) // 'LIST'
#define FOURCC_WAVL (0x6C766177) // 'wavl'
#define FOURCC_SLNT (0x746E6C73) // 'slnt'
// The supported wave formats.
#define WAVE_FORMAT_PCM (0x0001)
#define WAVE_FORMAT_IEEE_FLOAT (0x0003)
#define WAVE_FORMAT_EXTENSIBLE (0xFFFE)
// The maximum propagation delay value supported by OpenAL Soft.
#define MAX_HRTD (63.0)
// The OpenAL Soft HRTF format marker. It stands for minimum-phase head
// response protocol 02.
#define MHR_FORMAT ("MinPHR02")
// Sample and channel type enum values.
typedef enum SampleTypeT {
ST_S16 = 0,
ST_S24 = 1
} SampleTypeT;
// Certain iterations rely on these integer enum values.
typedef enum ChannelTypeT {
CT_NONE = -1,
CT_MONO = 0,
CT_STEREO = 1
} ChannelTypeT;
// Byte order for the serialization routines.
typedef enum ByteOrderT {
BO_NONE,
BO_LITTLE,
BO_BIG
} ByteOrderT;
// Source format for the references listed in the data set definition.
typedef enum SourceFormatT {
SF_NONE,
SF_WAVE, // RIFF/RIFX WAVE file.
SF_BIN_LE, // Little-endian binary file.
SF_BIN_BE, // Big-endian binary file.
SF_ASCII // ASCII text file.
} SourceFormatT;
// Element types for the references listed in the data set definition.
typedef enum ElementTypeT {
ET_NONE,
ET_INT, // Integer elements.
ET_FP // Floating-point elements.
} ElementTypeT;
// Head model used for calculating the impulse delays.
typedef enum HeadModelT {
HM_NONE,
HM_DATASET, // Measure the onset from the dataset.
HM_SPHERE // Calculate the onset using a spherical head model.
} HeadModelT;
// Unsigned integer type.
typedef unsigned int uint;
// Serialization types. The trailing digit indicates the number of bits.
typedef unsigned char uint8;
typedef int int32;
typedef unsigned int uint32;
typedef uint64_t uint64;
// Token reader state for parsing the data set definition.
typedef struct TokenReaderT {
FILE *mFile;
const char *mName;
uint mLine;
uint mColumn;
char mRing[TR_RING_SIZE];
size_t mIn;
size_t mOut;
} TokenReaderT;
// Source reference state used when loading sources.
typedef struct SourceRefT {
SourceFormatT mFormat;
ElementTypeT mType;
uint mSize;
int mBits;
uint mChannel;
uint mSkip;
uint mOffset;
char mPath[MAX_PATH_LEN+1];
} SourceRefT;
// Structured HRIR storage for stereo azimuth pairs, elevations, and fields.
typedef struct HrirAzT {
double mAzimuth;
uint mIndex;
double mDelays[2];
double *mIrs[2];
} HrirAzT;
typedef struct HrirEvT {
double mElevation;
uint mIrCount;
uint mAzCount;
HrirAzT *mAzs;
} HrirEvT;
typedef struct HrirFdT {
double mDistance;
uint mIrCount;
uint mEvCount;
uint mEvStart;
HrirEvT *mEvs;
} HrirFdT;
// The HRIR metrics and data set used when loading, processing, and storing
// the resulting HRTF.
typedef struct HrirDataT {
uint mIrRate;
SampleTypeT mSampleType;
ChannelTypeT mChannelType;
uint mIrPoints;
uint mFftSize;
uint mIrSize;
double mRadius;
uint mIrCount;
uint mFdCount;
HrirFdT *mFds;
} HrirDataT;
// The resampler metrics and FIR filter.
typedef struct ResamplerT {
uint mP, mQ, mM, mL;
double *mF;
} ResamplerT;
/****************************************
*** Complex number type and routines ***
****************************************/
typedef struct {
double Real, Imag;
} Complex;
static Complex MakeComplex(double r, double i)
{
Complex c = { r, i };
return c;
}
static Complex c_add(Complex a, Complex b)
{
Complex r;
r.Real = a.Real + b.Real;
r.Imag = a.Imag + b.Imag;
return r;
}
static Complex c_sub(Complex a, Complex b)
{
Complex r;
r.Real = a.Real - b.Real;
r.Imag = a.Imag - b.Imag;
return r;
}
static Complex c_mul(Complex a, Complex b)
{
Complex r;
r.Real = a.Real*b.Real - a.Imag*b.Imag;
r.Imag = a.Imag*b.Real + a.Real*b.Imag;
return r;
}
static Complex c_muls(Complex a, double s)
{
Complex r;
r.Real = a.Real * s;
r.Imag = a.Imag * s;
return r;
}
static double c_abs(Complex a)
{
return sqrt(a.Real*a.Real + a.Imag*a.Imag);
}
static Complex c_exp(Complex a)
{
Complex r;
double e = exp(a.Real);
r.Real = e * cos(a.Imag);
r.Imag = e * sin(a.Imag);
return r;
}
/*****************************
*** Token reader routines ***
*****************************/
/* Whitespace is not significant. It can process tokens as identifiers, numbers
* (integer and floating-point), strings, and operators. Strings must be
* encapsulated by double-quotes and cannot span multiple lines.
*/
// Setup the reader on the given file. The filename can be NULL if no error
// output is desired.
static void TrSetup(FILE *fp, const char *filename, TokenReaderT *tr)
{
const char *name = NULL;
if(filename)
{
const char *slash = strrchr(filename, '/');
if(slash)
{
const char *bslash = strrchr(slash+1, '\\');
if(bslash) name = bslash+1;
else name = slash+1;
}
else
{
const char *bslash = strrchr(filename, '\\');
if(bslash) name = bslash+1;
else name = filename;
}
}
tr->mFile = fp;
tr->mName = name;
tr->mLine = 1;
tr->mColumn = 1;
tr->mIn = 0;
tr->mOut = 0;
}
// Prime the reader's ring buffer, and return a result indicating that there
// is text to process.
static int TrLoad(TokenReaderT *tr)
{
size_t toLoad, in, count;
toLoad = TR_RING_SIZE - (tr->mIn - tr->mOut);
if(toLoad >= TR_LOAD_SIZE && !feof(tr->mFile))
{
// Load TR_LOAD_SIZE (or less if at the end of the file) per read.
toLoad = TR_LOAD_SIZE;
in = tr->mIn&TR_RING_MASK;
count = TR_RING_SIZE - in;
if(count < toLoad)
{
tr->mIn += fread(&tr->mRing[in], 1, count, tr->mFile);
tr->mIn += fread(&tr->mRing[0], 1, toLoad-count, tr->mFile);
}
else
tr->mIn += fread(&tr->mRing[in], 1, toLoad, tr->mFile);
if(tr->mOut >= TR_RING_SIZE)
{
tr->mOut -= TR_RING_SIZE;
tr->mIn -= TR_RING_SIZE;
}
}
if(tr->mIn > tr->mOut)
return 1;
return 0;
}
// Error display routine. Only displays when the base name is not NULL.
static void TrErrorVA(const TokenReaderT *tr, uint line, uint column, const char *format, va_list argPtr)
{
if(!tr->mName)
return;
fprintf(stderr, "Error (%s:%u:%u): ", tr->mName, line, column);
vfprintf(stderr, format, argPtr);
}
// Used to display an error at a saved line/column.
static void TrErrorAt(const TokenReaderT *tr, uint line, uint column, const char *format, ...)
{
va_list argPtr;
va_start(argPtr, format);
TrErrorVA(tr, line, column, format, argPtr);
va_end(argPtr);
}
// Used to display an error at the current line/column.
static void TrError(const TokenReaderT *tr, const char *format, ...)
{
va_list argPtr;
va_start(argPtr, format);
TrErrorVA(tr, tr->mLine, tr->mColumn, format, argPtr);
va_end(argPtr);
}
// Skips to the next line.
static void TrSkipLine(TokenReaderT *tr)
{
char ch;
while(TrLoad(tr))
{
ch = tr->mRing[tr->mOut&TR_RING_MASK];
tr->mOut++;
if(ch == '\n')
{
tr->mLine++;
tr->mColumn = 1;
break;
}
tr->mColumn ++;
}
}
// Skips to the next token.
static int TrSkipWhitespace(TokenReaderT *tr)
{
char ch;
while(TrLoad(tr))
{
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(isspace(ch))
{
tr->mOut++;
if(ch == '\n')
{
tr->mLine++;
tr->mColumn = 1;
}
else
tr->mColumn++;
}
else if(ch == '#')
TrSkipLine(tr);
else
return 1;
}
return 0;
}
// Get the line and/or column of the next token (or the end of input).
static void TrIndication(TokenReaderT *tr, uint *line, uint *column)
{
TrSkipWhitespace(tr);
if(line) *line = tr->mLine;
if(column) *column = tr->mColumn;
}
// Checks to see if a token is (likely to be) an identifier. It does not
// display any errors and will not proceed to the next token.
static int TrIsIdent(TokenReaderT *tr)
{
char ch;
if(!TrSkipWhitespace(tr))
return 0;
ch = tr->mRing[tr->mOut&TR_RING_MASK];
return ch == '_' || isalpha(ch);
}
// Checks to see if a token is the given operator. It does not display any
// errors and will not proceed to the next token.
static int TrIsOperator(TokenReaderT *tr, const char *op)
{
size_t out, len;
char ch;
if(!TrSkipWhitespace(tr))
return 0;
out = tr->mOut;
len = 0;
while(op[len] != '\0' && out < tr->mIn)
{
ch = tr->mRing[out&TR_RING_MASK];
if(ch != op[len]) break;
len++;
out++;
}
if(op[len] == '\0')
return 1;
return 0;
}
/* The TrRead*() routines obtain the value of a matching token type. They
* display type, form, and boundary errors and will proceed to the next
* token.
*/
// Reads and validates an identifier token.
static int TrReadIdent(TokenReaderT *tr, const uint maxLen, char *ident)
{
uint col, len;
char ch;
col = tr->mColumn;
if(TrSkipWhitespace(tr))
{
col = tr->mColumn;
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(ch == '_' || isalpha(ch))
{
len = 0;
do {
if(len < maxLen)
ident[len] = ch;
len++;
tr->mOut++;
if(!TrLoad(tr))
break;
ch = tr->mRing[tr->mOut&TR_RING_MASK];
} while(ch == '_' || isdigit(ch) || isalpha(ch));
tr->mColumn += len;
if(len < maxLen)
{
ident[len] = '\0';
return 1;
}
TrErrorAt(tr, tr->mLine, col, "Identifier is too long.\n");
return 0;
}
}
TrErrorAt(tr, tr->mLine, col, "Expected an identifier.\n");
return 0;
}
// Reads and validates (including bounds) an integer token.
static int TrReadInt(TokenReaderT *tr, const int loBound, const int hiBound, int *value)
{
uint col, digis, len;
char ch, temp[64+1];
col = tr->mColumn;
if(TrSkipWhitespace(tr))
{
col = tr->mColumn;
len = 0;
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(ch == '+' || ch == '-')
{
temp[len] = ch;
len++;
tr->mOut++;
}
digis = 0;
while(TrLoad(tr))
{
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(!isdigit(ch)) break;
if(len < 64)
temp[len] = ch;
len++;
digis++;
tr->mOut++;
}
tr->mColumn += len;
if(digis > 0 && ch != '.' && !isalpha(ch))
{
if(len > 64)
{
TrErrorAt(tr, tr->mLine, col, "Integer is too long.");
return 0;
}
temp[len] = '\0';
*value = strtol(temp, NULL, 10);
if(*value < loBound || *value > hiBound)
{
TrErrorAt(tr, tr->mLine, col, "Expected a value from %d to %d.\n", loBound, hiBound);
return 0;
}
return 1;
}
}
TrErrorAt(tr, tr->mLine, col, "Expected an integer.\n");
return 0;
}
// Reads and validates (including bounds) a float token.
static int TrReadFloat(TokenReaderT *tr, const double loBound, const double hiBound, double *value)
{
uint col, digis, len;
char ch, temp[64+1];
col = tr->mColumn;
if(TrSkipWhitespace(tr))
{
col = tr->mColumn;
len = 0;
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(ch == '+' || ch == '-')
{
temp[len] = ch;
len++;
tr->mOut++;
}
digis = 0;
while(TrLoad(tr))
{
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(!isdigit(ch)) break;
if(len < 64)
temp[len] = ch;
len++;
digis++;
tr->mOut++;
}
if(ch == '.')
{
if(len < 64)
temp[len] = ch;
len++;
tr->mOut++;
}
while(TrLoad(tr))
{
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(!isdigit(ch)) break;
if(len < 64)
temp[len] = ch;
len++;
digis++;
tr->mOut++;
}
if(digis > 0)
{
if(ch == 'E' || ch == 'e')
{
if(len < 64)
temp[len] = ch;
len++;
digis = 0;
tr->mOut++;
if(ch == '+' || ch == '-')
{
if(len < 64)
temp[len] = ch;
len++;
tr->mOut++;
}
while(TrLoad(tr))
{
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(!isdigit(ch)) break;
if(len < 64)
temp[len] = ch;
len++;
digis++;
tr->mOut++;
}
}
tr->mColumn += len;
if(digis > 0 && ch != '.' && !isalpha(ch))
{
if(len > 64)
{
TrErrorAt(tr, tr->mLine, col, "Float is too long.");
return 0;
}
temp[len] = '\0';
*value = strtod(temp, NULL);
if(*value < loBound || *value > hiBound)
{
TrErrorAt(tr, tr->mLine, col, "Expected a value from %f to %f.\n", loBound, hiBound);
return 0;
}
return 1;
}
}
else
tr->mColumn += len;
}
TrErrorAt(tr, tr->mLine, col, "Expected a float.\n");
return 0;
}
// Reads and validates a string token.
static int TrReadString(TokenReaderT *tr, const uint maxLen, char *text)
{
uint col, len;
char ch;
col = tr->mColumn;
if(TrSkipWhitespace(tr))
{
col = tr->mColumn;
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(ch == '\"')
{
tr->mOut++;
len = 0;
while(TrLoad(tr))
{
ch = tr->mRing[tr->mOut&TR_RING_MASK];
tr->mOut++;
if(ch == '\"')
break;
if(ch == '\n')
{
TrErrorAt(tr, tr->mLine, col, "Unterminated string at end of line.\n");
return 0;
}
if(len < maxLen)
text[len] = ch;
len++;
}
if(ch != '\"')
{
tr->mColumn += 1 + len;
TrErrorAt(tr, tr->mLine, col, "Unterminated string at end of input.\n");
return 0;
}
tr->mColumn += 2 + len;
if(len > maxLen)
{
TrErrorAt(tr, tr->mLine, col, "String is too long.\n");
return 0;
}
text[len] = '\0';
return 1;
}
}
TrErrorAt(tr, tr->mLine, col, "Expected a string.\n");
return 0;
}
// Reads and validates the given operator.
static int TrReadOperator(TokenReaderT *tr, const char *op)
{
uint col, len;
char ch;
col = tr->mColumn;
if(TrSkipWhitespace(tr))
{
col = tr->mColumn;
len = 0;
while(op[len] != '\0' && TrLoad(tr))
{
ch = tr->mRing[tr->mOut&TR_RING_MASK];
if(ch != op[len]) break;
len++;
tr->mOut++;
}
tr->mColumn += len;
if(op[len] == '\0')
return 1;
}
TrErrorAt(tr, tr->mLine, col, "Expected '%s' operator.\n", op);
return 0;
}
/* Performs a string substitution. Any case-insensitive occurrences of the
* pattern string are replaced with the replacement string. The result is
* truncated if necessary.
*/
static int StrSubst(const char *in, const char *pat, const char *rep, const size_t maxLen, char *out)
{
size_t inLen, patLen, repLen;
size_t si, di;
int truncated;
inLen = strlen(in);
patLen = strlen(pat);
repLen = strlen(rep);
si = 0;
di = 0;
truncated = 0;
while(si < inLen && di < maxLen)
{
if(patLen <= inLen-si)
{
if(strncasecmp(&in[si], pat, patLen) == 0)
{
if(repLen > maxLen-di)
{
repLen = maxLen - di;
truncated = 1;
}
strncpy(&out[di], rep, repLen);
si += patLen;
di += repLen;
}
}
out[di] = in[si];
si++;
di++;
}
if(si < inLen)
truncated = 1;
out[di] = '\0';
return !truncated;
}
/*********************
*** Math routines ***
*********************/
// Provide missing math routines for MSVC versions < 1800 (Visual Studio 2013).
#if defined(_MSC_VER) && _MSC_VER < 1800
static double round(double val)
{
if(val < 0.0)
return ceil(val-0.5);
return floor(val+0.5);
}
static double fmin(double a, double b)
{
return (a<b) ? a : b;
}
static double fmax(double a, double b)
{
return (a>b) ? a : b;
}
#endif
// Simple clamp routine.
static double Clamp(const double val, const double lower, const double upper)
{
return fmin(fmax(val, lower), upper);
}
// Performs linear interpolation.
static double Lerp(const double a, const double b, const double f)
{
return a + f * (b - a);
}
static inline uint dither_rng(uint *seed)
{
*seed = *seed * 96314165 + 907633515;
return *seed;
}
// Performs a triangular probability density function dither. The input samples
// should be normalized (-1 to +1).
static void TpdfDither(double *restrict out, const double *restrict in, const double scale,
const int count, const int step, uint *seed)
{
static const double PRNG_SCALE = 1.0 / UINT_MAX;
uint prn0, prn1;
int i;
for(i = 0;i < count;i++)
{
prn0 = dither_rng(seed);
prn1 = dither_rng(seed);
out[i*step] = round(in[i]*scale + (prn0*PRNG_SCALE - prn1*PRNG_SCALE));
}
}
// Allocates an array of doubles.
static double *CreateDoubles(size_t n)
{
double *a;
a = calloc(n?n:1, sizeof(*a));
if(a == NULL)
{
fprintf(stderr, "Error: Out of memory.\n");
exit(-1);
}
return a;
}
// Allocates an array of complex numbers.
static Complex *CreateComplexes(size_t n)
{
Complex *a;
a = calloc(n?n:1, sizeof(*a));
if(a == NULL)
{
fprintf(stderr, "Error: Out of memory.\n");
exit(-1);
}
return a;
}
/* Fast Fourier transform routines. The number of points must be a power of
* two. In-place operation is possible only if both the real and imaginary
* parts are in-place together.
*/
// Performs bit-reversal ordering.
static void FftArrange(const uint n, const Complex *in, Complex *out)
{
uint rk, k, m;
if(in == out)
{
// Handle in-place arrangement.
rk = 0;
for(k = 0;k < n;k++)
{
if(rk > k)
{
Complex temp = in[rk];
out[rk] = in[k];
out[k] = temp;
}
m = n;
while(rk&(m >>= 1))
rk &= ~m;
rk |= m;
}
}
else
{
// Handle copy arrangement.
rk = 0;
for(k = 0;k < n;k++)
{
out[rk] = in[k];
m = n;
while(rk&(m >>= 1))
rk &= ~m;
rk |= m;
}
}
}
// Performs the summation.
static void FftSummation(const int n, const double s, Complex *cplx)
{
double pi;
int m, m2;
int i, k, mk;
pi = s * M_PI;
for(m = 1, m2 = 2;m < n; m <<= 1, m2 <<= 1)
{
// v = Complex (-2.0 * sin (0.5 * pi / m) * sin (0.5 * pi / m), -sin (pi / m))
double sm = sin(0.5 * pi / m);
Complex v = MakeComplex(-2.0*sm*sm, -sin(pi / m));
Complex w = MakeComplex(1.0, 0.0);
for(i = 0;i < m;i++)
{
for(k = i;k < n;k += m2)
{
Complex t;
mk = k + m;
t = c_mul(w, cplx[mk]);
cplx[mk] = c_sub(cplx[k], t);
cplx[k] = c_add(cplx[k], t);
}
w = c_add(w, c_mul(v, w));
}
}
}
// Performs a forward FFT.
static void FftForward(const uint n, const Complex *in, Complex *out)
{
FftArrange(n, in, out);
FftSummation(n, 1.0, out);
}
// Performs an inverse FFT.
static void FftInverse(const uint n, const Complex *in, Complex *out)
{
double f;
uint i;
FftArrange(n, in, out);
FftSummation(n, -1.0, out);
f = 1.0 / n;
for(i = 0;i < n;i++)
out[i] = c_muls(out[i], f);
}
/* Calculate the complex helical sequence (or discrete-time analytical signal)
* of the given input using the Hilbert transform. Given the natural logarithm
* of a signal's magnitude response, the imaginary components can be used as
* the angles for minimum-phase reconstruction.
*/
static void Hilbert(const uint n, const Complex *in, Complex *out)
{
uint i;
if(in == out)
{
// Handle in-place operation.
for(i = 0;i < n;i++)
out[i].Imag = 0.0;
}
else
{
// Handle copy operation.
for(i = 0;i < n;i++)
out[i] = MakeComplex(in[i].Real, 0.0);
}
FftInverse(n, out, out);
for(i = 1;i < (n+1)/2;i++)
out[i] = c_muls(out[i], 2.0);
/* Increment i if n is even. */
i += (n&1)^1;
for(;i < n;i++)
out[i] = MakeComplex(0.0, 0.0);
FftForward(n, out, out);
}
/* Calculate the magnitude response of the given input. This is used in
* place of phase decomposition, since the phase residuals are discarded for
* minimum phase reconstruction. The mirrored half of the response is also
* discarded.
*/
static void MagnitudeResponse(const uint n, const Complex *in, double *out)
{
const uint m = 1 + (n / 2);
uint i;
for(i = 0;i < m;i++)
out[i] = fmax(c_abs(in[i]), EPSILON);
}
/* Apply a range limit (in dB) to the given magnitude response. This is used
* to adjust the effects of the diffuse-field average on the equalization
* process.
*/
static void LimitMagnitudeResponse(const uint n, const uint m, const double limit, const double *in, double *out)
{
double halfLim;
uint i, lower, upper;
double ave;
halfLim = limit / 2.0;
// Convert the response to dB.
for(i = 0;i < m;i++)
out[i] = 20.0 * log10(in[i]);
// Use six octaves to calculate the average magnitude of the signal.
lower = ((uint)ceil(n / pow(2.0, 8.0))) - 1;
upper = ((uint)floor(n / pow(2.0, 2.0))) - 1;
ave = 0.0;
for(i = lower;i <= upper;i++)
ave += out[i];
ave /= upper - lower + 1;
// Keep the response within range of the average magnitude.
for(i = 0;i < m;i++)
out[i] = Clamp(out[i], ave - halfLim, ave + halfLim);
// Convert the response back to linear magnitude.
for(i = 0;i < m;i++)
out[i] = pow(10.0, out[i] / 20.0);
}
/* Reconstructs the minimum-phase component for the given magnitude response
* of a signal. This is equivalent to phase recomposition, sans the missing
* residuals (which were discarded). The mirrored half of the response is
* reconstructed.
*/
static void MinimumPhase(const uint n, const double *in, Complex *out)
{
const uint m = 1 + (n / 2);
double *mags;
uint i;
mags = CreateDoubles(n);
for(i = 0;i < m;i++)
{
mags[i] = fmax(EPSILON, in[i]);
out[i] = MakeComplex(log(mags[i]), 0.0);
}
for(;i < n;i++)
{
mags[i] = mags[n - i];
out[i] = out[n - i];
}
Hilbert(n, out, out);
// Remove any DC offset the filter has.
mags[0] = EPSILON;
for(i = 0;i < n;i++)
{
Complex a = c_exp(MakeComplex(0.0, out[i].Imag));
out[i] = c_mul(MakeComplex(mags[i], 0.0), a);
}
free(mags);
}
/***************************
*** Resampler functions ***
***************************/
/* This is the normalized cardinal sine (sinc) function.
*
* sinc(x) = { 1, x = 0
* { sin(pi x) / (pi x), otherwise.
*/
static double Sinc(const double x)
{
if(fabs(x) < EPSILON)
return 1.0;
return sin(M_PI * x) / (M_PI * x);
}
/* The zero-order modified Bessel function of the first kind, used for the
* Kaiser window.
*
* I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
* = sum_{k=0}^inf ((x / 2)^k / k!)^2
*/
static double BesselI_0(const double x)
{
double term, sum, x2, y, last_sum;
int k;
// Start at k=1 since k=0 is trivial.
term = 1.0;
sum = 1.0;
x2 = x/2.0;
k = 1;
// Let the integration converge until the term of the sum is no longer
// significant.
do {
y = x2 / k;
k++;
last_sum = sum;
term *= y * y;
sum += term;
} while(sum != last_sum);
return sum;
}
/* Calculate a Kaiser window from the given beta value and a normalized k
* [-1, 1].
*
* w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B), -1 <= k <= 1
* { 0, elsewhere.
*
* Where k can be calculated as:
*
* k = i / l, where -l <= i <= l.
*
* or:
*
* k = 2 i / M - 1, where 0 <= i <= M.
*/
static double Kaiser(const double b, const double k)
{
if(!(k >= -1.0 && k <= 1.0))
return 0.0;
return BesselI_0(b * sqrt(1.0 - k*k)) / BesselI_0(b);
}
// Calculates the greatest common divisor of a and b.
static uint Gcd(uint x, uint y)
{
while(y > 0)
{
uint z = y;
y = x % y;
x = z;
}
return x;
}
/* Calculates the size (order) of the Kaiser window. Rejection is in dB and
* the transition width is normalized frequency (0.5 is nyquist).
*
* M = { ceil((r - 7.95) / (2.285 2 pi f_t)), r > 21
* { ceil(5.79 / 2 pi f_t), r <= 21.
*
*/
static uint CalcKaiserOrder(const double rejection, const double transition)
{
double w_t = 2.0 * M_PI * transition;
if(rejection > 21.0)
return (uint)ceil((rejection - 7.95) / (2.285 * w_t));
return (uint)ceil(5.79 / w_t);
}
// Calculates the beta value of the Kaiser window. Rejection is in dB.
static double CalcKaiserBeta(const double rejection)
{
if(rejection > 50.0)
return 0.1102 * (rejection - 8.7);
if(rejection >= 21.0)
return (0.5842 * pow(rejection - 21.0, 0.4)) +
(0.07886 * (rejection - 21.0));
return 0.0;
}
/* Calculates a point on the Kaiser-windowed sinc filter for the given half-
* width, beta, gain, and cutoff. The point is specified in non-normalized
* samples, from 0 to M, where M = (2 l + 1).
*
* w(k) 2 p f_t sinc(2 f_t x)
*
* x -- centered sample index (i - l)
* k -- normalized and centered window index (x / l)
* w(k) -- window function (Kaiser)
* p -- gain compensation factor when sampling
* f_t -- normalized center frequency (or cutoff; 0.5 is nyquist)
*/
static double SincFilter(const int l, const double b, const double gain, const double cutoff, const int i)
{
return Kaiser(b, (double)(i - l) / l) * 2.0 * gain * cutoff * Sinc(2.0 * cutoff * (i - l));
}
/* This is a polyphase sinc-filtered resampler.
*
* Upsample Downsample
*
* p/q = 3/2 p/q = 3/5
*
* M-+-+-+-> M-+-+-+->
* -------------------+ ---------------------+
* p s * f f f f|f| | p s * f f f f f |
* | 0 * 0 0 0|0|0 | | 0 * 0 0 0 0|0| |
* v 0 * 0 0|0|0 0 | v 0 * 0 0 0|0|0 |
* s * f|f|f f f | s * f f|f|f f |
* 0 * |0|0 0 0 0 | 0 * 0|0|0 0 0 |
* --------+=+--------+ 0 * |0|0 0 0 0 |
* d . d .|d|. d . d ----------+=+--------+
* d . . . .|d|. . . .
* q->
* q-+-+-+->
*
* P_f(i,j) = q i mod p + pj
* P_s(i,j) = floor(q i / p) - j
* d[i=0..N-1] = sum_{j=0}^{floor((M - 1) / p)} {
* { f[P_f(i,j)] s[P_s(i,j)], P_f(i,j) < M
* { 0, P_f(i,j) >= M. }
*/
// Calculate the resampling metrics and build the Kaiser-windowed sinc filter
// that's used to cut frequencies above the destination nyquist.
static void ResamplerSetup(ResamplerT *rs, const uint srcRate, const uint dstRate)
{
double cutoff, width, beta;
uint gcd, l;
int i;
gcd = Gcd(srcRate, dstRate);
rs->mP = dstRate / gcd;
rs->mQ = srcRate / gcd;
/* The cutoff is adjusted by half the transition width, so the transition
* ends before the nyquist (0.5). Both are scaled by the downsampling
* factor.
*/
if(rs->mP > rs->mQ)
{
cutoff = 0.475 / rs->mP;
width = 0.05 / rs->mP;
}
else
{
cutoff = 0.475 / rs->mQ;
width = 0.05 / rs->mQ;
}
// A rejection of -180 dB is used for the stop band. Round up when
// calculating the left offset to avoid increasing the transition width.
l = (CalcKaiserOrder(180.0, width)+1) / 2;
beta = CalcKaiserBeta(180.0);
rs->mM = l*2 + 1;
rs->mL = l;
rs->mF = CreateDoubles(rs->mM);
for(i = 0;i < ((int)rs->mM);i++)
rs->mF[i] = SincFilter((int)l, beta, rs->mP, cutoff, i);
}
// Clean up after the resampler.
static void ResamplerClear(ResamplerT *rs)
{
free(rs->mF);
rs->mF = NULL;
}
// Perform the upsample-filter-downsample resampling operation using a
// polyphase filter implementation.
static void ResamplerRun(ResamplerT *rs, const uint inN, const double *in, const uint outN, double *out)
{
const uint p = rs->mP, q = rs->mQ, m = rs->mM, l = rs->mL;
const double *f = rs->mF;
uint j_f, j_s;
double *work;
uint i;
if(outN == 0)
return;
// Handle in-place operation.
if(in == out)
work = CreateDoubles(outN);
else
work = out;
// Resample the input.
for(i = 0;i < outN;i++)
{
double r = 0.0;
// Input starts at l to compensate for the filter delay. This will
// drop any build-up from the first half of the filter.
j_f = (l + (q * i)) % p;
j_s = (l + (q * i)) / p;
while(j_f < m)
{
// Only take input when 0 <= j_s < inN. This single unsigned
// comparison catches both cases.
if(j_s < inN)
r += f[j_f] * in[j_s];
j_f += p;
j_s--;
}
work[i] = r;
}
// Clean up after in-place operation.
if(work != out)
{
for(i = 0;i < outN;i++)
out[i] = work[i];
free(work);
}
}
/*************************
*** File source input ***
*************************/
// Read a binary value of the specified byte order and byte size from a file,
// storing it as a 32-bit unsigned integer.
static int ReadBin4(FILE *fp, const char *filename, const ByteOrderT order, const uint bytes, uint32 *out)
{
uint8 in[4];
uint32 accum;
uint i;
if(fread(in, 1, bytes, fp) != bytes)
{
fprintf(stderr, "Error: Bad read from file '%s'.\n", filename);
return 0;
}
accum = 0;
switch(order)
{
case BO_LITTLE:
for(i = 0;i < bytes;i++)
accum = (accum<<8) | in[bytes - i - 1];
break;
case BO_BIG:
for(i = 0;i < bytes;i++)
accum = (accum<<8) | in[i];
break;
default:
break;
}
*out = accum;
return 1;
}
// Read a binary value of the specified byte order from a file, storing it as
// a 64-bit unsigned integer.
static int ReadBin8(FILE *fp, const char *filename, const ByteOrderT order, uint64 *out)
{
uint8 in [8];
uint64 accum;
uint i;
if(fread(in, 1, 8, fp) != 8)
{
fprintf(stderr, "Error: Bad read from file '%s'.\n", filename);
return 0;
}
accum = 0ULL;
switch(order)
{
case BO_LITTLE:
for(i = 0;i < 8;i++)
accum = (accum<<8) | in[8 - i - 1];
break;
case BO_BIG:
for(i = 0;i < 8;i++)
accum = (accum<<8) | in[i];
break;
default:
break;
}
*out = accum;
return 1;
}
/* Read a binary value of the specified type, byte order, and byte size from
* a file, converting it to a double. For integer types, the significant
* bits are used to normalize the result. The sign of bits determines
* whether they are padded toward the MSB (negative) or LSB (positive).
* Floating-point types are not normalized.
*/
static int ReadBinAsDouble(FILE *fp, const char *filename, const ByteOrderT order, const ElementTypeT type, const uint bytes, const int bits, double *out)
{
union {
uint32 ui;
int32 i;
float f;
} v4;
union {
uint64 ui;
double f;
} v8;
*out = 0.0;
if(bytes > 4)
{
if(!ReadBin8(fp, filename, order, &v8.ui))
return 0;
if(type == ET_FP)
*out = v8.f;
}
else
{
if(!ReadBin4(fp, filename, order, bytes, &v4.ui))
return 0;
if(type == ET_FP)
*out = v4.f;
else
{
if(bits > 0)
v4.ui >>= (8*bytes) - ((uint)bits);
else
v4.ui &= (0xFFFFFFFF >> (32+bits));
if(v4.ui&(uint)(1<<(abs(bits)-1)))
v4.ui |= (0xFFFFFFFF << abs (bits));
*out = v4.i / (double)(1<<(abs(bits)-1));
}
}
return 1;
}
/* Read an ascii value of the specified type from a file, converting it to a
* double. For integer types, the significant bits are used to normalize the
* result. The sign of the bits should always be positive. This also skips
* up to one separator character before the element itself.
*/
static int ReadAsciiAsDouble(TokenReaderT *tr, const char *filename, const ElementTypeT type, const uint bits, double *out)
{
if(TrIsOperator(tr, ","))
TrReadOperator(tr, ",");
else if(TrIsOperator(tr, ":"))
TrReadOperator(tr, ":");
else if(TrIsOperator(tr, ";"))
TrReadOperator(tr, ";");
else if(TrIsOperator(tr, "|"))
TrReadOperator(tr, "|");
if(type == ET_FP)
{
if(!TrReadFloat(tr, -HUGE_VAL, HUGE_VAL, out))
{
fprintf(stderr, "Error: Bad read from file '%s'.\n", filename);
return 0;
}
}
else
{
int v;
if(!TrReadInt(tr, -(1<<(bits-1)), (1<<(bits-1))-1, &v))
{
fprintf(stderr, "Error: Bad read from file '%s'.\n", filename);
return 0;
}
*out = v / (double)((1<<(bits-1))-1);
}
return 1;
}
// Read the RIFF/RIFX WAVE format chunk from a file, validating it against
// the source parameters and data set metrics.
static int ReadWaveFormat(FILE *fp, const ByteOrderT order, const uint hrirRate, SourceRefT *src)
{
uint32 fourCC, chunkSize;
uint32 format, channels, rate, dummy, block, size, bits;
chunkSize = 0;
do {
if(chunkSize > 0)
fseek (fp, (long) chunkSize, SEEK_CUR);
if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC) ||
!ReadBin4(fp, src->mPath, order, 4, &chunkSize))
return 0;
} while(fourCC != FOURCC_FMT);
if(!ReadBin4(fp, src->mPath, order, 2, &format) ||
!ReadBin4(fp, src->mPath, order, 2, &channels) ||
!ReadBin4(fp, src->mPath, order, 4, &rate) ||
!ReadBin4(fp, src->mPath, order, 4, &dummy) ||
!ReadBin4(fp, src->mPath, order, 2, &block))
return 0;
block /= channels;
if(chunkSize > 14)
{
if(!ReadBin4(fp, src->mPath, order, 2, &size))
return 0;
size /= 8;
if(block > size)
size = block;
}
else
size = block;
if(format == WAVE_FORMAT_EXTENSIBLE)
{
fseek(fp, 2, SEEK_CUR);
if(!ReadBin4(fp, src->mPath, order, 2, &bits))
return 0;
if(bits == 0)
bits = 8 * size;
fseek(fp, 4, SEEK_CUR);
if(!ReadBin4(fp, src->mPath, order, 2, &format))
return 0;
fseek(fp, (long)(chunkSize - 26), SEEK_CUR);
}
else
{
bits = 8 * size;
if(chunkSize > 14)
fseek(fp, (long)(chunkSize - 16), SEEK_CUR);
else
fseek(fp, (long)(chunkSize - 14), SEEK_CUR);
}
if(format != WAVE_FORMAT_PCM && format != WAVE_FORMAT_IEEE_FLOAT)
{
fprintf(stderr, "Error: Unsupported WAVE format in file '%s'.\n", src->mPath);
return 0;
}
if(src->mChannel >= channels)
{
fprintf(stderr, "Error: Missing source channel in WAVE file '%s'.\n", src->mPath);
return 0;
}
if(rate != hrirRate)
{
fprintf(stderr, "Error: Mismatched source sample rate in WAVE file '%s'.\n", src->mPath);
return 0;
}
if(format == WAVE_FORMAT_PCM)
{
if(size < 2 || size > 4)
{
fprintf(stderr, "Error: Unsupported sample size in WAVE file '%s'.\n", src->mPath);
return 0;
}
if(bits < 16 || bits > (8*size))
{
fprintf (stderr, "Error: Bad significant bits in WAVE file '%s'.\n", src->mPath);
return 0;
}
src->mType = ET_INT;
}
else
{
if(size != 4 && size != 8)
{
fprintf(stderr, "Error: Unsupported sample size in WAVE file '%s'.\n", src->mPath);
return 0;
}
src->mType = ET_FP;
}
src->mSize = size;
src->mBits = (int)bits;
src->mSkip = channels;
return 1;
}
// Read a RIFF/RIFX WAVE data chunk, converting all elements to doubles.
static int ReadWaveData(FILE *fp, const SourceRefT *src, const ByteOrderT order, const uint n, double *hrir)
{
int pre, post, skip;
uint i;
pre = (int)(src->mSize * src->mChannel);
post = (int)(src->mSize * (src->mSkip - src->mChannel - 1));
skip = 0;
for(i = 0;i < n;i++)
{
skip += pre;
if(skip > 0)
fseek(fp, skip, SEEK_CUR);
if(!ReadBinAsDouble(fp, src->mPath, order, src->mType, src->mSize, src->mBits, &hrir[i]))
return 0;
skip = post;
}
if(skip > 0)
fseek(fp, skip, SEEK_CUR);
return 1;
}
// Read the RIFF/RIFX WAVE list or data chunk, converting all elements to
// doubles.
static int ReadWaveList(FILE *fp, const SourceRefT *src, const ByteOrderT order, const uint n, double *hrir)
{
uint32 fourCC, chunkSize, listSize, count;
uint block, skip, offset, i;
double lastSample;
for(;;)
{
if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC) ||
!ReadBin4(fp, src->mPath, order, 4, &chunkSize))
return 0;
if(fourCC == FOURCC_DATA)
{
block = src->mSize * src->mSkip;
count = chunkSize / block;
if(count < (src->mOffset + n))
{
fprintf(stderr, "Error: Bad read from file '%s'.\n", src->mPath);
return 0;
}
fseek(fp, (long)(src->mOffset * block), SEEK_CUR);
if(!ReadWaveData(fp, src, order, n, &hrir[0]))
return 0;
return 1;
}
else if(fourCC == FOURCC_LIST)
{
if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC))
return 0;
chunkSize -= 4;
if(fourCC == FOURCC_WAVL)
break;
}
if(chunkSize > 0)
fseek(fp, (long)chunkSize, SEEK_CUR);
}
listSize = chunkSize;
block = src->mSize * src->mSkip;
skip = src->mOffset;
offset = 0;
lastSample = 0.0;
while(offset < n && listSize > 8)
{
if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC) ||
!ReadBin4(fp, src->mPath, order, 4, &chunkSize))
return 0;
listSize -= 8 + chunkSize;
if(fourCC == FOURCC_DATA)
{
count = chunkSize / block;
if(count > skip)
{
fseek(fp, (long)(skip * block), SEEK_CUR);
chunkSize -= skip * block;
count -= skip;
skip = 0;
if(count > (n - offset))
count = n - offset;
if(!ReadWaveData(fp, src, order, count, &hrir[offset]))
return 0;
chunkSize -= count * block;
offset += count;
lastSample = hrir [offset - 1];
}
else
{
skip -= count;
count = 0;
}
}
else if(fourCC == FOURCC_SLNT)
{
if(!ReadBin4(fp, src->mPath, order, 4, &count))
return 0;
chunkSize -= 4;
if(count > skip)
{
count -= skip;
skip = 0;
if(count > (n - offset))
count = n - offset;
for(i = 0; i < count; i ++)
hrir[offset + i] = lastSample;
offset += count;
}
else
{
skip -= count;
count = 0;
}
}
if(chunkSize > 0)
fseek(fp, (long)chunkSize, SEEK_CUR);
}
if(offset < n)
{
fprintf(stderr, "Error: Bad read from file '%s'.\n", src->mPath);
return 0;
}
return 1;
}
// Load a source HRIR from a RIFF/RIFX WAVE file.
static int LoadWaveSource(FILE *fp, SourceRefT *src, const uint hrirRate, const uint n, double *hrir)
{
uint32 fourCC, dummy;
ByteOrderT order;
if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC) ||
!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &dummy))
return 0;
if(fourCC == FOURCC_RIFF)
order = BO_LITTLE;
else if(fourCC == FOURCC_RIFX)
order = BO_BIG;
else
{
fprintf(stderr, "Error: No RIFF/RIFX chunk in file '%s'.\n", src->mPath);
return 0;
}
if(!ReadBin4(fp, src->mPath, BO_LITTLE, 4, &fourCC))
return 0;
if(fourCC != FOURCC_WAVE)
{
fprintf(stderr, "Error: Not a RIFF/RIFX WAVE file '%s'.\n", src->mPath);
return 0;
}
if(!ReadWaveFormat(fp, order, hrirRate, src))
return 0;
if(!ReadWaveList(fp, src, order, n, hrir))
return 0;
return 1;
}
// Load a source HRIR from a binary file.
static int LoadBinarySource(FILE *fp, const SourceRefT *src, const ByteOrderT order, const uint n, double *hrir)
{
uint i;
fseek(fp, (long)src->mOffset, SEEK_SET);
for(i = 0;i < n;i++)
{
if(!ReadBinAsDouble(fp, src->mPath, order, src->mType, src->mSize, src->mBits, &hrir[i]))
return 0;
if(src->mSkip > 0)
fseek(fp, (long)src->mSkip, SEEK_CUR);
}
return 1;
}
// Load a source HRIR from an ASCII text file containing a list of elements
// separated by whitespace or common list operators (',', ';', ':', '|').
static int LoadAsciiSource(FILE *fp, const SourceRefT *src, const uint n, double *hrir)
{
TokenReaderT tr;
uint i, j;
double dummy;
TrSetup(fp, NULL, &tr);
for(i = 0;i < src->mOffset;i++)
{
if(!ReadAsciiAsDouble(&tr, src->mPath, src->mType, (uint)src->mBits, &dummy))
return 0;
}
for(i = 0;i < n;i++)
{
if(!ReadAsciiAsDouble(&tr, src->mPath, src->mType, (uint)src->mBits, &hrir[i]))
return 0;
for(j = 0;j < src->mSkip;j++)
{
if(!ReadAsciiAsDouble(&tr, src->mPath, src->mType, (uint)src->mBits, &dummy))
return 0;
}
}
return 1;
}
// Load a source HRIR from a supported file type.
static int LoadSource(SourceRefT *src, const uint hrirRate, const uint n, double *hrir)
{
int result;
FILE *fp;
if(src->mFormat == SF_ASCII)
fp = fopen(src->mPath, "r");
else
fp = fopen(src->mPath, "rb");
if(fp == NULL)
{
fprintf(stderr, "Error: Could not open source file '%s'.\n", src->mPath);
return 0;
}
if(src->mFormat == SF_WAVE)
result = LoadWaveSource(fp, src, hrirRate, n, hrir);
else if(src->mFormat == SF_BIN_LE)
result = LoadBinarySource(fp, src, BO_LITTLE, n, hrir);
else if(src->mFormat == SF_BIN_BE)
result = LoadBinarySource(fp, src, BO_BIG, n, hrir);
else
result = LoadAsciiSource(fp, src, n, hrir);
fclose(fp);
return result;
}
/***************************
*** File storage output ***
***************************/
// Write an ASCII string to a file.
static int WriteAscii(const char *out, FILE *fp, const char *filename)
{
size_t len;
len = strlen(out);
if(fwrite(out, 1, len, fp) != len)
{
fclose(fp);
fprintf(stderr, "Error: Bad write to file '%s'.\n", filename);
return 0;
}
return 1;
}
// Write a binary value of the given byte order and byte size to a file,
// loading it from a 32-bit unsigned integer.
static int WriteBin4(const ByteOrderT order, const uint bytes, const uint32 in, FILE *fp, const char *filename)
{
uint8 out[4];
uint i;
switch(order)
{
case BO_LITTLE:
for(i = 0;i < bytes;i++)
out[i] = (in>>(i*8)) & 0x000000FF;
break;
case BO_BIG:
for(i = 0;i < bytes;i++)
out[bytes - i - 1] = (in>>(i*8)) & 0x000000FF;
break;
default:
break;
}
if(fwrite(out, 1, bytes, fp) != bytes)
{
fprintf(stderr, "Error: Bad write to file '%s'.\n", filename);
return 0;
}
return 1;
}
// Store the OpenAL Soft HRTF data set.
static int StoreMhr(const HrirDataT *hData, const char *filename)
{
uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
uint n = hData->mIrPoints;
FILE *fp;
uint fi, ei, ai, i;
uint dither_seed = 22222;
if((fp=fopen(filename, "wb")) == NULL)
{
fprintf(stderr, "Error: Could not open MHR file '%s'.\n", filename);
return 0;
}
if(!WriteAscii(MHR_FORMAT, fp, filename))
return 0;
if(!WriteBin4(BO_LITTLE, 4, (uint32)hData->mIrRate, fp, filename))
return 0;
if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mSampleType, fp, filename))
return 0;
if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mChannelType, fp, filename))
return 0;
if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mIrPoints, fp, filename))
return 0;
if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mFdCount, fp, filename))
return 0;
for(fi = 0;fi < hData->mFdCount;fi++)
{
if(!WriteBin4(BO_LITTLE, 2, (uint32)(1000.0 * hData->mFds[fi].mDistance), fp, filename))
return 0;
if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mFds[fi].mEvCount, fp, filename))
return 0;
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
if(!WriteBin4(BO_LITTLE, 1, (uint32)hData->mFds[fi].mEvs[ei].mAzCount, fp, filename))
return 0;
}
}
for(fi = 0;fi < hData->mFdCount;fi++)
{
const double scale = (hData->mSampleType == ST_S16) ? 32767.0 :
((hData->mSampleType == ST_S24) ? 8388607.0 : 0.0);
const int bps = (hData->mSampleType == ST_S16) ? 2 :
((hData->mSampleType == ST_S24) ? 3 : 0);
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
double out[2 * MAX_TRUNCSIZE];
TpdfDither(out, azd->mIrs[0], scale, n, channels, &dither_seed);
if(hData->mChannelType == CT_STEREO)
TpdfDither(out+1, azd->mIrs[1], scale, n, channels, &dither_seed);
for(i = 0;i < (channels * n);i++)
{
int v = (int)Clamp(out[i], -scale-1.0, scale);
if(!WriteBin4(BO_LITTLE, bps, (uint32)v, fp, filename))
return 0;
}
}
}
}
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
int v = (int)fmin(round(hData->mIrRate * azd->mDelays[0]), MAX_HRTD);
if(!WriteBin4(BO_LITTLE, 1, (uint32)v, fp, filename))
return 0;
if(hData->mChannelType == CT_STEREO)
{
v = (int)fmin(round(hData->mIrRate * azd->mDelays[1]), MAX_HRTD);
if(!WriteBin4(BO_LITTLE, 1, (uint32)v, fp, filename))
return 0;
}
}
}
}
fclose(fp);
return 1;
}
/***********************
*** HRTF processing ***
***********************/
// Calculate the onset time of an HRIR and average it with any existing
// timing for its field, elevation, azimuth, and ear.
static double AverageHrirOnset(const uint rate, const uint n, const double *hrir, const double f, const double onset)
{
double mag = 0.0;
uint i;
for(i = 0;i < n;i++)
mag = fmax(fabs(hrir[i]), mag);
mag *= 0.15;
for(i = 0;i < n;i++)
{
if(fabs(hrir[i]) >= mag)
break;
}
return Lerp(onset, (double)i / rate, f);
}
// Calculate the magnitude response of an HRIR and average it with any
// existing responses for its field, elevation, azimuth, and ear.
static void AverageHrirMagnitude(const uint points, const uint n, const double *hrir, const double f, double *mag)
{
uint m = 1 + (n / 2), i;
Complex *h = CreateComplexes(n);
double *r = CreateDoubles(n);
for(i = 0;i < points;i++)
h[i] = MakeComplex(hrir[i], 0.0);
for(;i < n;i++)
h[i] = MakeComplex(0.0, 0.0);
FftForward(n, h, h);
MagnitudeResponse(n, h, r);
for(i = 0;i < m;i++)
mag[i] = Lerp(mag[i], r[i], f);
free(r);
free(h);
}
/* Calculate the contribution of each HRIR to the diffuse-field average based
* on the area of its surface patch. All patches are centered at the HRIR
* coordinates on the unit sphere and are measured by solid angle.
*/
static void CalculateDfWeights(const HrirDataT *hData, double *weights)
{
double sum, evs, ev, upperEv, lowerEv, solidAngle;
uint fi, ei;
sum = 0.0;
for(fi = 0;fi < hData->mFdCount;fi++)
{
evs = M_PI / 2.0 / (hData->mFds[fi].mEvCount - 1);
for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvCount;ei++)
{
// For each elevation, calculate the upper and lower limits of
// the patch band.
ev = hData->mFds[fi].mEvs[ei].mElevation;
lowerEv = fmax(-M_PI / 2.0, ev - evs);
upperEv = fmin(M_PI / 2.0, ev + evs);
// Calculate the area of the patch band.
solidAngle = 2.0 * M_PI * (sin(upperEv) - sin(lowerEv));
// Each weight is the area of one patch.
weights[(fi * MAX_EV_COUNT) + ei] = solidAngle / hData->mFds[fi].mEvs[ei].mAzCount;
// Sum the total surface area covered by the HRIRs of all fields.
sum += solidAngle;
}
}
/* TODO: It may be interesting to experiment with how a volume-based
weighting performs compared to the existing distance-indepenent
surface patches.
*/
for(fi = 0;fi < hData->mFdCount;fi++)
{
// Normalize the weights given the total surface coverage for all
// fields.
for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvCount;ei++)
weights[(fi * MAX_EV_COUNT) + ei] /= sum;
}
}
/* Calculate the diffuse-field average from the given magnitude responses of
* the HRIR set. Weighting can be applied to compensate for the varying
* surface area covered by each HRIR. The final average can then be limited
* by the specified magnitude range (in positive dB; 0.0 to skip).
*/
static void CalculateDiffuseFieldAverage(const HrirDataT *hData, const uint channels, const uint m, const int weighted, const double limit, double *dfa)
{
double *weights = CreateDoubles(hData->mFdCount * MAX_EV_COUNT);
uint count, ti, fi, ei, i, ai;
if(weighted)
{
// Use coverage weighting to calculate the average.
CalculateDfWeights(hData, weights);
}
else
{
double weight;
// If coverage weighting is not used, the weights still need to be
// averaged by the number of existing HRIRs.
count = hData->mIrCount;
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvStart;ei++)
count -= hData->mFds[fi].mEvs[ei].mAzCount;
}
weight = 1.0 / count;
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvCount;ei++)
weights[(fi * MAX_EV_COUNT) + ei] = weight;
}
}
for(ti = 0;ti < channels;ti++)
{
for(i = 0;i < m;i++)
dfa[(ti * m) + i] = 0.0;
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
// Get the weight for this HRIR's contribution.
double weight = weights[(fi * MAX_EV_COUNT) + ei];
// Add this HRIR's weighted power average to the total.
for(i = 0;i < m;i++)
dfa[(ti * m) + i] += weight * azd->mIrs[ti][i] * azd->mIrs[ti][i];
}
}
}
// Finish the average calculation and keep it from being too small.
for(i = 0;i < m;i++)
dfa[(ti * m) + i] = fmax(sqrt(dfa[(ti * m) + i]), EPSILON);
// Apply a limit to the magnitude range of the diffuse-field average
// if desired.
if(limit > 0.0)
LimitMagnitudeResponse(hData->mFftSize, m, limit, &dfa[ti * m], &dfa[ti * m]);
}
free(weights);
}
// Perform diffuse-field equalization on the magnitude responses of the HRIR
// set using the given average response.
static void DiffuseFieldEqualize(const uint channels, const uint m, const double *dfa, const HrirDataT *hData)
{
uint ti, fi, ei, ai, i;
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
for(i = 0;i < m;i++)
azd->mIrs[ti][i] /= dfa[(ti * m) + i];
}
}
}
}
}
// Perform minimum-phase reconstruction using the magnitude responses of the
// HRIR set.
static void ReconstructHrirs(const HrirDataT *hData)
{
uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
uint n = hData->mFftSize;
uint ti, fi, ei, ai, i;
Complex *h = CreateComplexes(n);
uint total, count, pcdone, lastpc;
total = hData->mIrCount;
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvStart;ei++)
total -= hData->mFds[fi].mEvs[ei].mAzCount;
}
total *= channels;
count = pcdone = lastpc = 0;
printf("%3d%% done.", pcdone);
fflush(stdout);
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
MinimumPhase(n, azd->mIrs[ti], h);
FftInverse(n, h, h);
for(i = 0;i < hData->mIrPoints;i++)
azd->mIrs[ti][i] = h[i].Real;
pcdone = ++count * 100 / total;
if(pcdone != lastpc)
{
lastpc = pcdone;
printf("\r%3d%% done.", pcdone);
fflush(stdout);
}
}
}
}
}
printf("\n");
free(h);
}
// Resamples the HRIRs for use at the given sampling rate.
static void ResampleHrirs(const uint rate, HrirDataT *hData)
{
uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
uint n = hData->mIrPoints;
uint ti, fi, ei, ai;
ResamplerT rs;
ResamplerSetup(&rs, hData->mIrRate, rate);
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = hData->mFds[fi].mEvStart;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
ResamplerRun(&rs, n, azd->mIrs[ti], n, azd->mIrs[ti]);
}
}
}
}
hData->mIrRate = rate;
ResamplerClear(&rs);
}
/* Given field and elevation indices and an azimuth, calculate the indices of
* the two HRIRs that bound the coordinate along with a factor for
* calculating the continuous HRIR using interpolation.
*/
static void CalcAzIndices(const HrirDataT *hData, const uint fi, const uint ei, const double az, uint *a0, uint *a1, double *af)
{
double f = (2.0*M_PI + az) * hData->mFds[fi].mEvs[ei].mAzCount / (2.0*M_PI);
uint i = (uint)f % hData->mFds[fi].mEvs[ei].mAzCount;
f -= floor(f);
*a0 = i;
*a1 = (i + 1) % hData->mFds[fi].mEvs[ei].mAzCount;
*af = f;
}
// Synthesize any missing onset timings at the bottom elevations of each
// field. This just blends between slightly exaggerated known onsets (not
// an accurate model).
static void SynthesizeOnsets(HrirDataT *hData)
{
uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
uint ti, fi, oi, ai, ei, a0, a1;
double t, of, af;
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
if(hData->mFds[fi].mEvStart <= 0)
continue;
oi = hData->mFds[fi].mEvStart;
t = 0.0;
for(ai = 0;ai < hData->mFds[fi].mEvs[oi].mAzCount;ai++)
t += hData->mFds[fi].mEvs[oi].mAzs[ai].mDelays[ti];
hData->mFds[fi].mEvs[0].mAzs[0].mDelays[ti] = 1.32e-4 + (t / hData->mFds[fi].mEvs[oi].mAzCount);
for(ei = 1;ei < hData->mFds[fi].mEvStart;ei++)
{
of = (double)ei / hData->mFds[fi].mEvStart;
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
CalcAzIndices(hData, fi, oi, hData->mFds[fi].mEvs[ei].mAzs[ai].mAzimuth, &a0, &a1, &af);
hData->mFds[fi].mEvs[ei].mAzs[ai].mDelays[ti] = Lerp(hData->mFds[fi].mEvs[0].mAzs[0].mDelays[ti], Lerp(hData->mFds[fi].mEvs[oi].mAzs[a0].mDelays[ti], hData->mFds[fi].mEvs[oi].mAzs[a1].mDelays[ti], af), of);
}
}
}
}
}
/* Attempt to synthesize any missing HRIRs at the bottom elevations of each
* field. Right now this just blends the lowest elevation HRIRs together and
* applies some attenuation and high frequency damping. It is a simple, if
* inaccurate model.
*/
static void SynthesizeHrirs(HrirDataT *hData)
{
uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
uint n = hData->mIrPoints;
uint ti, fi, oi, ai, ei, i;
double lp[4], s0, s1;
double of, b;
uint a0, a1;
double af;
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
if(hData->mFds[fi].mEvStart <= 0)
continue;
oi = hData->mFds[fi].mEvStart;
for(i = 0;i < n;i++)
hData->mFds[fi].mEvs[0].mAzs[0].mIrs[ti][i] = 0.0;
for(ai = 0;ai < hData->mFds[fi].mEvs[oi].mAzCount;ai++)
{
for(i = 0;i < n;i++)
hData->mFds[fi].mEvs[0].mAzs[0].mIrs[ti][i] += hData->mFds[fi].mEvs[oi].mAzs[ai].mIrs[ti][i] / hData->mFds[fi].mEvs[oi].mAzCount;
}
for(ei = 1;ei < hData->mFds[fi].mEvStart;ei++)
{
of = (double)ei / hData->mFds[fi].mEvStart;
b = (1.0 - of) * (3.5e-6 * hData->mIrRate);
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
CalcAzIndices(hData, fi, oi, hData->mFds[fi].mEvs[ei].mAzs[ai].mAzimuth, &a0, &a1, &af);
lp[0] = 0.0;
lp[1] = 0.0;
lp[2] = 0.0;
lp[3] = 0.0;
for(i = 0;i < n;i++)
{
s0 = hData->mFds[fi].mEvs[0].mAzs[0].mIrs[ti][i];
s1 = Lerp(hData->mFds[fi].mEvs[oi].mAzs[a0].mIrs[ti][i], hData->mFds[fi].mEvs[oi].mAzs[a1].mIrs[ti][i], af);
s0 = Lerp(s0, s1, of);
lp[0] = Lerp(s0, lp[0], b);
lp[1] = Lerp(lp[0], lp[1], b);
lp[2] = Lerp(lp[1], lp[2], b);
lp[3] = Lerp(lp[2], lp[3], b);
hData->mFds[fi].mEvs[ei].mAzs[ai].mIrs[ti][i] = lp[3];
}
}
}
b = 3.5e-6 * hData->mIrRate;
lp[0] = 0.0;
lp[1] = 0.0;
lp[2] = 0.0;
lp[3] = 0.0;
for(i = 0;i < n;i++)
{
s0 = hData->mFds[fi].mEvs[0].mAzs[0].mIrs[ti][i];
lp[0] = Lerp(s0, lp[0], b);
lp[1] = Lerp(lp[0], lp[1], b);
lp[2] = Lerp(lp[1], lp[2], b);
lp[3] = Lerp(lp[2], lp[3], b);
hData->mFds[fi].mEvs[0].mAzs[0].mIrs[ti][i] = lp[3];
}
hData->mFds[fi].mEvStart = 0;
}
}
}
// The following routines assume a full set of HRIRs for all elevations.
// Normalize the HRIR set and slightly attenuate the result.
static void NormalizeHrirs(const HrirDataT *hData)
{
uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
uint n = hData->mIrPoints;
uint ti, fi, ei, ai, i;
double maxLevel = 0.0;
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
for(i = 0;i < n;i++)
maxLevel = fmax(fabs(azd->mIrs[ti][i]), maxLevel);
}
}
}
}
maxLevel = 1.01 * maxLevel;
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
for(i = 0;i < n;i++)
azd->mIrs[ti][i] /= maxLevel;
}
}
}
}
}
// Calculate the left-ear time delay using a spherical head model.
static double CalcLTD(const double ev, const double az, const double rad, const double dist)
{
double azp, dlp, l, al;
azp = asin(cos(ev) * sin(az));
dlp = sqrt((dist*dist) + (rad*rad) + (2.0*dist*rad*sin(azp)));
l = sqrt((dist*dist) - (rad*rad));
al = (0.5 * M_PI) + azp;
if(dlp > l)
dlp = l + (rad * (al - acos(rad / dist)));
return dlp / 343.3;
}
// Calculate the effective head-related time delays for each minimum-phase
// HRIR.
static void CalculateHrtds(const HeadModelT model, const double radius, HrirDataT *hData)
{
uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
double minHrtd = INFINITY, maxHrtd = -INFINITY;
uint ti, fi, ei, ai;
double t;
if(model == HM_DATASET)
{
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
t = azd->mDelays[ti] * radius / hData->mRadius;
azd->mDelays[ti] = t;
maxHrtd = fmax(t, maxHrtd);
minHrtd = fmin(t, minHrtd);
}
}
}
}
}
else
{
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
HrirEvT *evd = &hData->mFds[fi].mEvs[ei];
for(ai = 0;ai < evd->mAzCount;ai++)
{
HrirAzT *azd = &evd->mAzs[ai];
t = CalcLTD(evd->mElevation, azd->mAzimuth, radius, hData->mFds[fi].mDistance);
azd->mDelays[ti] = t;
maxHrtd = fmax(t, maxHrtd);
minHrtd = fmin(t, minHrtd);
}
}
}
}
}
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
hData->mFds[fi].mEvs[ei].mAzs[ai].mDelays[ti] -= minHrtd;
}
}
}
}
// Clear the initial HRIR data state.
static void ResetHrirData(HrirDataT *hData)
{
hData->mIrRate = 0;
hData->mSampleType = ST_S24;
hData->mChannelType = CT_NONE;
hData->mIrPoints = 0;
hData->mFftSize = 0;
hData->mIrSize = 0;
hData->mRadius = 0.0;
hData->mIrCount = 0;
hData->mFdCount = 0;
hData->mFds = NULL;
}
// Allocate and configure dynamic HRIR structures.
static int PrepareHrirData(const uint fdCount, const double distances[MAX_FD_COUNT], const uint evCounts[MAX_FD_COUNT], const uint azCounts[MAX_FD_COUNT * MAX_EV_COUNT], HrirDataT *hData)
{
uint evTotal = 0, azTotal = 0, fi, ei, ai;
for(fi = 0;fi < fdCount;fi++)
{
evTotal += evCounts[fi];
for(ei = 0;ei < evCounts[fi];ei++)
azTotal += azCounts[(fi * MAX_EV_COUNT) + ei];
}
if(!fdCount || !evTotal || !azTotal)
return 0;
hData->mFds = calloc(fdCount, sizeof(*hData->mFds));
if(hData->mFds == NULL)
return 0;
hData->mFds[0].mEvs = calloc(evTotal, sizeof(*hData->mFds[0].mEvs));
if(hData->mFds[0].mEvs == NULL)
return 0;
hData->mFds[0].mEvs[0].mAzs = calloc(azTotal, sizeof(*hData->mFds[0].mEvs[0].mAzs));
if(hData->mFds[0].mEvs[0].mAzs == NULL)
return 0;
hData->mIrCount = azTotal;
hData->mFdCount = fdCount;
evTotal = 0;
azTotal = 0;
for(fi = 0;fi < fdCount;fi++)
{
hData->mFds[fi].mDistance = distances[fi];
hData->mFds[fi].mEvCount = evCounts[fi];
hData->mFds[fi].mEvStart = 0;
hData->mFds[fi].mEvs = &hData->mFds[0].mEvs[evTotal];
evTotal += evCounts[fi];
for(ei = 0;ei < evCounts[fi];ei++)
{
uint azCount = azCounts[(fi * MAX_EV_COUNT) + ei];
hData->mFds[fi].mIrCount += azCount;
hData->mFds[fi].mEvs[ei].mElevation = -M_PI / 2.0 + M_PI * ei / (evCounts[fi] - 1);
hData->mFds[fi].mEvs[ei].mIrCount += azCount;
hData->mFds[fi].mEvs[ei].mAzCount = azCount;
hData->mFds[fi].mEvs[ei].mAzs = &hData->mFds[0].mEvs[0].mAzs[azTotal];
for(ai = 0;ai < azCount;ai++)
{
hData->mFds[fi].mEvs[ei].mAzs[ai].mAzimuth = 2.0 * M_PI * ai / azCount;
hData->mFds[fi].mEvs[ei].mAzs[ai].mIndex = azTotal + ai;
hData->mFds[fi].mEvs[ei].mAzs[ai].mDelays[0] = 0.0;
hData->mFds[fi].mEvs[ei].mAzs[ai].mDelays[1] = 0.0;
hData->mFds[fi].mEvs[ei].mAzs[ai].mIrs[0] = NULL;
hData->mFds[fi].mEvs[ei].mAzs[ai].mIrs[1] = NULL;
}
azTotal += azCount;
}
}
return 1;
}
// Clean up HRIR data.
static void FreeHrirData(HrirDataT *hData)
{
if(hData->mFds != NULL)
{
if(hData->mFds[0].mEvs != NULL)
{
if(hData->mFds[0].mEvs[0].mAzs)
{
if(hData->mFds[0].mEvs[0].mAzs[0].mIrs[0] != NULL)
free(hData->mFds[0].mEvs[0].mAzs[0].mIrs[0]);
free(hData->mFds[0].mEvs[0].mAzs);
}
free(hData->mFds[0].mEvs);
}
free(hData->mFds);
hData->mFds = NULL;
}
}
// Match the channel type from a given identifier.
static ChannelTypeT MatchChannelType(const char *ident)
{
if(strcasecmp(ident, "mono") == 0)
return CT_MONO;
if(strcasecmp(ident, "stereo") == 0)
return CT_STEREO;
return CT_NONE;
}
// Process the data set definition to read and validate the data set metrics.
static int ProcessMetrics(TokenReaderT *tr, const uint fftSize, const uint truncSize, HrirDataT *hData)
{
int hasRate = 0, hasType = 0, hasPoints = 0, hasRadius = 0;
int hasDistance = 0, hasAzimuths = 0;
char ident[MAX_IDENT_LEN+1];
uint line, col;
double fpVal;
uint points;
int intVal;
double distances[MAX_FD_COUNT];
uint fdCount = 0;
uint evCounts[MAX_FD_COUNT];
uint *azCounts = calloc(MAX_FD_COUNT * MAX_EV_COUNT, sizeof(*azCounts));
if(azCounts == NULL)
{
fprintf(stderr, "Error: Out of memory.\n");
exit(-1);
}
TrIndication(tr, &line, &col);
while(TrIsIdent(tr))
{
TrIndication(tr, &line, &col);
if(!TrReadIdent(tr, MAX_IDENT_LEN, ident))
goto error;
if(strcasecmp(ident, "rate") == 0)
{
if(hasRate)
{
TrErrorAt(tr, line, col, "Redefinition of 'rate'.\n");
goto error;
}
if(!TrReadOperator(tr, "="))
goto error;
if(!TrReadInt(tr, MIN_RATE, MAX_RATE, &intVal))
goto error;
hData->mIrRate = (uint)intVal;
hasRate = 1;
}
else if(strcasecmp(ident, "type") == 0)
{
char type[MAX_IDENT_LEN+1];
if(hasType)
{
TrErrorAt(tr, line, col, "Redefinition of 'type'.\n");
goto error;
}
if(!TrReadOperator(tr, "="))
goto error;
if(!TrReadIdent(tr, MAX_IDENT_LEN, type))
goto error;
hData->mChannelType = MatchChannelType(type);
if(hData->mChannelType == CT_NONE)
{
TrErrorAt(tr, line, col, "Expected a channel type.\n");
goto error;
}
hasType = 1;
}
else if(strcasecmp(ident, "points") == 0)
{
if(hasPoints)
{
TrErrorAt(tr, line, col, "Redefinition of 'points'.\n");
goto error;
}
if(!TrReadOperator(tr, "="))
goto error;
TrIndication(tr, &line, &col);
if(!TrReadInt(tr, MIN_POINTS, MAX_POINTS, &intVal))
goto error;
points = (uint)intVal;
if(fftSize > 0 && points > fftSize)
{
TrErrorAt(tr, line, col, "Value exceeds the overridden FFT size.\n");
goto error;
}
if(points < truncSize)
{
TrErrorAt(tr, line, col, "Value is below the truncation size.\n");
goto error;
}
hData->mIrPoints = points;
if(fftSize <= 0)
{
hData->mFftSize = DEFAULT_FFTSIZE;
hData->mIrSize = 1 + (DEFAULT_FFTSIZE / 2);
}
else
{
hData->mFftSize = fftSize;
hData->mIrSize = 1 + (fftSize / 2);
if(points > hData->mIrSize)
hData->mIrSize = points;
}
hasPoints = 1;
}
else if(strcasecmp(ident, "radius") == 0)
{
if(hasRadius)
{
TrErrorAt(tr, line, col, "Redefinition of 'radius'.\n");
goto error;
}
if(!TrReadOperator(tr, "="))
goto error;
if(!TrReadFloat(tr, MIN_RADIUS, MAX_RADIUS, &fpVal))
goto error;
hData->mRadius = fpVal;
hasRadius = 1;
}
else if(strcasecmp(ident, "distance") == 0)
{
uint count = 0;
if(hasDistance)
{
TrErrorAt(tr, line, col, "Redefinition of 'distance'.\n");
goto error;
}
if(!TrReadOperator(tr, "="))
goto error;
for(;;)
{
if(!TrReadFloat(tr, MIN_DISTANCE, MAX_DISTANCE, &fpVal))
goto error;
if(count > 0 && fpVal <= distances[count - 1])
{
TrError(tr, "Distances are not ascending.\n");
goto error;
}
distances[count++] = fpVal;
if(!TrIsOperator(tr, ","))
break;
if(count >= MAX_FD_COUNT)
{
TrError(tr, "Exceeded the maximum of %d fields.\n", MAX_FD_COUNT);
goto error;
}
TrReadOperator(tr, ",");
}
if(fdCount != 0 && count != fdCount)
{
TrError(tr, "Did not match the specified number of %d fields.\n", fdCount);
goto error;
}
fdCount = count;
hasDistance = 1;
}
else if(strcasecmp(ident, "azimuths") == 0)
{
uint count = 0;
if(hasAzimuths)
{
TrErrorAt(tr, line, col, "Redefinition of 'azimuths'.\n");
goto error;
}
if(!TrReadOperator(tr, "="))
goto error;
evCounts[0] = 0;
for(;;)
{
if(!TrReadInt(tr, MIN_AZ_COUNT, MAX_AZ_COUNT, &intVal))
goto error;
azCounts[(count * MAX_EV_COUNT) + evCounts[count]++] = (uint)intVal;
if(TrIsOperator(tr, ","))
{
if(evCounts[count] >= MAX_EV_COUNT)
{
TrError(tr, "Exceeded the maximum of %d elevations.\n", MAX_EV_COUNT);
goto error;
}
TrReadOperator(tr, ",");
}
else
{
if(evCounts[count] < MIN_EV_COUNT)
{
TrErrorAt(tr, line, col, "Did not reach the minimum of %d azimuth counts.\n", MIN_EV_COUNT);
goto error;
}
if(azCounts[count * MAX_EV_COUNT] != 1 || azCounts[(count * MAX_EV_COUNT) + evCounts[count] - 1] != 1)
{
TrError(tr, "Poles are not singular for field %d.\n", count - 1);
goto error;
}
count++;
if(TrIsOperator(tr, ";"))
{
if(count >= MAX_FD_COUNT)
{
TrError(tr, "Exceeded the maximum number of %d fields.\n", MAX_FD_COUNT);
goto error;
}
evCounts[count] = 0;
TrReadOperator(tr, ";");
}
else
{
break;
}
}
}
if(fdCount != 0 && count != fdCount)
{
TrError(tr, "Did not match the specified number of %d fields.\n", fdCount);
goto error;
}
fdCount = count;
hasAzimuths = 1;
}
else
{
TrErrorAt(tr, line, col, "Expected a metric name.\n");
goto error;
}
TrSkipWhitespace(tr);
}
if(!(hasRate && hasPoints && hasRadius && hasDistance && hasAzimuths))
{
TrErrorAt(tr, line, col, "Expected a metric name.\n");
goto error;
}
if(distances[0] < hData->mRadius)
{
TrError(tr, "Distance cannot start below head radius.\n");
goto error;
}
if(hData->mChannelType == CT_NONE)
hData->mChannelType = CT_MONO;
if(!PrepareHrirData(fdCount, distances, evCounts, azCounts, hData))
{
fprintf(stderr, "Error: Out of memory.\n");
exit(-1);
}
free(azCounts);
return 1;
error:
free(azCounts);
return 0;
}
// Parse an index triplet from the data set definition.
static int ReadIndexTriplet(TokenReaderT *tr, const HrirDataT *hData, uint *fi, uint *ei, uint *ai)
{
int intVal;
if(hData->mFdCount > 1)
{
if(!TrReadInt(tr, 0, (int)hData->mFdCount - 1, &intVal))
return 0;
*fi = (uint)intVal;
if(!TrReadOperator(tr, ","))
return 0;
}
else
{
*fi = 0;
}
if(!TrReadInt(tr, 0, (int)hData->mFds[*fi].mEvCount - 1, &intVal))
return 0;
*ei = (uint)intVal;
if(!TrReadOperator(tr, ","))
return 0;
if(!TrReadInt(tr, 0, (int)hData->mFds[*fi].mEvs[*ei].mAzCount - 1, &intVal))
return 0;
*ai = (uint)intVal;
return 1;
}
// Match the source format from a given identifier.
static SourceFormatT MatchSourceFormat(const char *ident)
{
if(strcasecmp(ident, "wave") == 0)
return SF_WAVE;
if(strcasecmp(ident, "bin_le") == 0)
return SF_BIN_LE;
if(strcasecmp(ident, "bin_be") == 0)
return SF_BIN_BE;
if(strcasecmp(ident, "ascii") == 0)
return SF_ASCII;
return SF_NONE;
}
// Match the source element type from a given identifier.
static ElementTypeT MatchElementType(const char *ident)
{
if(strcasecmp(ident, "int") == 0)
return ET_INT;
if(strcasecmp(ident, "fp") == 0)
return ET_FP;
return ET_NONE;
}
// Parse and validate a source reference from the data set definition.
static int ReadSourceRef(TokenReaderT *tr, SourceRefT *src)
{
char ident[MAX_IDENT_LEN+1];
uint line, col;
int intVal;
TrIndication(tr, &line, &col);
if(!TrReadIdent(tr, MAX_IDENT_LEN, ident))
return 0;
src->mFormat = MatchSourceFormat(ident);
if(src->mFormat == SF_NONE)
{
TrErrorAt(tr, line, col, "Expected a source format.\n");
return 0;
}
if(!TrReadOperator(tr, "("))
return 0;
if(src->mFormat == SF_WAVE)
{
if(!TrReadInt(tr, 0, MAX_WAVE_CHANNELS, &intVal))
return 0;
src->mType = ET_NONE;
src->mSize = 0;
src->mBits = 0;
src->mChannel = (uint)intVal;
src->mSkip = 0;
}
else
{
TrIndication(tr, &line, &col);
if(!TrReadIdent(tr, MAX_IDENT_LEN, ident))
return 0;
src->mType = MatchElementType(ident);
if(src->mType == ET_NONE)
{
TrErrorAt(tr, line, col, "Expected a source element type.\n");
return 0;
}
if(src->mFormat == SF_BIN_LE || src->mFormat == SF_BIN_BE)
{
if(!TrReadOperator(tr, ","))
return 0;
if(src->mType == ET_INT)
{
if(!TrReadInt(tr, MIN_BIN_SIZE, MAX_BIN_SIZE, &intVal))
return 0;
src->mSize = (uint)intVal;
if(!TrIsOperator(tr, ","))
src->mBits = (int)(8*src->mSize);
else
{
TrReadOperator(tr, ",");
TrIndication(tr, &line, &col);
if(!TrReadInt(tr, -2147483647-1, 2147483647, &intVal))
return 0;
if(abs(intVal) < MIN_BIN_BITS || (uint)abs(intVal) > (8*src->mSize))
{
TrErrorAt(tr, line, col, "Expected a value of (+/-) %d to %d.\n", MIN_BIN_BITS, 8*src->mSize);
return 0;
}
src->mBits = intVal;
}
}
else
{
TrIndication(tr, &line, &col);
if(!TrReadInt(tr, -2147483647-1, 2147483647, &intVal))
return 0;
if(intVal != 4 && intVal != 8)
{
TrErrorAt(tr, line, col, "Expected a value of 4 or 8.\n");
return 0;
}
src->mSize = (uint)intVal;
src->mBits = 0;
}
}
else if(src->mFormat == SF_ASCII && src->mType == ET_INT)
{
if(!TrReadOperator(tr, ","))
return 0;
if(!TrReadInt(tr, MIN_ASCII_BITS, MAX_ASCII_BITS, &intVal))
return 0;
src->mSize = 0;
src->mBits = intVal;
}
else
{
src->mSize = 0;
src->mBits = 0;
}
if(!TrIsOperator(tr, ";"))
src->mSkip = 0;
else
{
TrReadOperator(tr, ";");
if(!TrReadInt(tr, 0, 0x7FFFFFFF, &intVal))
return 0;
src->mSkip = (uint)intVal;
}
}
if(!TrReadOperator(tr, ")"))
return 0;
if(TrIsOperator(tr, "@"))
{
TrReadOperator(tr, "@");
if(!TrReadInt(tr, 0, 0x7FFFFFFF, &intVal))
return 0;
src->mOffset = (uint)intVal;
}
else
src->mOffset = 0;
if(!TrReadOperator(tr, ":"))
return 0;
if(!TrReadString(tr, MAX_PATH_LEN, src->mPath))
return 0;
return 1;
}
// Match the target ear (index) from a given identifier.
static int MatchTargetEar(const char *ident)
{
if(strcasecmp(ident, "left") == 0)
return 0;
if(strcasecmp(ident, "right") == 0)
return 1;
return -1;
}
// Process the list of sources in the data set definition.
static int ProcessSources(const HeadModelT model, TokenReaderT *tr, HrirDataT *hData)
{
uint channels = (hData->mChannelType == CT_STEREO) ? 2 : 1;
double *hrirs = CreateDoubles(channels * hData->mIrCount * hData->mIrSize);
double *hrir = CreateDoubles(hData->mIrPoints);
uint line, col, fi, ei, ai, ti;
int count;
printf("Loading sources...");
fflush(stdout);
count = 0;
while(TrIsOperator(tr, "["))
{
double factor[2] = { 1.0, 1.0 };
TrIndication(tr, &line, &col);
TrReadOperator(tr, "[");
if(!ReadIndexTriplet(tr, hData, &fi, &ei, &ai))
goto error;
if(!TrReadOperator(tr, "]"))
goto error;
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
if(azd->mIrs[0] != NULL)
{
TrErrorAt(tr, line, col, "Redefinition of source.\n");
goto error;
}
if(!TrReadOperator(tr, "="))
goto error;
for(;;)
{
SourceRefT src;
uint ti = 0;
if(!ReadSourceRef(tr, &src))
goto error;
// TODO: Would be nice to display 'x of y files', but that would
// require preparing the source refs first to get a total count
// before loading them.
++count;
printf("\rLoading sources... %d file%s", count, (count==1)?"":"s");
fflush(stdout);
if(!LoadSource(&src, hData->mIrRate, hData->mIrPoints, hrir))
goto error;
if(hData->mChannelType == CT_STEREO)
{
char ident[MAX_IDENT_LEN+1];
if(!TrReadIdent(tr, MAX_IDENT_LEN, ident))
goto error;
ti = MatchTargetEar(ident);
if((int)ti < 0)
{
TrErrorAt(tr, line, col, "Expected a target ear.\n");
goto error;
}
}
azd->mIrs[ti] = &hrirs[hData->mIrSize * (ti * hData->mIrCount + azd->mIndex)];
if(model == HM_DATASET)
azd->mDelays[ti] = AverageHrirOnset(hData->mIrRate, hData->mIrPoints, hrir, 1.0 / factor[ti], azd->mDelays[ti]);
AverageHrirMagnitude(hData->mIrPoints, hData->mFftSize, hrir, 1.0 / factor[ti], azd->mIrs[ti]);
factor[ti] += 1.0;
if(!TrIsOperator(tr, "+"))
break;
TrReadOperator(tr, "+");
}
if(hData->mChannelType == CT_STEREO)
{
if(azd->mIrs[0] == NULL)
{
TrErrorAt(tr, line, col, "Missing left ear source reference(s).\n");
goto error;
}
else if(azd->mIrs[1] == NULL)
{
TrErrorAt(tr, line, col, "Missing right ear source reference(s).\n");
goto error;
}
}
}
printf("\n");
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
if(azd->mIrs[0] != NULL)
break;
}
if(ai < hData->mFds[fi].mEvs[ei].mAzCount)
break;
}
if(ei >= hData->mFds[fi].mEvCount)
{
TrError(tr, "Missing source references [ %d, *, * ].\n", fi);
goto error;
}
hData->mFds[fi].mEvStart = ei;
for(;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
if(azd->mIrs[0] == NULL)
{
TrError(tr, "Missing source reference [ %d, %d, %d ].\n", fi, ei, ai);
goto error;
}
}
}
}
for(ti = 0;ti < channels;ti++)
{
for(fi = 0;fi < hData->mFdCount;fi++)
{
for(ei = 0;ei < hData->mFds[fi].mEvCount;ei++)
{
for(ai = 0;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
{
HrirAzT *azd = &hData->mFds[fi].mEvs[ei].mAzs[ai];
azd->mIrs[ti] = &hrirs[hData->mIrSize * (ti * hData->mIrCount + azd->mIndex)];
}
}
}
}
if(!TrLoad(tr))
{
free(hrir);
return 1;
}
TrError(tr, "Errant data at end of source list.\n");
error:
free(hrir);
return 0;
}
/* Parse the data set definition and process the source data, storing the
* resulting data set as desired. If the input name is NULL it will read
* from standard input.
*/
static int ProcessDefinition(const char *inName, const uint outRate, const uint fftSize, const int equalize, const int surface, const double limit, const uint truncSize, const HeadModelT model, const double radius, const char *outName)
{
char rateStr[8+1], expName[MAX_PATH_LEN];
TokenReaderT tr;
HrirDataT hData;
FILE *fp;
int ret;
ResetHrirData(&hData);
fprintf(stdout, "Reading HRIR definition from %s...\n", inName?inName:"stdin");
if(inName != NULL)
{
fp = fopen(inName, "r");
if(fp == NULL)
{
fprintf(stderr, "Error: Could not open definition file '%s'\n", inName);
return 0;
}
TrSetup(fp, inName, &tr);
}
else
{
fp = stdin;
TrSetup(fp, "<stdin>", &tr);
}
if(!ProcessMetrics(&tr, fftSize, truncSize, &hData))
{
if(inName != NULL)
fclose(fp);
return 0;
}
if(!ProcessSources(model, &tr, &hData))
{
FreeHrirData(&hData);
if(inName != NULL)
fclose(fp);
return 0;
}
if(fp != stdin)
fclose(fp);
if(equalize)
{
uint c = (hData.mChannelType == CT_STEREO) ? 2 : 1;
uint m = 1 + hData.mFftSize / 2;
double *dfa = CreateDoubles(c * m);
fprintf(stdout, "Calculating diffuse-field average...\n");
CalculateDiffuseFieldAverage(&hData, c, m, surface, limit, dfa);
fprintf(stdout, "Performing diffuse-field equalization...\n");
DiffuseFieldEqualize(c, m, dfa, &hData);
free(dfa);
}
fprintf(stdout, "Performing minimum phase reconstruction...\n");
ReconstructHrirs(&hData);
if(outRate != 0 && outRate != hData.mIrRate)
{
fprintf(stdout, "Resampling HRIRs...\n");
ResampleHrirs(outRate, &hData);
}
fprintf(stdout, "Truncating minimum-phase HRIRs...\n");
hData.mIrPoints = truncSize;
fprintf(stdout, "Synthesizing missing elevations...\n");
if(model == HM_DATASET)
SynthesizeOnsets(&hData);
SynthesizeHrirs(&hData);
fprintf(stdout, "Normalizing final HRIRs...\n");
NormalizeHrirs(&hData);
fprintf(stdout, "Calculating impulse delays...\n");
CalculateHrtds(model, (radius > DEFAULT_CUSTOM_RADIUS) ? radius : hData.mRadius, &hData);
snprintf(rateStr, 8, "%u", hData.mIrRate);
StrSubst(outName, "%r", rateStr, MAX_PATH_LEN, expName);
fprintf(stdout, "Creating MHR data set %s...\n", expName);
ret = StoreMhr(&hData, expName);
FreeHrirData(&hData);
return ret;
}
static void PrintHelp(const char *argv0, FILE *ofile)
{
fprintf(ofile, "Usage: %s [<option>...]\n\n", argv0);
fprintf(ofile, "Options:\n");
fprintf(ofile, " -m Ignored for compatibility.\n");
fprintf(ofile, " -r <rate> Change the data set sample rate to the specified value and\n");
fprintf(ofile, " resample the HRIRs accordingly.\n");
fprintf(ofile, " -f <points> Override the FFT window size (default: %u).\n", DEFAULT_FFTSIZE);
fprintf(ofile, " -e {on|off} Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE ? "on" : "off"));
fprintf(ofile, " -s {on|off} Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE ? "on" : "off"));
fprintf(ofile, " -l {<dB>|none} Specify a limit to the magnitude range of the diffuse-field\n");
fprintf(ofile, " average (default: %.2f).\n", DEFAULT_LIMIT);
fprintf(ofile, " -w <points> Specify the size of the truncation window that's applied\n");
fprintf(ofile, " after minimum-phase reconstruction (default: %u).\n", DEFAULT_TRUNCSIZE);
fprintf(ofile, " -d {dataset| Specify the model used for calculating the head-delay timing\n");
fprintf(ofile, " sphere} values (default: %s).\n", ((DEFAULT_HEAD_MODEL == HM_DATASET) ? "dataset" : "sphere"));
fprintf(ofile, " -c <size> Use a customized head radius measured ear-to-ear in meters.\n");
fprintf(ofile, " -i <filename> Specify an HRIR definition file to use (defaults to stdin).\n");
fprintf(ofile, " -o <filename> Specify an output file. Use of '%%r' will be substituted with\n");
fprintf(ofile, " the data set sample rate.\n");
}
// Standard command line dispatch.
int main(int argc, char *argv[])
{
const char *inName = NULL, *outName = NULL;
uint outRate, fftSize;
int equalize, surface;
char *end = NULL;
HeadModelT model;
uint truncSize;
double radius;
double limit;
int opt;
GET_UNICODE_ARGS(&argc, &argv);
if(argc < 2)
{
fprintf(stdout, "HRTF Processing and Composition Utility\n\n");
PrintHelp(argv[0], stdout);
exit(EXIT_SUCCESS);
}
outName = "./oalsoft_hrtf_%r.mhr";
outRate = 0;
fftSize = 0;
equalize = DEFAULT_EQUALIZE;
surface = DEFAULT_SURFACE;
limit = DEFAULT_LIMIT;
truncSize = DEFAULT_TRUNCSIZE;
model = DEFAULT_HEAD_MODEL;
radius = DEFAULT_CUSTOM_RADIUS;
while((opt=getopt(argc, argv, "mr:f:e:s:l:w:d:c:e:i:o:h")) != -1)
{
switch(opt)
{
case 'm':
fprintf(stderr, "Ignoring unused command '-m'.\n");
break;
case 'r':
outRate = strtoul(optarg, &end, 10);
if(end[0] != '\0' || outRate < MIN_RATE || outRate > MAX_RATE)
{
fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected between %u to %u.\n", optarg, opt, MIN_RATE, MAX_RATE);
exit(EXIT_FAILURE);
}
break;
case 'f':
fftSize = strtoul(optarg, &end, 10);
if(end[0] != '\0' || (fftSize&(fftSize-1)) || fftSize < MIN_FFTSIZE || fftSize > MAX_FFTSIZE)
{
fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected a power-of-two between %u to %u.\n", optarg, opt, MIN_FFTSIZE, MAX_FFTSIZE);
exit(EXIT_FAILURE);
}
break;
case 'e':
if(strcmp(optarg, "on") == 0)
equalize = 1;
else if(strcmp(optarg, "off") == 0)
equalize = 0;
else
{
fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected on or off.\n", optarg, opt);
exit(EXIT_FAILURE);
}
break;
case 's':
if(strcmp(optarg, "on") == 0)
surface = 1;
else if(strcmp(optarg, "off") == 0)
surface = 0;
else
{
fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected on or off.\n", optarg, opt);
exit(EXIT_FAILURE);
}
break;
case 'l':
if(strcmp(optarg, "none") == 0)
limit = 0.0;
else
{
limit = strtod(optarg, &end);
if(end[0] != '\0' || limit < MIN_LIMIT || limit > MAX_LIMIT)
{
fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected between %.0f to %.0f.\n", optarg, opt, MIN_LIMIT, MAX_LIMIT);
exit(EXIT_FAILURE);
}
}
break;
case 'w':
truncSize = strtoul(optarg, &end, 10);
if(end[0] != '\0' || truncSize < MIN_TRUNCSIZE || truncSize > MAX_TRUNCSIZE || (truncSize%MOD_TRUNCSIZE))
{
fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected multiple of %u between %u to %u.\n", optarg, opt, MOD_TRUNCSIZE, MIN_TRUNCSIZE, MAX_TRUNCSIZE);
exit(EXIT_FAILURE);
}
break;
case 'd':
if(strcmp(optarg, "dataset") == 0)
model = HM_DATASET;
else if(strcmp(optarg, "sphere") == 0)
model = HM_SPHERE;
else
{
fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected dataset or sphere.\n", optarg, opt);
exit(EXIT_FAILURE);
}
break;
case 'c':
radius = strtod(optarg, &end);
if(end[0] != '\0' || radius < MIN_CUSTOM_RADIUS || radius > MAX_CUSTOM_RADIUS)
{
fprintf(stderr, "Error: Got unexpected value \"%s\" for option -%c, expected between %.2f to %.2f.\n", optarg, opt, MIN_CUSTOM_RADIUS, MAX_CUSTOM_RADIUS);
exit(EXIT_FAILURE);
}
break;
case 'i':
inName = optarg;
break;
case 'o':
outName = optarg;
break;
case 'h':
PrintHelp(argv[0], stdout);
exit(EXIT_SUCCESS);
default: /* '?' */
PrintHelp(argv[0], stderr);
exit(EXIT_FAILURE);
}
}
if(!ProcessDefinition(inName, outRate, fftSize, equalize, surface, limit,
truncSize, model, radius, outName))
return -1;
fprintf(stdout, "Operation completed.\n");
return EXIT_SUCCESS;
}
Reference in New Issue View Git Blame Copy Permalink