2749 lines
88 KiB
C
2749 lines
88 KiB
C
/*
|
|
* HRTF utility for producing and demonstrating the process of creating an
|
|
* OpenAL Soft compatible HRIR data set.
|
|
*
|
|
* Copyright (C) 2011-2014 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"
|
|
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <stdarg.h>
|
|
#include <string.h>
|
|
#include <ctype.h>
|
|
#include <math.h>
|
|
#ifdef HAVE_STRINGS_H
|
|
#include <strings.h>
|
|
#endif
|
|
|
|
// Rely (if naively) on OpenAL's header for the types used for serialization.
|
|
#include "AL/al.h"
|
|
#include "AL/alext.h"
|
|
|
|
#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-15)
|
|
|
|
// 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 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 in the definition
|
|
// file.
|
|
#define MIN_DISTANCE (0.5)
|
|
#define MAX_DISTANCE (2.5)
|
|
|
|
// 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 (512)
|
|
#define MAX_FFTSIZE (16384)
|
|
|
|
// 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 (8)
|
|
#define MAX_TRUNCSIZE (128)
|
|
|
|
// 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_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 01.
|
|
#define MHR_FORMAT ("MinPHR01")
|
|
|
|
// Byte order for the serialization routines.
|
|
enum ByteOrderT {
|
|
BO_NONE = 0,
|
|
BO_LITTLE ,
|
|
BO_BIG
|
|
};
|
|
|
|
// Source format for the references listed in the data set definition.
|
|
enum SourceFormatT {
|
|
SF_NONE = 0,
|
|
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.
|
|
};
|
|
|
|
// Element types for the references listed in the data set definition.
|
|
enum ElementTypeT {
|
|
ET_NONE = 0,
|
|
ET_INT , // Integer elements.
|
|
ET_FP // Floating-point elements.
|
|
};
|
|
|
|
// Head model used for calculating the impulse delays.
|
|
enum HeadModelT {
|
|
HM_NONE = 0,
|
|
HM_DATASET , // Measure the onset from the dataset.
|
|
HM_SPHERE // Calculate the onset using a spherical head model.
|
|
};
|
|
|
|
// Desired output format from the command line.
|
|
enum OutputFormatT {
|
|
OF_NONE = 0,
|
|
OF_MHR , // OpenAL Soft MHR data set file.
|
|
OF_TABLE // OpenAL Soft built-in table file (used when compiling).
|
|
};
|
|
|
|
// Unsigned integer type.
|
|
typedef unsigned int uint;
|
|
|
|
// Serialization types. The trailing digit indicates the number of bits.
|
|
typedef ALubyte uint8;
|
|
|
|
typedef ALint int32;
|
|
typedef ALuint uint32;
|
|
typedef ALuint64SOFT uint64;
|
|
|
|
typedef enum ByteOrderT ByteOrderT;
|
|
typedef enum SourceFormatT SourceFormatT;
|
|
typedef enum ElementTypeT ElementTypeT;
|
|
typedef enum HeadModelT HeadModelT;
|
|
typedef enum OutputFormatT OutputFormatT;
|
|
|
|
typedef struct TokenReaderT TokenReaderT;
|
|
typedef struct SourceRefT SourceRefT;
|
|
typedef struct HrirDataT HrirDataT;
|
|
typedef struct ResamplerT ResamplerT;
|
|
|
|
// Token reader state for parsing the data set definition.
|
|
struct TokenReaderT {
|
|
FILE * mFile;
|
|
const char * mName;
|
|
uint mLine,
|
|
mColumn;
|
|
char mRing [TR_RING_SIZE];
|
|
size_t mIn,
|
|
mOut;
|
|
};
|
|
|
|
// Source reference state used when loading sources.
|
|
struct SourceRefT {
|
|
SourceFormatT mFormat;
|
|
ElementTypeT mType;
|
|
uint mSize;
|
|
int mBits;
|
|
uint mChannel,
|
|
mSkip,
|
|
mOffset;
|
|
char mPath [MAX_PATH_LEN + 1];
|
|
};
|
|
|
|
// The HRIR metrics and data set used when loading, processing, and storing
|
|
// the resulting HRTF.
|
|
struct HrirDataT {
|
|
uint mIrRate,
|
|
mIrCount,
|
|
mIrSize,
|
|
mIrPoints,
|
|
mFftSize,
|
|
mEvCount,
|
|
mEvStart,
|
|
mAzCount [MAX_EV_COUNT],
|
|
mEvOffset [MAX_EV_COUNT];
|
|
double mRadius,
|
|
mDistance,
|
|
* mHrirs,
|
|
* mHrtds,
|
|
mMaxHrtd;
|
|
};
|
|
|
|
// The resampler metrics and FIR filter.
|
|
struct ResamplerT {
|
|
uint mP,
|
|
mQ,
|
|
mM,
|
|
mL;
|
|
double * mF;
|
|
};
|
|
|
|
/* Token reader routines for parsing text files. 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;
|
|
char ch;
|
|
|
|
tr -> mFile = fp;
|
|
name = filename;
|
|
// If a filename was given, store a pointer to the base name.
|
|
if (filename != NULL) {
|
|
while ((ch = (* filename)) != '\0') {
|
|
if ((ch == '/') || (ch == '\\'))
|
|
name = filename + 1;
|
|
filename ++;
|
|
}
|
|
}
|
|
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 != NULL) {
|
|
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 != NULL)
|
|
(* line) = tr -> mLine;
|
|
if (column != NULL)
|
|
(* column) = tr -> mColumn;
|
|
}
|
|
|
|
// 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) {
|
|
TrErrorAt (tr, tr -> mLine, col, "Identifier is too long.\n");
|
|
return (0);
|
|
}
|
|
ident [len] = '\0';
|
|
return (1);
|
|
}
|
|
}
|
|
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);
|
|
}
|
|
|
|
// 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)));
|
|
}
|
|
|
|
// Performs a high-passed triangular probability density function dither from
|
|
// a double to an integer. It assumes the input sample is already scaled.
|
|
static int HpTpdfDither (const double in, int * hpHist) {
|
|
const double PRNG_SCALE = 1.0 / (RAND_MAX + 1.0);
|
|
int prn;
|
|
double out;
|
|
|
|
prn = rand ();
|
|
out = round (in + (PRNG_SCALE * (prn - (* hpHist))));
|
|
(* hpHist) = prn;
|
|
return ((int) out);
|
|
}
|
|
|
|
// Allocates an array of doubles.
|
|
static double *CreateArray(const size_t n)
|
|
{
|
|
double *a;
|
|
|
|
a = calloc(n, sizeof(double));
|
|
if(a == NULL)
|
|
{
|
|
fprintf(stderr, "Error: Out of memory.\n");
|
|
exit(-1);
|
|
}
|
|
return a;
|
|
}
|
|
|
|
// Frees an array of doubles.
|
|
static void DestroyArray(double *a)
|
|
{ free(a); }
|
|
|
|
// Complex number routines. All outputs must be non-NULL.
|
|
|
|
// Magnitude/absolute value.
|
|
static double ComplexAbs (const double r, const double i) {
|
|
return (sqrt ((r * r) + (i * i)));
|
|
}
|
|
|
|
// Multiply.
|
|
static void ComplexMul (const double aR, const double aI, const double bR, const double bI, double * outR, double * outI) {
|
|
(* outR) = (aR * bR) - (aI * bI);
|
|
(* outI) = (aI * bR) + (aR * bI);
|
|
}
|
|
|
|
// Base-e exponent.
|
|
static void ComplexExp (const double inR, const double inI, double * outR, double * outI) {
|
|
double e;
|
|
|
|
e = exp (inR);
|
|
(* outR) = e * cos (inI);
|
|
(* outI) = e * sin (inI);
|
|
}
|
|
|
|
/* 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 double * inR, const double * inI, double * outR, double * outI) {
|
|
uint rk, k, m;
|
|
double tempR, tempI;
|
|
|
|
if ((inR == outR) && (inI == outI)) {
|
|
// Handle in-place arrangement.
|
|
rk = 0;
|
|
for (k = 0; k < n; k ++) {
|
|
if (rk > k) {
|
|
tempR = inR [rk];
|
|
tempI = inI [rk];
|
|
outR [rk] = inR [k];
|
|
outI [rk] = inI [k];
|
|
outR [k] = tempR;
|
|
outI [k] = tempI;
|
|
}
|
|
m = n;
|
|
while (rk & (m >>= 1))
|
|
rk &= ~m;
|
|
rk |= m;
|
|
}
|
|
} else {
|
|
// Handle copy arrangement.
|
|
rk = 0;
|
|
for (k = 0; k < n; k ++) {
|
|
outR [rk] = inR [k];
|
|
outI [rk] = inI [k];
|
|
m = n;
|
|
while (rk & (m >>= 1))
|
|
rk &= ~m;
|
|
rk |= m;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Performs the summation.
|
|
static void FftSummation (const uint n, const double s, double * re, double * im) {
|
|
double pi;
|
|
uint m, m2;
|
|
double vR, vI, wR, wI;
|
|
uint i, k, mk;
|
|
double tR, tI;
|
|
|
|
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))
|
|
vR = sin (0.5 * pi / m);
|
|
vR = -2.0 * vR * vR;
|
|
vI = -sin (pi / m);
|
|
// w = Complex (1.0, 0.0)
|
|
wR = 1.0;
|
|
wI = 0.0;
|
|
for (i = 0; i < m; i ++) {
|
|
for (k = i; k < n; k += m2) {
|
|
mk = k + m;
|
|
// t = ComplexMul (w, out [km2])
|
|
tR = (wR * re [mk]) - (wI * im [mk]);
|
|
tI = (wR * im [mk]) + (wI * re [mk]);
|
|
// out [mk] = ComplexSub (out [k], t)
|
|
re [mk] = re [k] - tR;
|
|
im [mk] = im [k] - tI;
|
|
// out [k] = ComplexAdd (out [k], t)
|
|
re [k] += tR;
|
|
im [k] += tI;
|
|
}
|
|
// t = ComplexMul (v, w)
|
|
tR = (vR * wR) - (vI * wI);
|
|
tI = (vR * wI) + (vI * wR);
|
|
// w = ComplexAdd (w, t)
|
|
wR += tR;
|
|
wI += tI;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Performs a forward FFT.
|
|
static void FftForward (const uint n, const double * inR, const double * inI, double * outR, double * outI) {
|
|
FftArrange (n, inR, inI, outR, outI);
|
|
FftSummation (n, 1.0, outR, outI);
|
|
}
|
|
|
|
// Performs an inverse FFT.
|
|
static void FftInverse (const uint n, const double * inR, const double * inI, double * outR, double * outI) {
|
|
double f;
|
|
uint i;
|
|
|
|
FftArrange (n, inR, inI, outR, outI);
|
|
FftSummation (n, -1.0, outR, outI);
|
|
f = 1.0 / n;
|
|
for (i = 0; i < n; i ++) {
|
|
outR [i] *= f;
|
|
outI [i] *= f;
|
|
}
|
|
}
|
|
|
|
/* Calculate the complex helical sequence (or discrete-time analytical
|
|
* signal) of the given input using the Hilbert transform. Given the
|
|
* negative 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 double * in, double * outR, double * outI) {
|
|
uint i;
|
|
|
|
if (in == outR) {
|
|
// Handle in-place operation.
|
|
for (i = 0; i < n; i ++)
|
|
outI [i] = 0.0;
|
|
} else {
|
|
// Handle copy operation.
|
|
for (i = 0; i < n; i ++) {
|
|
outR [i] = in [i];
|
|
outI [i] = 0.0;
|
|
}
|
|
}
|
|
FftForward (n, outR, outI, outR, outI);
|
|
/* Currently the Fourier routines operate only on point counts that are
|
|
* powers of two. If that changes and n is odd, the following conditional
|
|
* should be: i < (n + 1) / 2.
|
|
*/
|
|
for (i = 1; i < (n / 2); i ++) {
|
|
outR [i] *= 2.0;
|
|
outI [i] *= 2.0;
|
|
}
|
|
// If n is odd, the following increment should be skipped.
|
|
i ++;
|
|
for (; i < n; i ++) {
|
|
outR [i] = 0.0;
|
|
outI [i] = 0.0;
|
|
}
|
|
FftInverse (n, outR, outI, outR, outI);
|
|
}
|
|
|
|
/* 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 double * inR, const double * inI, double * out) {
|
|
const uint m = 1 + (n / 2);
|
|
uint i;
|
|
|
|
for (i = 0; i < m; i ++)
|
|
out [i] = fmax (ComplexAbs (inR [i], inI [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 double limit, const double * in, double * out) {
|
|
const uint m = 1 + (n / 2);
|
|
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, double * outR, double * outI) {
|
|
const uint m = 1 + (n / 2);
|
|
double * mags = NULL;
|
|
uint i;
|
|
double aR, aI;
|
|
|
|
mags = CreateArray (n);
|
|
for (i = 0; i < m; i ++) {
|
|
mags [i] = fmax (in [i], EPSILON);
|
|
outR [i] = -log (mags [i]);
|
|
}
|
|
for (; i < n; i ++) {
|
|
mags [i] = mags [n - i];
|
|
outR [i] = outR [n - i];
|
|
}
|
|
Hilbert (n, outR, outR, outI);
|
|
// Remove any DC offset the filter has.
|
|
outR [0] = 0.0;
|
|
outI [0] = 0.0;
|
|
for (i = 1; i < n; i ++) {
|
|
ComplexExp (0.0, outI [i], & aR, & aI);
|
|
ComplexMul (mags [i], 0.0, aR, aI, & outR [i], & outI [i]);
|
|
}
|
|
DestroyArray (mags);
|
|
}
|
|
|
|
/* 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) {
|
|
double k2;
|
|
|
|
k2 = Clamp (k, -1.0, 1.0);
|
|
if ((k < -1.0) || (k > 1.0))
|
|
return (0.0);
|
|
k2 *= k2;
|
|
return (BesselI_0 (b * sqrt (1.0 - k2)) / BesselI_0 (b));
|
|
}
|
|
|
|
// Calculates the greatest common divisor of a and b.
|
|
static uint Gcd (const uint a, const uint b) {
|
|
uint x, y, z;
|
|
|
|
x = a;
|
|
y = b;
|
|
while (y > 0) {
|
|
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;
|
|
|
|
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));
|
|
else if (rejection >= 21.0)
|
|
return ((0.5842 * pow (rejection - 21.0, 0.4)) +
|
|
(0.07886 * (rejection - 21.0)));
|
|
else
|
|
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) {
|
|
uint gcd, l;
|
|
double cutoff, width, beta;
|
|
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.45 / rs -> mP;
|
|
width = 0.1 / rs -> mP;
|
|
} else {
|
|
cutoff = 0.45 / rs -> mQ;
|
|
width = 0.1 / rs -> mQ;
|
|
}
|
|
// A rejection of -180 dB is used for the stop band.
|
|
l = CalcKaiserOrder (180.0, width) / 2;
|
|
beta = CalcKaiserBeta (180.0);
|
|
rs -> mM = (2 * l) + 1;
|
|
rs -> mL = l;
|
|
rs -> mF = CreateArray (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) {
|
|
DestroyArray (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;
|
|
double * work = NULL;
|
|
uint i;
|
|
double r;
|
|
uint j_f, j_s;
|
|
|
|
if (outN == 0)
|
|
return;
|
|
|
|
// Handle in-place operation.
|
|
if (in == out)
|
|
work = CreateArray (outN);
|
|
else
|
|
work = out;
|
|
// Resample the input.
|
|
for (i = 0; i < outN; i ++) {
|
|
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 (in == out) {
|
|
for (i = 0; i < outN; i ++)
|
|
out [i] = work [i];
|
|
DestroyArray (work);
|
|
}
|
|
}
|
|
|
|
// 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);
|
|
}
|
|
|
|
// 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);
|
|
}
|
|
|
|
/* 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) = (double) 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) {
|
|
int v;
|
|
|
|
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 {
|
|
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) {
|
|
FILE * fp = NULL;
|
|
int result;
|
|
|
|
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);
|
|
}
|
|
|
|
// Calculate the onset time of an HRIR and average it with any existing
|
|
// timing for its elevation and azimuth.
|
|
static void AverageHrirOnset (const double * hrir, const double f, const uint ei, const uint ai, const HrirDataT * hData) {
|
|
double mag;
|
|
uint n, i, j;
|
|
|
|
mag = 0.0;
|
|
n = hData -> mIrPoints;
|
|
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;
|
|
}
|
|
j = hData -> mEvOffset [ei] + ai;
|
|
hData -> mHrtds [j] = Lerp (hData -> mHrtds [j], ((double) i) / hData -> mIrRate, f);
|
|
}
|
|
|
|
// Calculate the magnitude response of an HRIR and average it with any
|
|
// existing responses for its elevation and azimuth.
|
|
static void AverageHrirMagnitude (const double * hrir, const double f, const uint ei, const uint ai, const HrirDataT * hData) {
|
|
double * re = NULL, * im = NULL;
|
|
uint n, m, i, j;
|
|
|
|
n = hData -> mFftSize;
|
|
re = CreateArray (n);
|
|
im = CreateArray (n);
|
|
for (i = 0; i < hData -> mIrPoints; i ++) {
|
|
re [i] = hrir [i];
|
|
im [i] = 0.0;
|
|
}
|
|
for (; i < n; i ++) {
|
|
re [i] = 0.0;
|
|
im [i] = 0.0;
|
|
}
|
|
FftForward (n, re, im, re, im);
|
|
MagnitudeResponse (n, re, im, re);
|
|
m = 1 + (n / 2);
|
|
j = (hData -> mEvOffset [ei] + ai) * hData -> mIrSize;
|
|
for (i = 0; i < m; i ++)
|
|
hData -> mHrirs [j + i] = Lerp (hData -> mHrirs [j + i], re [i], f);
|
|
DestroyArray (im);
|
|
DestroyArray (re);
|
|
}
|
|
|
|
/* 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) {
|
|
uint ei;
|
|
double evs, sum, ev, up_ev, down_ev, solidAngle;
|
|
|
|
evs = 90.0 / (hData -> mEvCount - 1);
|
|
sum = 0.0;
|
|
for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++) {
|
|
// For each elevation, calculate the upper and lower limits of the
|
|
// patch band.
|
|
ev = -90.0 + (ei * 2.0 * evs);
|
|
if (ei < (hData -> mEvCount - 1))
|
|
up_ev = (ev + evs) * M_PI / 180.0;
|
|
else
|
|
up_ev = M_PI / 2.0;
|
|
if (ei > 0)
|
|
down_ev = (ev - evs) * M_PI / 180.0;
|
|
else
|
|
down_ev = -M_PI / 2.0;
|
|
// Calculate the area of the patch band.
|
|
solidAngle = 2.0 * M_PI * (sin (up_ev) - sin (down_ev));
|
|
// Each weight is the area of one patch.
|
|
weights [ei] = solidAngle / hData -> mAzCount [ei];
|
|
// Sum the total surface area covered by the HRIRs.
|
|
sum += solidAngle;
|
|
}
|
|
// Normalize the weights given the total surface coverage.
|
|
for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++)
|
|
weights [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 int weighted, const double limit, double * dfa) {
|
|
double * weights = NULL;
|
|
uint ei, ai, count, step, start, end, m, j, i;
|
|
double weight;
|
|
|
|
weights = CreateArray (hData -> mEvCount);
|
|
if (weighted) {
|
|
// Use coverage weighting to calculate the average.
|
|
CalculateDfWeights (hData, weights);
|
|
} else {
|
|
// If coverage weighting is not used, the weights still need to be
|
|
// averaged by the number of HRIRs.
|
|
count = 0;
|
|
for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++)
|
|
count += hData -> mAzCount [ei];
|
|
for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++)
|
|
weights [ei] = 1.0 / count;
|
|
}
|
|
ei = hData -> mEvStart;
|
|
ai = 0;
|
|
step = hData -> mIrSize;
|
|
start = hData -> mEvOffset [ei] * step;
|
|
end = hData -> mIrCount * step;
|
|
m = 1 + (hData -> mFftSize / 2);
|
|
for (i = 0; i < m; i ++)
|
|
dfa [i] = 0.0;
|
|
for (j = start; j < end; j += step) {
|
|
// Get the weight for this HRIR's contribution.
|
|
weight = weights [ei];
|
|
// Add this HRIR's weighted power average to the total.
|
|
for (i = 0; i < m; i ++)
|
|
dfa [i] += weight * hData -> mHrirs [j + i] * hData -> mHrirs [j + i];
|
|
// Determine the next weight to use.
|
|
ai ++;
|
|
if (ai >= hData -> mAzCount [ei]) {
|
|
ei ++;
|
|
ai = 0;
|
|
}
|
|
}
|
|
// Finish the average calculation and keep it from being too small.
|
|
for (i = 0; i < m; i ++)
|
|
dfa [i] = fmax (sqrt (dfa [i]), EPSILON);
|
|
// Apply a limit to the magnitude range of the diffuse-field average if
|
|
// desired.
|
|
if (limit > 0.0)
|
|
LimitMagnitudeResponse (hData -> mFftSize, limit, dfa, dfa);
|
|
DestroyArray (weights);
|
|
}
|
|
|
|
// Perform diffuse-field equalization on the magnitude responses of the HRIR
|
|
// set using the given average response.
|
|
static void DiffuseFieldEqualize (const double * dfa, const HrirDataT * hData) {
|
|
uint step, start, end, m, j, i;
|
|
|
|
step = hData -> mIrSize;
|
|
start = hData -> mEvOffset [hData -> mEvStart] * step;
|
|
end = hData -> mIrCount * step;
|
|
m = 1 + (hData -> mFftSize / 2);
|
|
for (j = start; j < end; j += step) {
|
|
for (i = 0; i < m; i ++)
|
|
hData -> mHrirs [j + i] /= dfa [i];
|
|
}
|
|
}
|
|
|
|
// Perform minimum-phase reconstruction using the magnitude responses of the
|
|
// HRIR set.
|
|
static void ReconstructHrirs (const HrirDataT * hData) {
|
|
double * re = NULL, * im = NULL;
|
|
uint step, start, end, n, j, i;
|
|
|
|
step = hData -> mIrSize;
|
|
start = hData -> mEvOffset [hData -> mEvStart] * step;
|
|
end = hData -> mIrCount * step;
|
|
n = hData -> mFftSize;
|
|
re = CreateArray (n);
|
|
im = CreateArray (n);
|
|
for (j = start; j < end; j += step) {
|
|
MinimumPhase (n, & hData -> mHrirs [j], re, im);
|
|
FftInverse (n, re, im, re, im);
|
|
for (i = 0; i < hData -> mIrPoints; i ++)
|
|
hData -> mHrirs [j + i] = re [i];
|
|
}
|
|
DestroyArray (im);
|
|
DestroyArray (re);
|
|
}
|
|
|
|
// Resamples the HRIRs for use at the given sampling rate.
|
|
static void ResampleHrirs (const uint rate, HrirDataT * hData) {
|
|
ResamplerT rs;
|
|
uint n, step, start, end, j;
|
|
|
|
ResamplerSetup (& rs, hData -> mIrRate, rate);
|
|
n = hData -> mIrPoints;
|
|
step = hData -> mIrSize;
|
|
start = hData -> mEvOffset [hData -> mEvStart] * step;
|
|
end = hData -> mIrCount * step;
|
|
for (j = start; j < end; j += step)
|
|
ResamplerRun (& rs, n, & hData -> mHrirs [j], n, & hData -> mHrirs [j]);
|
|
ResamplerClear (& rs);
|
|
hData -> mIrRate = rate;
|
|
}
|
|
|
|
/* Given an elevation index and an azimuth, calculate the indices of the two
|
|
* HRIRs that bound the coordinate along with a factor for calculating the
|
|
* continous HRIR using interpolation.
|
|
*/
|
|
static void CalcAzIndices (const HrirDataT * hData, const uint ei, const double az, uint * j0, uint * j1, double * jf) {
|
|
double af;
|
|
uint ai;
|
|
|
|
af = ((2.0 * M_PI) + az) * hData -> mAzCount [ei] / (2.0 * M_PI);
|
|
ai = ((uint) af) % hData -> mAzCount [ei];
|
|
af -= floor (af);
|
|
(* j0) = hData -> mEvOffset [ei] + ai;
|
|
(* j1) = hData -> mEvOffset [ei] + ((ai + 1) % hData -> mAzCount [ei]);
|
|
(* jf) = af;
|
|
}
|
|
|
|
// Synthesize any missing onset timings at the bottom elevations. This just
|
|
// blends between slightly exaggerated known onsets. Not an accurate model.
|
|
static void SynthesizeOnsets (HrirDataT * hData) {
|
|
uint oi, e, a, j0, j1;
|
|
double t, of, jf;
|
|
|
|
oi = hData -> mEvStart;
|
|
t = 0.0;
|
|
for (a = 0; a < hData -> mAzCount [oi]; a ++)
|
|
t += hData -> mHrtds [hData -> mEvOffset [oi] + a];
|
|
hData -> mHrtds [0] = 1.32e-4 + (t / hData -> mAzCount [oi]);
|
|
for (e = 1; e < hData -> mEvStart; e ++) {
|
|
of = ((double) e) / hData -> mEvStart;
|
|
for (a = 0; a < hData -> mAzCount [e]; a ++) {
|
|
CalcAzIndices (hData, oi, a * 2.0 * M_PI / hData -> mAzCount [e], & j0, & j1, & jf);
|
|
hData -> mHrtds [hData -> mEvOffset [e] + a] = Lerp (hData -> mHrtds [0], Lerp (hData -> mHrtds [j0], hData -> mHrtds [j1], jf), of);
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Attempt to synthesize any missing HRIRs at the bottom elevations. 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 oi, a, e, step, n, i, j;
|
|
double of, b;
|
|
uint j0, j1;
|
|
double jf;
|
|
double lp [4], s0, s1;
|
|
|
|
if (hData -> mEvStart <= 0)
|
|
return;
|
|
step = hData -> mIrSize;
|
|
oi = hData -> mEvStart;
|
|
n = hData -> mIrPoints;
|
|
for (i = 0; i < n; i ++)
|
|
hData -> mHrirs [i] = 0.0;
|
|
for (a = 0; a < hData -> mAzCount [oi]; a ++) {
|
|
j = (hData -> mEvOffset [oi] + a) * step;
|
|
for (i = 0; i < n; i ++)
|
|
hData -> mHrirs [i] += hData -> mHrirs [j + i] / hData -> mAzCount [oi];
|
|
}
|
|
for (e = 1; e < hData -> mEvStart; e ++) {
|
|
of = ((double) e) / hData -> mEvStart;
|
|
b = (1.0 - of) * (3.5e-6 * hData -> mIrRate);
|
|
for (a = 0; a < hData -> mAzCount [e]; a ++) {
|
|
j = (hData -> mEvOffset [e] + a) * step;
|
|
CalcAzIndices (hData, oi, a * 2.0 * M_PI / hData -> mAzCount [e], & j0, & j1, & jf);
|
|
j0 *= step;
|
|
j1 *= step;
|
|
lp [0] = 0.0;
|
|
lp [1] = 0.0;
|
|
lp [2] = 0.0;
|
|
lp [3] = 0.0;
|
|
for (i = 0; i < n; i ++) {
|
|
s0 = hData -> mHrirs [i];
|
|
s1 = Lerp (hData -> mHrirs [j0 + i], hData -> mHrirs [j1 + i], jf);
|
|
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 -> mHrirs [j + 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 -> mHrirs [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 -> mHrirs [i] = lp [3];
|
|
}
|
|
hData -> 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 step, end, n, j, i;
|
|
double maxLevel;
|
|
|
|
step = hData -> mIrSize;
|
|
end = hData -> mIrCount * step;
|
|
n = hData -> mIrPoints;
|
|
maxLevel = 0.0;
|
|
for (j = 0; j < end; j += step) {
|
|
for (i = 0; i < n; i ++)
|
|
maxLevel = fmax (fabs (hData -> mHrirs [j + i]), maxLevel);
|
|
}
|
|
maxLevel = 1.01 * maxLevel;
|
|
for (j = 0; j < end; j += step) {
|
|
for (i = 0; i < n; i ++)
|
|
hData -> mHrirs [j + 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) {
|
|
double minHrtd, maxHrtd;
|
|
uint e, a, j;
|
|
double t;
|
|
|
|
minHrtd = 1000.0;
|
|
maxHrtd = -1000.0;
|
|
for (e = 0; e < hData -> mEvCount; e ++) {
|
|
for (a = 0; a < hData -> mAzCount [e]; a ++) {
|
|
j = hData -> mEvOffset [e] + a;
|
|
if (model == HM_DATASET) {
|
|
t = hData -> mHrtds [j] * radius / hData -> mRadius;
|
|
} else {
|
|
t = CalcLTD ((-90.0 + (e * 180.0 / (hData -> mEvCount - 1))) * M_PI / 180.0,
|
|
(a * 360.0 / hData -> mAzCount [e]) * M_PI / 180.0,
|
|
radius, hData -> mDistance);
|
|
}
|
|
hData -> mHrtds [j] = t;
|
|
maxHrtd = fmax (t, maxHrtd);
|
|
minHrtd = fmin (t, minHrtd);
|
|
}
|
|
}
|
|
maxHrtd -= minHrtd;
|
|
for (j = 0; j < hData -> mIrCount; j ++)
|
|
hData -> mHrtds [j] -= minHrtd;
|
|
hData -> mMaxHrtd = maxHrtd;
|
|
}
|
|
|
|
// Store the OpenAL Soft HRTF data set.
|
|
static int StoreMhr (const HrirDataT * hData, const char * filename) {
|
|
FILE * fp = NULL;
|
|
uint e, step, end, n, j, i;
|
|
int hpHist, v;
|
|
|
|
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 -> mIrPoints, fp, filename))
|
|
return (0);
|
|
if (! WriteBin4 (BO_LITTLE, 1, (uint32) hData -> mEvCount, fp, filename))
|
|
return (0);
|
|
for (e = 0; e < hData -> mEvCount; e ++) {
|
|
if (! WriteBin4 (BO_LITTLE, 1, (uint32) hData -> mAzCount [e], fp, filename))
|
|
return (0);
|
|
}
|
|
step = hData -> mIrSize;
|
|
end = hData -> mIrCount * step;
|
|
n = hData -> mIrPoints;
|
|
srand (0x31DF840C);
|
|
for (j = 0; j < end; j += step) {
|
|
hpHist = 0;
|
|
for (i = 0; i < n; i ++) {
|
|
v = HpTpdfDither (32767.0 * hData -> mHrirs [j + i], & hpHist);
|
|
if (! WriteBin4 (BO_LITTLE, 2, (uint32) v, fp, filename))
|
|
return (0);
|
|
}
|
|
}
|
|
for (j = 0; j < hData -> mIrCount; j ++) {
|
|
v = (int) fmin (round (hData -> mIrRate * hData -> mHrtds [j]), MAX_HRTD);
|
|
if (! WriteBin4 (BO_LITTLE, 1, (uint32) v, fp, filename))
|
|
return (0);
|
|
}
|
|
fclose (fp);
|
|
return (1);
|
|
}
|
|
|
|
// Store the OpenAL Soft built-in table.
|
|
static int StoreTable (const HrirDataT * hData, const char * filename) {
|
|
FILE * fp = NULL;
|
|
uint step, end, n, j, i;
|
|
int hpHist, v;
|
|
char text [128 + 1];
|
|
|
|
if ((fp = fopen (filename, "wb")) == NULL) {
|
|
fprintf (stderr, "Error: Could not open table file '%s'.\n", filename);
|
|
return (0);
|
|
}
|
|
snprintf (text, 128, "/* Elevation metrics */\n"
|
|
"static const ALubyte defaultAzCount[%u] = { ", hData -> mEvCount);
|
|
if (! WriteAscii (text, fp, filename))
|
|
return (0);
|
|
for (i = 0; i < hData -> mEvCount; i ++) {
|
|
snprintf (text, 128, "%u, ", hData -> mAzCount [i]);
|
|
if (! WriteAscii (text, fp, filename))
|
|
return (0);
|
|
}
|
|
snprintf (text, 128, "};\n"
|
|
"static const ALushort defaultEvOffset[%u] = { ", hData -> mEvCount);
|
|
if (! WriteAscii (text, fp, filename))
|
|
return (0);
|
|
for (i = 0; i < hData -> mEvCount; i ++) {
|
|
snprintf (text, 128, "%u, ", hData -> mEvOffset [i]);
|
|
if (! WriteAscii (text, fp, filename))
|
|
return (0);
|
|
}
|
|
step = hData -> mIrSize;
|
|
end = hData -> mIrCount * step;
|
|
n = hData -> mIrPoints;
|
|
snprintf (text, 128, "};\n\n"
|
|
"/* HRIR Coefficients */\n"
|
|
"static const ALshort defaultCoeffs[%u] =\n{\n", hData -> mIrCount * n);
|
|
if (! WriteAscii (text, fp, filename))
|
|
return (0);
|
|
srand (0x31DF840C);
|
|
for (j = 0; j < end; j += step) {
|
|
if (! WriteAscii (" ", fp, filename))
|
|
return (0);
|
|
hpHist = 0;
|
|
for (i = 0; i < n; i ++) {
|
|
v = HpTpdfDither (32767.0 * hData -> mHrirs [j + i], & hpHist);
|
|
snprintf (text, 128, " %+d,", v);
|
|
if (! WriteAscii (text, fp, filename))
|
|
return (0);
|
|
}
|
|
if (! WriteAscii ("\n", fp, filename))
|
|
return (0);
|
|
}
|
|
snprintf (text, 128, "};\n\n"
|
|
"/* HRIR Delays */\n"
|
|
"static const ALubyte defaultDelays[%u] =\n{\n"
|
|
" ", hData -> mIrCount);
|
|
if (! WriteAscii (text, fp, filename))
|
|
return (0);
|
|
for (j = 0; j < hData -> mIrCount; j ++) {
|
|
v = (int) fmin (round (hData -> mIrRate * hData -> mHrtds [j]), MAX_HRTD);
|
|
snprintf (text, 128, " %d,", v);
|
|
if (! WriteAscii (text, fp, filename))
|
|
return (0);
|
|
}
|
|
if (! WriteAscii ("\n};\n\n"
|
|
"/* Default HRTF Definition */\n", fp, filename))
|
|
return (0);
|
|
snprintf (text, 128, "static const struct Hrtf DefaultHrtf = {\n"
|
|
" %u, %u, %u, defaultAzCount, defaultEvOffset,\n",
|
|
hData -> mIrRate, hData -> mIrPoints, hData -> mEvCount);
|
|
if (! WriteAscii (text, fp, filename))
|
|
return (0);
|
|
if (! WriteAscii (" defaultCoeffs, defaultDelays, NULL\n"
|
|
"};\n", fp, filename))
|
|
return (0);
|
|
fclose (fp);
|
|
return (1);
|
|
}
|
|
|
|
// 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) {
|
|
char ident [MAX_IDENT_LEN + 1];
|
|
uint line, col;
|
|
int intVal;
|
|
uint points;
|
|
double fpVal;
|
|
int hasRate = 0, hasPoints = 0, hasAzimuths = 0;
|
|
int hasRadius = 0, hasDistance = 0;
|
|
|
|
while (! (hasRate && hasPoints && hasAzimuths && hasRadius && hasDistance)) {
|
|
TrIndication (tr, & line, & col);
|
|
if (! TrReadIdent (tr, MAX_IDENT_LEN, ident))
|
|
return (0);
|
|
if (strcasecmp (ident, "rate") == 0) {
|
|
if (hasRate) {
|
|
TrErrorAt (tr, line, col, "Redefinition of 'rate'.\n");
|
|
return (0);
|
|
}
|
|
if (! TrReadOperator (tr, "="))
|
|
return (0);
|
|
if (! TrReadInt (tr, MIN_RATE, MAX_RATE, & intVal))
|
|
return (0);
|
|
hData -> mIrRate = (uint) intVal;
|
|
hasRate = 1;
|
|
} else if (strcasecmp (ident, "points") == 0) {
|
|
if (hasPoints) {
|
|
TrErrorAt (tr, line, col, "Redefinition of 'points'.\n");
|
|
return (0);
|
|
}
|
|
if (! TrReadOperator (tr, "="))
|
|
return (0);
|
|
TrIndication (tr, & line, & col);
|
|
if (! TrReadInt (tr, MIN_POINTS, MAX_POINTS, & intVal))
|
|
return (0);
|
|
points = (uint) intVal;
|
|
if ((fftSize > 0) && (points > fftSize)) {
|
|
TrErrorAt (tr, line, col, "Value exceeds the overriden FFT size.\n");
|
|
return (0);
|
|
}
|
|
if (points < truncSize) {
|
|
TrErrorAt (tr, line, col, "Value is below the truncation size.\n");
|
|
return (0);
|
|
}
|
|
hData -> mIrPoints = points;
|
|
hData -> mFftSize = fftSize;
|
|
if (fftSize <= 0) {
|
|
points = 1;
|
|
while (points < (4 * hData -> mIrPoints))
|
|
points <<= 1;
|
|
hData -> mFftSize = points;
|
|
hData -> mIrSize = 1 + (points / 2);
|
|
} else {
|
|
hData -> mFftSize = fftSize;
|
|
hData -> mIrSize = 1 + (fftSize / 2);
|
|
if (points > hData -> mIrSize)
|
|
hData -> mIrSize = points;
|
|
}
|
|
hasPoints = 1;
|
|
} else if (strcasecmp (ident, "azimuths") == 0) {
|
|
if (hasAzimuths) {
|
|
TrErrorAt (tr, line, col, "Redefinition of 'azimuths'.\n");
|
|
return (0);
|
|
}
|
|
if (! TrReadOperator (tr, "="))
|
|
return (0);
|
|
hData -> mIrCount = 0;
|
|
hData -> mEvCount = 0;
|
|
hData -> mEvOffset [0] = 0;
|
|
for (;;) {
|
|
if (! TrReadInt (tr, MIN_AZ_COUNT, MAX_AZ_COUNT, & intVal))
|
|
return (0);
|
|
hData -> mAzCount [hData -> mEvCount] = (uint) intVal;
|
|
hData -> mIrCount += (uint) intVal;
|
|
hData -> mEvCount ++;
|
|
if (! TrIsOperator (tr, ","))
|
|
break;
|
|
if (hData -> mEvCount >= MAX_EV_COUNT) {
|
|
TrError (tr, "Exceeded the maximum of %d elevations.\n", MAX_EV_COUNT);
|
|
return (0);
|
|
}
|
|
hData -> mEvOffset [hData -> mEvCount] = hData -> mEvOffset [hData -> mEvCount - 1] + ((uint) intVal);
|
|
TrReadOperator (tr, ",");
|
|
}
|
|
if (hData -> mEvCount < MIN_EV_COUNT) {
|
|
TrErrorAt (tr, line, col, "Did not reach the minimum of %d azimuth counts.\n", MIN_EV_COUNT);
|
|
return (0);
|
|
}
|
|
hasAzimuths = 1;
|
|
} else if (strcasecmp (ident, "radius") == 0) {
|
|
if (hasRadius) {
|
|
TrErrorAt (tr, line, col, "Redefinition of 'radius'.\n");
|
|
return (0);
|
|
}
|
|
if (! TrReadOperator (tr, "="))
|
|
return (0);
|
|
if (! TrReadFloat (tr, MIN_RADIUS, MAX_RADIUS, & fpVal))
|
|
return (0);
|
|
hData -> mRadius = fpVal;
|
|
hasRadius = 1;
|
|
} else if (strcasecmp (ident, "distance") == 0) {
|
|
if (hasDistance) {
|
|
TrErrorAt (tr, line, col, "Redefinition of 'distance'.\n");
|
|
return (0);
|
|
}
|
|
if (! TrReadOperator (tr, "="))
|
|
return (0);
|
|
if (! TrReadFloat (tr, MIN_DISTANCE, MAX_DISTANCE, & fpVal))
|
|
return (0);
|
|
hData -> mDistance = fpVal;
|
|
hasDistance = 1;
|
|
} else {
|
|
TrErrorAt (tr, line, col, "Expected a metric name.\n");
|
|
return (0);
|
|
}
|
|
TrSkipWhitespace (tr);
|
|
}
|
|
return (1);
|
|
}
|
|
|
|
// Parse an index pair from the data set definition.
|
|
static int ReadIndexPair (TokenReaderT * tr, const HrirDataT * hData, uint * ei, uint * ai) {
|
|
int intVal;
|
|
|
|
if (! TrReadInt (tr, 0, (int) hData -> mEvCount, & intVal))
|
|
return (0);
|
|
(* ei) = (uint) intVal;
|
|
if (! TrReadOperator (tr, ","))
|
|
return (0);
|
|
if (! TrReadInt (tr, 0, (int) hData -> mAzCount [(* ei)], & 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);
|
|
else if (strcasecmp (ident, "bin_le") == 0)
|
|
return (SF_BIN_LE);
|
|
else if (strcasecmp (ident, "bin_be") == 0)
|
|
return (SF_BIN_BE);
|
|
else 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);
|
|
else 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) {
|
|
uint line, col;
|
|
char ident [MAX_IDENT_LEN + 1];
|
|
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, ",")) {
|
|
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 {
|
|
src -> mBits = (int) (8 * src -> mSize);
|
|
}
|
|
} 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, ";")) {
|
|
TrReadOperator (tr, ";");
|
|
if (! TrReadInt (tr, 0, 0x7FFFFFFF, & intVal))
|
|
return (0);
|
|
src -> mSkip = (uint) intVal;
|
|
} else {
|
|
src -> mSkip = 0;
|
|
}
|
|
}
|
|
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);
|
|
}
|
|
|
|
// Process the list of sources in the data set definition.
|
|
static int ProcessSources (const HeadModelT model, TokenReaderT * tr, HrirDataT * hData) {
|
|
uint * setCount = NULL, * setFlag = NULL;
|
|
double * hrir = NULL;
|
|
uint line, col, ei, ai;
|
|
SourceRefT src;
|
|
double factor;
|
|
|
|
setCount = (uint *) calloc (hData -> mEvCount, sizeof (uint));
|
|
setFlag = (uint *) calloc (hData -> mIrCount, sizeof (uint));
|
|
hrir = CreateArray (hData -> mIrPoints);
|
|
while (TrIsOperator (tr, "[")) {
|
|
TrIndication (tr, & line, & col);
|
|
TrReadOperator (tr, "[");
|
|
if (ReadIndexPair (tr, hData, & ei, & ai)) {
|
|
if (TrReadOperator (tr, "]")) {
|
|
if (! setFlag [hData -> mEvOffset [ei] + ai]) {
|
|
if (TrReadOperator (tr, "=")) {
|
|
factor = 1.0;
|
|
for (;;) {
|
|
if (ReadSourceRef (tr, & src)) {
|
|
if (LoadSource (& src, hData -> mIrRate, hData -> mIrPoints, hrir)) {
|
|
if (model == HM_DATASET)
|
|
AverageHrirOnset (hrir, 1.0 / factor, ei, ai, hData);
|
|
AverageHrirMagnitude (hrir, 1.0 / factor, ei, ai, hData);
|
|
factor += 1.0;
|
|
if (! TrIsOperator (tr, "+"))
|
|
break;
|
|
TrReadOperator (tr, "+");
|
|
continue;
|
|
}
|
|
}
|
|
DestroyArray (hrir);
|
|
free (setFlag);
|
|
free (setCount);
|
|
return (0);
|
|
}
|
|
setFlag [hData -> mEvOffset [ei] + ai] = 1;
|
|
setCount [ei] ++;
|
|
continue;
|
|
}
|
|
} else {
|
|
TrErrorAt (tr, line, col, "Redefinition of source.\n");
|
|
}
|
|
}
|
|
}
|
|
DestroyArray (hrir);
|
|
free (setFlag);
|
|
free (setCount);
|
|
return (0);
|
|
}
|
|
ei = 0;
|
|
while ((ei < hData -> mEvCount) && (setCount [ei] < 1))
|
|
ei ++;
|
|
if (ei < hData -> mEvCount) {
|
|
hData -> mEvStart = ei;
|
|
while ((ei < hData -> mEvCount) && (setCount [ei] == hData -> mAzCount [ei]))
|
|
ei ++;
|
|
if (ei >= hData -> mEvCount) {
|
|
if (! TrLoad (tr)) {
|
|
DestroyArray (hrir);
|
|
free (setFlag);
|
|
free (setCount);
|
|
return (1);
|
|
} else {
|
|
TrError (tr, "Errant data at end of source list.\n");
|
|
}
|
|
} else {
|
|
TrError (tr, "Missing sources for elevation index %d.\n", ei);
|
|
}
|
|
} else {
|
|
TrError (tr, "Missing source references.\n");
|
|
}
|
|
DestroyArray (hrir);
|
|
free (setFlag);
|
|
free (setCount);
|
|
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 OutputFormatT outFormat, const char * outName) {
|
|
FILE * fp = NULL;
|
|
TokenReaderT tr;
|
|
HrirDataT hData;
|
|
double * dfa = NULL;
|
|
char rateStr [8 + 1], expName [MAX_PATH_LEN];
|
|
|
|
hData . mIrRate = 0;
|
|
hData . mIrPoints = 0;
|
|
hData . mFftSize = 0;
|
|
hData . mIrSize = 0;
|
|
hData . mIrCount = 0;
|
|
hData . mEvCount = 0;
|
|
hData . mRadius = 0;
|
|
hData . mDistance = 0;
|
|
fprintf (stdout, "Reading HRIR definition...\n");
|
|
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);
|
|
}
|
|
hData . mHrirs = CreateArray (hData . mIrCount * hData . mIrSize);
|
|
hData . mHrtds = CreateArray (hData . mIrCount);
|
|
if (! ProcessSources (model, & tr, & hData)) {
|
|
DestroyArray (hData . mHrtds);
|
|
DestroyArray (hData . mHrirs);
|
|
if (inName != NULL)
|
|
fclose (fp);
|
|
return (0);
|
|
}
|
|
if (inName != NULL)
|
|
fclose (fp);
|
|
if (equalize) {
|
|
dfa = CreateArray (1 + (hData . mFftSize / 2));
|
|
fprintf (stdout, "Calculating diffuse-field average...\n");
|
|
CalculateDiffuseFieldAverage (& hData, surface, limit, dfa);
|
|
fprintf (stdout, "Performing diffuse-field equalization...\n");
|
|
DiffuseFieldEqualize (dfa, & hData);
|
|
DestroyArray (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);
|
|
switch (outFormat) {
|
|
case OF_MHR :
|
|
fprintf (stdout, "Creating MHR data set file...\n");
|
|
if (! StoreMhr (& hData, expName))
|
|
return (0);
|
|
break;
|
|
case OF_TABLE :
|
|
fprintf (stderr, "Creating OpenAL Soft table file...\n");
|
|
if (! StoreTable (& hData, expName))
|
|
return (0);
|
|
break;
|
|
default :
|
|
break;
|
|
}
|
|
DestroyArray (hData . mHrtds);
|
|
DestroyArray (hData . mHrirs);
|
|
return (1);
|
|
}
|
|
|
|
// Standard command line dispatch.
|
|
int main (const int argc, const char * argv []) {
|
|
const char * inName = NULL, * outName = NULL;
|
|
OutputFormatT outFormat;
|
|
int argi;
|
|
uint outRate, fftSize;
|
|
int equalize, surface;
|
|
double limit;
|
|
uint truncSize;
|
|
HeadModelT model;
|
|
double radius;
|
|
char * end = NULL;
|
|
|
|
if (argc < 2) {
|
|
fprintf (stderr, "Error: No command specified. See '%s -h' for help.\n", argv [0]);
|
|
return (-1);
|
|
}
|
|
if ((strcmp (argv [1], "--help") == 0) || (strcmp (argv [1], "-h") == 0)) {
|
|
fprintf (stdout, "HRTF Processing and Composition Utility\n\n");
|
|
fprintf (stdout, "Usage: %s <command> [<option>...]\n\n", argv [0]);
|
|
fprintf (stdout, "Commands:\n");
|
|
fprintf (stdout, " -m, --make-mhr Makes an OpenAL Soft compatible HRTF data set.\n");
|
|
fprintf (stdout, " Defaults output to: ./oalsoft_hrtf_%%r.mhr\n");
|
|
fprintf (stdout, " -t, --make-tab Makes the built-in table used when compiling OpenAL Soft.\n");
|
|
fprintf (stdout, " Defaults output to: ./hrtf_tables.inc\n");
|
|
fprintf (stdout, " -h, --help Displays this help information.\n\n");
|
|
fprintf (stdout, "Options:\n");
|
|
fprintf (stdout, " -r=<rate> Change the data set sample rate to the specified value and\n");
|
|
fprintf (stdout, " resample the HRIRs accordingly.\n");
|
|
fprintf (stdout, " -f=<points> Override the FFT window size (defaults to the first power-\n");
|
|
fprintf (stdout, " of-two that fits four times the number of HRIR points).\n");
|
|
fprintf (stdout, " -e={on|off} Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE ? "on" : "off"));
|
|
fprintf (stdout, " -s={on|off} Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE ? "on" : "off"));
|
|
fprintf (stdout, " -l={<dB>|none} Specify a limit to the magnitude range of the diffuse-field\n");
|
|
fprintf (stdout, " average (default: %.2f).\n", DEFAULT_LIMIT);
|
|
fprintf (stdout, " -w=<points> Specify the size of the truncation window that's applied\n");
|
|
fprintf (stdout, " after minimum-phase reconstruction (default: %u).\n", DEFAULT_TRUNCSIZE);
|
|
fprintf (stdout, " -d={dataset| Specify the model used for calculating the head-delay timing\n");
|
|
fprintf (stdout, " sphere} values (default: %s).\n", ((DEFAULT_HEAD_MODEL == HM_DATASET) ? "dataset" : "sphere"));
|
|
fprintf (stdout, " -c=<size> Use a customized head radius measured ear-to-ear in meters.\n");
|
|
fprintf (stdout, " -i=<filename> Specify an HRIR definition file to use (defaults to stdin).\n");
|
|
fprintf (stdout, " -o=<filename> Specify an output file. Overrides command-selected default.\n");
|
|
fprintf (stdout, " Use of '%%r' will be substituted with the data set sample rate.\n");
|
|
return (0);
|
|
}
|
|
if ((strcmp (argv [1], "--make-mhr") == 0) || (strcmp (argv [1], "-m") == 0)) {
|
|
if (argc > 3)
|
|
outName = argv [3];
|
|
else
|
|
outName = "./oalsoft_hrtf_%r.mhr";
|
|
outFormat = OF_MHR;
|
|
} else if ((strcmp (argv [1], "--make-tab") == 0) || (strcmp (argv [1], "-t") == 0)) {
|
|
if (argc > 3)
|
|
outName = argv [3];
|
|
else
|
|
outName = "./hrtf_tables.inc";
|
|
outFormat = OF_TABLE;
|
|
} else {
|
|
fprintf (stderr, "Error: Invalid command '%s'.\n", argv [1]);
|
|
return (-1);
|
|
}
|
|
argi = 2;
|
|
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 (argi < argc) {
|
|
if (strncmp (argv [argi], "-r=", 3) == 0) {
|
|
outRate = strtoul (& argv [argi] [3], & end, 10);
|
|
if ((end [0] != '\0') || (outRate < MIN_RATE) || (outRate > MAX_RATE)) {
|
|
fprintf (stderr, "Error: Expected a value from %u to %u for '-r'.\n", MIN_RATE, MAX_RATE);
|
|
return (-1);
|
|
}
|
|
} else if (strncmp (argv [argi], "-f=", 3) == 0) {
|
|
fftSize = strtoul (& argv [argi] [3], & end, 10);
|
|
if ((end [0] != '\0') || (fftSize & (fftSize - 1)) || (fftSize < MIN_FFTSIZE) || (fftSize > MAX_FFTSIZE)) {
|
|
fprintf (stderr, "Error: Expected a power-of-two value from %u to %u for '-f'.\n", MIN_FFTSIZE, MAX_FFTSIZE);
|
|
return (-1);
|
|
}
|
|
} else if (strncmp (argv [argi], "-e=", 3) == 0) {
|
|
if (strcmp (& argv [argi] [3], "on") == 0) {
|
|
equalize = 1;
|
|
} else if (strcmp (& argv [argi] [3], "off") == 0) {
|
|
equalize = 0;
|
|
} else {
|
|
fprintf (stderr, "Error: Expected 'on' or 'off' for '-e'.\n");
|
|
return (-1);
|
|
}
|
|
} else if (strncmp (argv [argi], "-s=", 3) == 0) {
|
|
if (strcmp (& argv [argi] [3], "on") == 0) {
|
|
surface = 1;
|
|
} else if (strcmp (& argv [argi] [3], "off") == 0) {
|
|
surface = 0;
|
|
} else {
|
|
fprintf (stderr, "Error: Expected 'on' or 'off' for '-s'.\n");
|
|
return (-1);
|
|
}
|
|
} else if (strncmp (argv [argi], "-l=", 3) == 0) {
|
|
if (strcmp (& argv [argi] [3], "none") == 0) {
|
|
limit = 0.0;
|
|
} else {
|
|
limit = strtod (& argv [argi] [3], & end);
|
|
if ((end [0] != '\0') || (limit < MIN_LIMIT) || (limit > MAX_LIMIT)) {
|
|
fprintf (stderr, "Error: Expected 'none' or a value from %.2f to %.2f for '-l'.\n", MIN_LIMIT, MAX_LIMIT);
|
|
return (-1);
|
|
}
|
|
}
|
|
} else if (strncmp (argv [argi], "-w=", 3) == 0) {
|
|
truncSize = strtoul (& argv [argi] [3], & end, 10);
|
|
if ((end [0] != '\0') || (truncSize < MIN_TRUNCSIZE) || (truncSize > MAX_TRUNCSIZE) || (truncSize % MOD_TRUNCSIZE)) {
|
|
fprintf (stderr, "Error: Expected a value from %u to %u in multiples of %u for '-w'.\n", MIN_TRUNCSIZE, MAX_TRUNCSIZE, MOD_TRUNCSIZE);
|
|
return (-1);
|
|
}
|
|
} else if (strncmp (argv [argi], "-d=", 3) == 0) {
|
|
if (strcmp (& argv [argi] [3], "dataset") == 0) {
|
|
model = HM_DATASET;
|
|
} else if (strcmp (& argv [argi] [3], "sphere") == 0) {
|
|
model = HM_SPHERE;
|
|
} else {
|
|
fprintf (stderr, "Error: Expected 'dataset' or 'sphere' for '-d'.\n");
|
|
return (-1);
|
|
}
|
|
} else if (strncmp (argv [argi], "-c=", 3) == 0) {
|
|
radius = strtod (& argv [argi] [3], & end);
|
|
if ((end [0] != '\0') || (radius < MIN_CUSTOM_RADIUS) || (radius > MAX_CUSTOM_RADIUS)) {
|
|
fprintf (stderr, "Error: Expected a value from %.2f to %.2f for '-c'.\n", MIN_CUSTOM_RADIUS, MAX_CUSTOM_RADIUS);
|
|
return (-1);
|
|
}
|
|
} else if (strncmp (argv [argi], "-i=", 3) == 0) {
|
|
inName = & argv [argi] [3];
|
|
} else if (strncmp (argv [argi], "-o=", 3) == 0) {
|
|
outName = & argv [argi] [3];
|
|
} else {
|
|
fprintf (stderr, "Error: Invalid option '%s'.\n", argv [argi]);
|
|
return (-1);
|
|
}
|
|
argi ++;
|
|
}
|
|
if (! ProcessDefinition (inName, outRate, fftSize, equalize, surface, limit, truncSize, model, radius, outFormat, outName))
|
|
return (-1);
|
|
fprintf (stdout, "Operation completed.\n");
|
|
return (0);
|
|
}
|
|
|
|
|
|
|
|
|