ISO 532 B / DIN 45 631 c++ codes (=?iso-8859-1?Q?Hannes_L=F6schke?= )

Subject: ISO 532 B / DIN 45 631 c++ codes From: =?iso-8859-1?Q?Hannes_L=F6schke?= <hannes_loeschke@xxxxxxxx> Date: Tue, 5 Jun 2007 06:57:24 +0000 List-Archive:<http://lists.mcgill.ca/scripts/wa.exe?LIST=AUDITORY> --0-1037622900-1181026644=:66844 Content-Type: text/plain; charset=ascii 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; } ___________________________________________________________ Telefonate ohne weitere Kosten vom PC zum PC: http://messenger.yahoo.de --0-1037622900-1181026644=:66844 Content-Type: text/html; charset=ascii <html><head><style type="text/css"></style></head><body><div style="font-family:arial, helvetica, sans-serif;font-size:10pt"><DIV>Dear List,</DIV> <DIV> </DIV> <DIV>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. </DIV> <DIV>The second function is the transcription of the basic code found in ISO 532 B/ DIN 45 631.</DIV> <DIV> </DIV> <DIV>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.</DIV> <DIV> </DIV> <DIV>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:</DIV> <DIV> </DIV> <DIV> // 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; }</DIV></div> <hr size=1>Yahoo! Messenger - <a href=http://de.rd.yahoo.com/evt=39058/*http://de.messenger.yahoo.com >kostenlos* mit Familie und Freunden von PC zu PC telefonieren </a>. </body></html> --0-1037622900-1181026644=:66844--

This message came from the mail archive
http://www.auditory.org/postings/2007/
maintained by:

DAn Ellis <dpwe@ee.columbia.edu>
Electrical Engineering Dept., Columbia University