ISO 532 B / DIN 45 631 c++ codes

Dear List,

on popular demand i'll post my c++ transcription of the basic code for loudness calculation. There's two versions of the calculation scheme. The first one is the transcribed basic code released by ZWICKER as found in Acustica 55. I found there were two glitches in the loudness slopes that I took the freedom to correct.

The second function is the transcription of the basic code found in ISO 532 B/ DIN 45 631.

Both are code snipets from my project so the output datatypes might be a little confusing. Basically input is an array of 1/3rd octave levels ranging from 25Hz to 12.5kHz. the Output is an array of specific loudness values scaled to the bark scale from 0 to 24 bark in 0.1 bark steps.

The code is tested and gives reasonable results but I will not take responsibility for incorrectly copied coefficients so please refer to the original sources before using it. Alreight, here it goes:

// Calculation according to Acustica Vol.55

tModularData* calcZwickerLoudness(tModularData* input, double* N, long firstBin, bool freeField )

{

// Hearing thresholds for first three 1/3rd octave Bands

double lowFrequencyThresholds[] = {63.0, 54.0, 47.0, 40.0, 34.0, 28.0};

double thresholdInQuiet[] = {36.0, 21.0, 12.5, 9.0, 7.3, 6.0, 5.0, 4.40, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0}; //% Threshold due to internal noise

//% Hearing thresholds for the excitation levels (each number corresponds to a critical band 12.5kHz is not included)

//% Attenuation representing transmission between freefield and our hearing system

// AO

double outerEarTransmission[] = {.0, .0, .0, .0, .0, .0, .0, .0, .0, .0, -.5, -1.6, -3.2, -5.4, -5.6, -4.0, -1.5, 2.0, 5.0, 12.0};

// Attenuation due to transmission in the middle ear

// Moore et al disagrees with this being flat for low frequencies

//% Level correction to convert from a free field to a diffuse field (last critical band 12.5kHz is not included)

// DDF

double diffuseFieldCorrection[] = {0.0, 0.0, .5, .9, 1.2, 1.6, 2.3, 2.8, 3.0, 2.0, 0.0, -1.4, -2.0, -1.9, -1.0, .5, 3.0, 4.0, 4.3, 4.0};

// Correction factor because using third octave band levels (rather than critical bands)

// DCB

double criticalBandCorrection[] = {-.25, -.6, -.8, -.8, -.5, 0.0, .5, 1.1, 1.5, 1.7, 1.8, 1.8, 1.7, 1.6, 1.4, 1.2, .8, .5, 0.0, -.5};

// Upper limits of the approximated critical bands

// ZUP

double upperLimitOfCriticalBands[] = {.9, 1.8, 2.8, 3.5, 4.4, 5.4, 6.6, 7.9, 9.2, 10.6, 12.3, 13.8, 15.2, 16.7, 18.1, 19.3, 20.6, 21.8, 22.7, 23.6, 24.0};

//% Range of specific loudness for the determination of the steepness of the upper slopes in the specific loudness

//% - critical band rate pattern (used to plot the correct USL curve)

// RNS

double slopeThresholds[] = {23.5, 19.0, 15.1, 11.9, 9.0, 6.6, 4.6, 3.2, 2.13, 1.36, .82, .43, .21, .08, .03, .00};

//double slopeThresholds[] = {21.5, 18.0, 15.1, 11.5, 9.0, 6.1, 4.4, 3.1, 2.13, 1.36, .82, .42, .30, .22, .15, .10, .035, 0.0};

// This is used to design the right hand slope of the loudness

// USL

double upperSlopes[16][8] =

{{13.0, 8.2, /*5.7*/ 6.5, 5.0, 5.0, 5.0, 5.0, 5.0},

{9.0, 7.5, 6.0, 5.1, 4.5, 4.5, 4.5, 4.5},

{7.8, 6.7, 5.6, 4.9, 4.4, 3.9, 3.9, 3.9},

{6.4, 5.5, 4.7, /*4.1*/4.5, 3.6, 3.2, 3.2, 3.2},

{5.6, 5.0, 4.5, 4.3, 3.5, 2.9, 2.9, 2.9},

{4.2, 3.9, 3.7, 3.3, 2.9, 2.42, 2.42, 2.42},

{3.2, 2.8, 2.5, 2.3, 2.2, 2.2, 2.2, 2.02},

{2.8, 2.1, 1.9, 1.8, 1.7, 1.6, 1.6, 1.41},

{1.6, 1.5, 1.4, 1.3, 1.2, 1.1, 1.1, 1.02},

{1.5, 1.2, .94, .77, .77, .77, .77, .77},

{.72, .66, .61, .54, .54, .54, .54, .54},

{.44, .41, .40, .39, .39, .39, .39, .39},

{.29, .25, .22, .22, .22, .22, .22, .22},

{.15, .13, .13, .13, .13, .13, .13, .13},

{.06, .05, .05, .05, .05, .05, .05, .05},

{.04, .04, .04, .04, .04, .04, .04, .04}};

int i;

int j;

double Xp;

double Le;

double factor;

double dz;

// read in 1/3rd octaves

if (input->type != MDT_DOUBLE || (input->length + firstBin < 28)) return NULL;

double Ti[11];

double Lg[3];

double Kern[21];

double Ni;

double Gi;

tModularData* Ns = (tModularData*)malloc(sizeof(tModularData));

Ns->length = 240;

Ns->type = MDT_DOUBLE;

Ns->lengthBytes = Ns->length * sizeof(double);

Ns->valuesDouble = (double*)calloc(Ns->length, sizeof(double));

// correct bands for the first 3 critical bands

double c80 = 0.64* pow(10.0, 0.025*thresholdInQuiet[0]);

for (i = 0; i < 3; i++)

{

Lg[i] = thresholdInQuiet[i];

factor = 0.064 * pow(10.0, 0.025 * lowFrequencyThresholds[i]);

Ni = factor * (pow(1+ 0.25*pow(10.0, 0.1*(input->valuesDouble[i+firstBin]-lowFrequencyThresholds[i])), 0.25)-1);

Gi = 4*(pow(Ni/c80+1,4) - 1);

Ti[i] = 0;

if (Gi > 0)

{

Xp = log10(Gi) + 0.1 * thresholdInQuiet[0];

Ti[i] = pow(10.0, Xp);

}

for ( i = 3 ; i <= 10 ; i++)

{

Ti[i] = pow(10.0,(0.1 * input->valuesDouble[i+firstBin]));

};

Ti[0] = Ti[0]+Ti[1]+Ti[2]+Ti[3]+Ti[4]+Ti[5];

Ti[1] = Ti[6]+Ti[7]+Ti[8];

Ti[2] = Ti[9]+Ti[10];

if (Ti[0] > 0) Ti[0] = 10.0 * log10(Ti[0]);

else Ti[0] = thresholdInQuiet[0];

if (Ti[1] > 0) Ti[1] = 10.0 * log10(Ti[1]);

else Ti[1] = thresholdInQuiet[1];

if (Ti[2] > 0) Ti[2] = 10.0 * log10(Ti[2]);

else Ti[2] = thresholdInQuiet[2];

// determination of main loudness

for (i = 0; i < 20; i++)

{

Le = input->valuesDouble[i+firstBin+8]; // excitation level = level of 1/3rd octave

if (i < 3) // first 11 bands where combined to 3 bands stored

Le = Ti[i]; // in Ti[] so we take them instead of 1/3rd octaves

Le = Le - outerEarTransmission[i]; // apply correction for outer ear transmission

Kern[i] = 0; // initialize loudness kernels

if (!freeField) // apply correction for diffuse field

Le = Le + diffuseFieldCorrection[i];

if (Le <= thresholdInQuiet[i]) // kernel in critical band stays 0

continue; // if excitation level is equal or

// below threshold in quiet

Le = Le - criticalBandCorrection[i];// correction for using 1/3rd octave instead of

// true critical bands

// calculate kernel excitation

factor = 0.064 * pow(10.0, 0.025 * thresholdInQuiet[i]);

Kern[i] = factor * (pow(1 + 0.25 * pow(10.0, 0.1 * (Le - thresholdInQuiet[i])), 0.25)-1);

}

Kern[20] = 0;

// initialize counter variables

*N = 0.2;

double z = 0.1; // position of calculation

double N1 = 0.0; // specific loudness on lower edge of segment

double N2 = 0.0; // specific loudness on upper edge of segment

double z1 = 0.0; // lower edge of section

double z2 = 0.0; // upper edge of section

//double Ns[240];

j = 15; // indes

int iz = 0; // index for specific loudness in z-domain

int ig; // index for slopeselection depending on z

for (i = 0; i < 21; i++) // for each loudness kernel index

{

ig = i-1; // init index

if (ig > 7)

ig = 7; // limit index to 7

while (z1 < upperLimitOfCriticalBands[i])

{

if (N1 <= Kern[i]) // Kernel will affect loudness

{

if (N1 < Kern[i])

{

j = 0; // select index for applicable slope

while((Kern[i] < slopeThresholds[j]) && j < 16) j++;

}

z2 = upperLimitOfCriticalBands[i];

N2 = Kern[i];

*N = *N + N2 * (z2-z1); // integrate for overall loudness

for ( ; z < z2; z+=0.1)

{ // kernel loudness evokes equal excitation

Ns->valuesDouble[iz] = N2; // over full critical band

iz++;

}

} else

{

N2 = slopeThresholds[j];

if (N2 < Kern[i])

N2 = Kern[i];

dz = (N1-N2) / upperSlopes[j][ig]; // calculate delta z from selected slope

z2 = z1 + dz; // apply delta z

if (z2 > upperLimitOfCriticalBands[i])

{ // upper limit of section ran over band boundary

z2 = upperLimitOfCriticalBands[i]; // correct limits

dz = z2-z1; // correct delta

N2 = N1 - dz*upperSlopes[j][ig]; // apply to temp loudness

}

*N = *N + dz * (N2 + N1) * 0.5; // integrate overall loudness

for ( ; z < z2; z+= 0.1)

{

Ns->valuesDouble[iz] = N1 - (z - z1) * upperSlopes[j][ig];

iz++;

}

if (N2 == slopeThresholds[j]) j++;

if (j > 15) j = 15;

N1 = N2;

z1 = z2;

}

return (tModularData*) Ns;

}/**/

// calculation of loudness according to ISO 532

tModularData* calcISO532Loudness(tModularData* input, double* N, long firstBin, bool freeField )

{

// Ranges of 1/3 Oct bands for correction at low frequencies according to equal loudness contours

// RAP

double lowFreqThresholds[] = {45.0, 55.0, 65.0, 71.0, 80.0, 90.0, 100.0, 120.0};

// Reduction of 1/3 Oct Band levels at low frequencies according to equal loudness contours

// within the eight ranges defined by RAP (DLL)

// DLL

double lowFreqCorrection[8][11] =

{{-32.0, -24.0, -16.0, -10.0, -5.0, 0.0, -7.0, -3.0, 0.0, -2.0, 0},

{-29.0, -22.0, -15.0, -10.0, -4.0, 0.0, -7.0, -2.0, 0.0, -2.0, 0},

{-27.0, -19.0, -14.0, -9.0, -4.0, 0.0, -6.0, -2.0, 0.0, -2.0, 0},

{-25.0, -17.0, -12.0, -9.0, -3.0, 0.0, -5.0, -2.0, 0.0, -2.0, 0},

{-23.0, -16.0, -11.0, -7.0, -3.0, 0.0, -4.0, -1.0, 0.0, -1.0, 0},

{-20.0, -14.0, -10.0, -6.0, -3.0, 0.0, -4.0, -1.0, 0.0, -1.0, 0},

{-18.0, -12.0, -9.0, -6.0, -2.0, 0.0, -3.0, -1.0, 0.0, -1.0, 0},

{-15.0, -10.0, -8.0, -4.0, -2.0, 0.0, -3.0, -1.0, 0.0, -1.0, 0}};

//Critical band level at absolute threshold without taking into account the

//transmission characteristics of the ear

// LTQ

double thresholdInQuiet[] = {30.0, 18.0, 12.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3}; //% Threshold due to internal noise

//% Hearing thresholds for the excitation levels (each number corresponds to a critical band 12.5kHz is not included)

//% Attenuation representing transmission between freefield and our hearing system

// AO

double outerEarTransmission[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -.5, -1.6, -3.2, -5.4, -5.6, -4.0, -1.5, 2.0, 5.0, 12.0};

// Attenuation due to transmission in the middle ear

// Moore et al disagrees with this being flat for low frequencies

//% Level correction to convert from a free field to a diffuse field (last critical band 12.5kHz is not included)

// DDF

double diffuseFieldCorrection[] = {0.0, 0.0, .5, .9, 1.2, 1.6, 2.3, 2.8, 3.0, 2.0, 0.0, -1.4, -2.0, -1.9, -1.0, .5, 3.0, 4.0, 4.3, 4.0};

// Correction factor because using third octave band levels (rather than critical bands)

// DCB

double criticalBandCorrection[] = {-.25, -.6, -.8, -.8, -.5, 0.0, .5, 1.1, 1.5, 1.7, 1.8, 1.8, 1.7, 1.6, 1.4, 1.2, .8, .5, 0.0, -.5};

// Upper limits of the approximated critical bands

// ZUP

double upperLimitOfCriticalBands[] = {.9, 1.8, 2.8, 3.5, 4.4, 5.4, 6.6, 7.9, 9.2, 10.6, 12.3, 13.8, 15.2, 16.7, 18.1, 19.3, 20.6, 21.8, 22.7, 23.6, 24};

//% Range of specific loudness for the determination of the steepness of the upper slopes in the specific loudness

//% - critical band rate pattern (used to plot the correct USL curve)

// RNS

double slopeThresholds[] = {21.5, 18.0, 15.1, 11.5, 9.0, 6.1, 4.4, 3.1, 2.13, 1.36, .82, .42, .30, .22, .15, .10, .035, 0.0};

// This is used to design the right hand slope of the loudness

// USL

double upperSlopes[18][8] =

{{13.0, 8.2, 6.3, 5.5, 5.5, 5.5, 5.5, 5.5},

{9.0, 7.5, 6.0, 5.1, 4.5, 4.5, 4.5, 4.5},

{7.8, 6.7, 5.6, 4.9, 4.4, 3.9, 3.9, 3.9},

{6.2, 5.4, 4.6, 4.0, 3.5, 3.2, 3.2, 3.2},

{4.5, 3.8, 3.6, 3.2, 2.9, 2.7, 2.7, 2.7},

{3.7, 3.0, 2.8, 2.35, 2.2, 2.2, 2.2, 2.2},

{2.9, 2.3, 2.1, 1.9, 1.8, 1.7, 1.7, 1.7},

{2.4, 1.7, 1.5, 1.35, 1.3, 1.3, 1.3, 1.3},

{1.95, 1.45, 1.3, 1.15, 1.1, 1.1, 1.1, 1.1},

{1.5, 1.2, .94, .86, .82, .82, .82, .82},

{.72, .67, .64, .63, .62, .62, .62, .62},

{.59, .53, .51, .50, .42, .42, .42, .42},

{.40, .33, .26, .24, .24, .22, .22, .22},

{.27, .21, .20, .18, .17, .17, .17, .17},

{.16, .15, .14, .12, .11, .11, .11, .11},

{.12, .11, .10, .08, .08, .08, .08, .08},

{.09, .08, .07, .06, .06, .06, .06, .05},

{.06, .05, .03, .02, .02, .02, .02, .02}};

int i;

int j;

double Xp;

double Le;

double factor;

double dz;

// read in 1/3rd octaves

if (input->type != MDT_DOUBLE || (input->length + firstBin < 28)) return NULL;

double Ti[11];

double Kern[21];

tModularData* Ns = (tModularData*)malloc(sizeof(tModularData));

Ns->length = 240;

Ns->type = MDT_DOUBLE;

Ns->lengthBytes = Ns->length * sizeof(double);

Ns->valuesDouble = (double*)calloc(Ns->length, sizeof(double));

// correct bands for the first 3 critical bands

for ( i = 0 ; i <= 10 ; i++)

{

j=1;

while (input->valuesDouble[i+firstBin] > (lowFreqThresholds[j]-lowFreqCorrection[i][j]) && (j < 7))

{

j++;

}

Xp = input->valuesDouble[i+firstBin] + lowFreqCorrection[j][i];

Ti[i] = pow(10.0,(0.1 * Xp));

};

Ti[0] = Ti[0]+Ti[1]+Ti[2]+Ti[3]+Ti[4]+Ti[5];

Ti[1] = Ti[6]+Ti[7]+Ti[8];

Ti[2] = Ti[9]+Ti[10];

if (Ti[0] > 0) Ti[0] = 10.0 * log10(Ti[0]);

else Ti[0] = thresholdInQuiet[0];

if (Ti[1] > 0) Ti[1] = 10.0 * log10(Ti[1]);

else Ti[1] = thresholdInQuiet[1];

if (Ti[2] > 0) Ti[2] = 10.0 * log10(Ti[2]);

else Ti[2] = thresholdInQuiet[2];

// determination of main loudness

for (i = 0; i < 20; i++)

{

Le = input->valuesDouble[i+firstBin+8]; // excitation level = level of 1/3rd octave

if (i < 3) // first 11 bands where combined to 3 bands stored

Le = Ti[i]; // in Ti[] so we take them instead of 1/3rd octaves

Le = Le - outerEarTransmission[i]; // apply correction for outer ear transmission

Kern[i] = 0; // initialize loudness kernels

if (!freeField) // apply correction for diffuse field

Le = Le + diffuseFieldCorrection[i];

if (Le <= thresholdInQuiet[i]) // kernel in critical band stays 0

continue; // if excitation level is equal or

// below threshold in quiet

Le = Le - criticalBandCorrection[i];// correction for using 1/3rd octave instead of

// true critical bands

// calculate kernel excitation

factor = 0.0635 * pow(10.0, 0.025 * thresholdInQuiet[i]);

Kern[i] = factor * (pow(0.75 + 0.25 * pow(10.0, 0.1 * (Le - thresholdInQuiet[i])), 0.25)-1);

if (Kern[i] < 0) Kern [i] = 0;

}

Kern[20] = 0;

// correct first critical Band

double corr = 0.4 + 0.32 * pow(Kern[0], 0.2);

if (corr > 1.0) corr = 1.0;

Kern[0] = Kern[0] * corr;

// initialize counter variables

*N = 0.0;

double z = 0.1;

double N1 = 0.0;

double N2 = 0.0;

double z1 = 0.0;

double z2 = 0.0;

//double Ns[240];

j = 17;

int iz = 0; // index in z-domain for specific Loudness

int ig;

for (i = 0; i < 21; i++) // for each loudness kernel index

{

//upperLimitOfCriticalBands[i] += 0.0001;

ig = i-1; // init index

if (ig > 7) ig = 7; // limit index to 7

while (z1 < upperLimitOfCriticalBands[i])

{

if (N1 <= Kern[i]) // Kernel will affect loudness

{

if (N1 < Kern[i])

{

j = 0; // select index for applicable slope

while((Kern[i] < slopeThresholds[j]) && j < 18) j++;

}

z2 = upperLimitOfCriticalBands[i];

N2 = Kern[i];

*N = *N + N2 * (z2-z1); // integrate for overall loudness

for ( ; z < z2; z+=0.1)

{ // kernel loudness evokes equal excitation

Ns->valuesDouble[iz] = N2; // over full critical band

iz++;

}

} else

{

N2 = slopeThresholds[j];

if (N2 < Kern[i])

N2 = Kern[i];

dz = (N1-N2) / upperSlopes[j][ig]; // calculate delta z from selected slope

z2 = z1 + dz; // apply delta z

if (z2 > upperLimitOfCriticalBands[i])

{ // upper limit of section ran over band boundary

z2 = upperLimitOfCriticalBands[i]; // correct limits

dz = z2-z1; // correct delta

N2 = N1 - dz*upperSlopes[j][ig]; // apply to temp loudness

}

*N = *N + dz * (N2 + N1) * 0.5; // integrate overall loudness

for ( ; z < z2; z+= 0.1)

{

Ns->valuesDouble[iz] = N1 - (z - z1) * upperSlopes[j][ig];

iz++;

}

while (N2 <= slopeThresholds[j] && j < 18) j++;

if (j > 17) j = 17;

N1 = N2;

z1 = z2;

}

return (tModularData*) Ns;

}