|
|
/*
|
|
|
* Musepack audio compression
|
|
|
* Copyright (c) 2005-2009, The Musepack Development Team
|
|
|
* Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf
|
|
|
*
|
|
|
* This library is free software; you can redistribute it and/or
|
|
|
* modify it under the terms of the GNU Lesser General Public
|
|
|
* License as published by the Free Software Foundation; either
|
|
|
* version 2.1 of the License, or (at your option) any later version.
|
|
|
*
|
|
|
* This library 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
|
|
|
* Lesser General Public License for more details.
|
|
|
*
|
|
|
* You should have received a copy of the GNU Lesser General Public
|
|
|
* License along with this library; if not, write to the Free Software
|
|
|
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
|
|
*/
|
|
|
|
|
|
#include "libmpcpsy.h"
|
|
|
#include <mpc/minimax.h>
|
|
|
#include <mpc/mpcmath.h>
|
|
|
|
|
|
/*
|
|
|
* Klemm 1994 and 1997. Experimental data. Sorry, data looks a little bit
|
|
|
* dodderly. Data below 30 Hz is extrapolated from other material, above 18
|
|
|
* kHz the ATH is limited due to the original purpose (too much noise at
|
|
|
* ATH is not good even if it's theoretically inaudible).
|
|
|
*/
|
|
|
|
|
|
// FIXME : move this in a struct, to eliminate multiple declarations
|
|
|
// w_low for long 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
|
|
|
const int wl [PART_LONG] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 38, 41, 44, 47, 50, 54, 58, 62, 67, 72, 78, 84, 91, 98, 106, 115, 124, 134, 145, 157, 170, 184, 199, 216, 234, 254, 276, 301, 329, 360, 396, 437, 485 };
|
|
|
const int wh [PART_LONG] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 37, 40, 43, 46, 49, 53, 57, 61, 66, 71, 77, 83, 90, 97, 105, 114, 123, 133, 144, 156, 169, 183, 198, 215, 233, 253, 275, 300, 328, 359, 395, 436, 484, 511 };
|
|
|
// Width: 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 4 4 4 5 5 6 6 7 7 8 9 9 10 11 12 13 14 15 17 18 20 22 25 28 31 36 41 48 27
|
|
|
|
|
|
// inverse partition-width for long
|
|
|
const float iw [PART_LONG] = { 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/3, 1.f/3, 1.f/3, 1.f/3, 1.f/3, 1.f/4, 1.f/4, 1.f/4, 1.f/5, 1.f/5, 1.f/6, 1.f/6, 1.f/7, 1.f/7, 1.f/8, 1.f/9, 1.f/9, 1.f/10, 1.f/11, 1.f/12, 1.f/13, 1.f/14, 1.f/15, 1.f/17, 1.f/18, 1.f/20, 1.f/22, 1.f/25, 1.f/28, 1.f/31, 1.f/36, 1.f/41, 1.f/48, 1.f/27 };
|
|
|
|
|
|
// w_low for short 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
|
|
|
const int wl_short [PART_SHORT] = { 0, 1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 18, 23, 29, 36, 46, 59, 75, 99 };
|
|
|
const int wh_short [PART_SHORT] = { 0, 1, 2, 3, 5, 6, 7, 9, 12, 14, 18, 23, 29, 36, 46, 58, 75, 99, 127 };
|
|
|
|
|
|
// inverse partition-width for short
|
|
|
const float iw_short [PART_SHORT] = { 1.f, 1.f, 1.f, 1.f, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/3, 1.f/3, 1.f/4, 1.f/6, 1.f/7, 1.f/8, 1.f/11, 1.f/13, 1.f/17, 1.f/25, 1.f/29 };
|
|
|
|
|
|
float MinVal [PART_LONG]; // contains minimum tonality soffsets
|
|
|
float Loudness [PART_LONG]; // weighting factors for loudness calculation
|
|
|
float SPRD [PART_LONG] [PART_LONG]; // tabulated spreading function
|
|
|
float O_MAX;
|
|
|
float O_MIN;
|
|
|
float FAC1;
|
|
|
float FAC2; // constants for offset calculation
|
|
|
float partLtq [PART_LONG]; // threshold in quiet (partitions)
|
|
|
float invLtq [PART_LONG]; // inverse threshold in quiet (partitions, long)
|
|
|
float fftLtq [512]; // threshold in quiet (FFT)
|
|
|
|
|
|
static float
|
|
|
ATHformula_Frank ( float freq )
|
|
|
{
|
|
|
/*
|
|
|
* one value per 100 cent = 1
|
|
|
* semitone = 1/4
|
|
|
* third = 1/12
|
|
|
* octave = 1/40 decade
|
|
|
* rest is linear interpolated, values are currently in millibel rel. 20 Pa
|
|
|
*/
|
|
|
static short tab [] = {
|
|
|
/* 10.0 */ 9669, 9669, 9626, 9512,
|
|
|
/* 12.6 */ 9353, 9113, 8882, 8676,
|
|
|
/* 15.8 */ 8469, 8243, 7997, 7748,
|
|
|
/* 20.0 */ 7492, 7239, 7000, 6762,
|
|
|
/* 25.1 */ 6529, 6302, 6084, 5900,
|
|
|
/* 31.6 */ 5717, 5534, 5351, 5167,
|
|
|
/* 39.8 */ 5004, 4812, 4638, 4466,
|
|
|
/* 50.1 */ 4310, 4173, 4050, 3922,
|
|
|
/* 63.1 */ 3723, 3577, 3451, 3281,
|
|
|
/* 79.4 */ 3132, 3036, 2902, 2760,
|
|
|
/* 100.0 */ 2658, 2591, 2441, 2301,
|
|
|
/* 125.9 */ 2212, 2125, 2018, 1900,
|
|
|
/* 158.5 */ 1770, 1682, 1594, 1512,
|
|
|
/* 199.5 */ 1430, 1341, 1260, 1198,
|
|
|
/* 251.2 */ 1136, 1057, 998, 943,
|
|
|
/* 316.2 */ 887, 846, 744, 712,
|
|
|
/* 398.1 */ 693, 668, 637, 606,
|
|
|
/* 501.2 */ 580, 555, 529, 502,
|
|
|
/* 631.0 */ 475, 448, 422, 398,
|
|
|
/* 794.3 */ 375, 351, 327, 322,
|
|
|
/* 1000.0 */ 312, 301, 291, 268,
|
|
|
/* 1258.9 */ 246, 215, 182, 146,
|
|
|
/* 1584.9 */ 107, 61, 13, -35,
|
|
|
/* 1995.3 */ -96, -156, -179, -235,
|
|
|
/* 2511.9 */ -295, -350, -401, -421,
|
|
|
/* 3162.3 */ -446, -499, -532, -535,
|
|
|
/* 3981.1 */ -513, -476, -431, -313,
|
|
|
/* 5011.9 */ -179, 8, 203, 403,
|
|
|
/* 6309.6 */ 580, 736, 881, 1022,
|
|
|
/* 7943.3 */ 1154, 1251, 1348, 1421,
|
|
|
/* 10000.0 */ 1479, 1399, 1285, 1193,
|
|
|
/* 12589.3 */ 1287, 1519, 1914, 2369,
|
|
|
#if 0
|
|
|
/* 15848.9 */ 3352, 4865, 5942, 6177,
|
|
|
/* 19952.6 */ 6385, 6604, 6833, 7009,
|
|
|
/* 25118.9 */ 7066, 7127, 7191, 7260,
|
|
|
#else
|
|
|
/* 15848.9 */ 3352, 4352, 5352, 6352,
|
|
|
/* 19952.6 */ 7352, 8352, 9352, 9999,
|
|
|
/* 25118.9 */ 9999, 9999, 9999, 9999,
|
|
|
#endif
|
|
|
};
|
|
|
double freq_log;
|
|
|
unsigned index;
|
|
|
|
|
|
if ( freq < 10. ) freq = 10.;
|
|
|
if ( freq > 29853. ) freq = 29853.;
|
|
|
|
|
|
freq_log = 40. * log10 (0.1 * freq); /* 4 steps per third, starting at 10 Hz */
|
|
|
index = (unsigned) freq_log;
|
|
|
return 0.01 * (tab [index] * (1 + index - freq_log) + tab [index+1] * (freq_log - index));
|
|
|
}
|
|
|
|
|
|
|
|
|
/* F U N C T I O N S */
|
|
|
// calculation of the threshold in quiet in FFT-resolution
|
|
|
static void
|
|
|
Ruhehoerschwelle ( PsyModel* m,
|
|
|
unsigned int EarModelFlag,
|
|
|
int Ltq_offset,
|
|
|
int Ltq_max )
|
|
|
{
|
|
|
int n;
|
|
|
int k;
|
|
|
float f;
|
|
|
float erg;
|
|
|
double tmp;
|
|
|
float absLtq [512];
|
|
|
|
|
|
for ( n = 0; n < 512; n++ ) {
|
|
|
f = (float) ( (n+1) * (float)(m->SampleFreq / 2000.) / 512 ); // Frequency in kHz
|
|
|
|
|
|
switch ( EarModelFlag / 100 ) {
|
|
|
case 0: // ISO-threshold in quiet
|
|
|
tmp = 3.64*pow (f,-0.8) - 6.5*exp (-0.6*(f-3.3)*(f-3.3)) + 0.001*pow (f, 4.0);
|
|
|
break;
|
|
|
default:
|
|
|
case 1: // measured threshold in quiet (Nick Berglmeir, Andree Buschmann, Kopfh<66>er)
|
|
|
tmp = 3.00*pow (f,-0.8) - 5.0*exp (-0.1*(f-3.0)*(f-3.0)) + 0.0000015022693846297*pow (f, 6.0) + 10.*exp (-(f-0.1)*(f-0.1));
|
|
|
break;
|
|
|
case 2: // measured threshold in quiet (Filburt, Kopfh<66>er)
|
|
|
tmp = 9.00*pow (f,-0.5) - 15.0*exp (-0.1*(f-4.0)*(f-4.0)) + 0.0341796875*pow (f, 2.5) + 15.*exp (-(f-0.1)*(f-0.1)) - 18;
|
|
|
tmp = mind ( tmp, Ltq_max - 18 );
|
|
|
break;
|
|
|
case 3:
|
|
|
tmp = ATHformula_Frank ( 1.e3 * f );
|
|
|
break;
|
|
|
case 4:
|
|
|
tmp = ATHformula_Frank ( 1.e3 * f );
|
|
|
if ( f > 4.8 ) {
|
|
|
tmp += 3.00*pow (f,-0.8) - 5.0*exp (-0.1*(f-3.0)*(f-3.0)) + 0.0000015022693846297*pow (f, 6.0) + 10.*exp (-(f-0.1)*(f-0.1));
|
|
|
tmp *= 0.5 ;
|
|
|
}
|
|
|
break;
|
|
|
case 5:
|
|
|
tmp = ATHformula_Frank ( 1.e3 * f );
|
|
|
if ( f > 4.8 ) {
|
|
|
tmp = 3.00*pow (f,-0.8) - 5.0*exp (-0.1*(f-3.0)*(f-3.0)) + 0.0000015022693846297*pow (f, 6.0) + 10.*exp (-(f-0.1)*(f-0.1));
|
|
|
}
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
tmp -= f * f * (int)(EarModelFlag % 100 - 50) * 0.0015; // 00: +30 dB, 100: -30 dB @20 kHz
|
|
|
|
|
|
tmp = mind ( tmp, Ltq_max ); // Limit ATH
|
|
|
tmp += Ltq_offset - 23; // Add chosen Offset
|
|
|
fftLtq[n] = absLtq[n] = POW10 ( 0.1 * tmp); // conversion into power
|
|
|
}
|
|
|
|
|
|
// threshold in quiet in partitions (long)
|
|
|
for ( n = 0; n < PART_LONG; n++ ) {
|
|
|
erg = 1.e20f;
|
|
|
for ( k = wl[n]; k <= wh[n]; k++ )
|
|
|
erg = minf (erg, absLtq[k]);
|
|
|
|
|
|
partLtq[n] = erg; // threshold in quiet
|
|
|
invLtq [n] = 1.f / partLtq[n]; // Inverse
|
|
|
}
|
|
|
}
|
|
|
|
|
|
#ifdef _MSC_VER
|
|
|
static double
|
|
|
asinh ( double x )
|
|
|
{
|
|
|
return x >= 0 ? log (sqrt (x*x+1) + x) : -log (sqrt (x*x+1) - x);
|
|
|
}
|
|
|
#endif
|
|
|
|
|
|
|
|
|
static double
|
|
|
Freq2Bark ( double Hz ) // Klemm 2002
|
|
|
{
|
|
|
return 9.97074*asinh (1.1268e-3 * Hz) - 6.25817*asinh (0.197193e-3 * Hz) ;
|
|
|
}
|
|
|
|
|
|
// static double
|
|
|
// Bark2Freq ( double Bark ) // Klemm 2002
|
|
|
// {
|
|
|
// return 956.86 * sinh (0.101561*Bark) + 11.7296 * sinh (0.304992*Bark) + 6.33622e-3*sinh (0.538621*Bark);
|
|
|
// }
|
|
|
|
|
|
static double
|
|
|
LongPart2Bark ( PsyModel* m, int Part )
|
|
|
{
|
|
|
return Freq2Bark ((wl [Part] + wh [Part]) * m->SampleFreq / 2048.);
|
|
|
}
|
|
|
|
|
|
// calculating the table for loudness calculation based on absLtq = ank
|
|
|
static void
|
|
|
Loudness_Tabelle (PsyModel* m)
|
|
|
{
|
|
|
int n;
|
|
|
float midfreq;
|
|
|
float tmp;
|
|
|
|
|
|
// ca. dB(A)
|
|
|
for ( n = 0; n < PART_LONG; n++ ){
|
|
|
midfreq = (wh[n] + wl[n] + 3) * (0.25 * m->SampleFreq / 512); // center frequency in kHz, why +3 ???
|
|
|
tmp = LOG10 (midfreq) - 3.5f; // dB(A)
|
|
|
tmp = -10 * tmp * tmp + 3 - midfreq/3000;
|
|
|
Loudness [n] = POW10 ( 0.1 * tmp ); // conversion into power
|
|
|
}
|
|
|
}
|
|
|
|
|
|
|
|
|
static double
|
|
|
Bass ( float f, float TMN, float NMT, float bass )
|
|
|
{
|
|
|
static unsigned char lfe [11] = { 120, 100, 80, 60, 50, 40, 30, 20, 15, 10, 5 };
|
|
|
int tmp = (int) ( 1024/44100. * f + 0.5 );
|
|
|
|
|
|
switch ( tmp ) {
|
|
|
case 0:
|
|
|
case 1:
|
|
|
case 2:
|
|
|
case 3:
|
|
|
case 4:
|
|
|
case 5:
|
|
|
case 6:
|
|
|
case 7:
|
|
|
case 8:
|
|
|
case 9:
|
|
|
case 10:
|
|
|
return TMN + bass * lfe [tmp];
|
|
|
case 11:
|
|
|
case 12:
|
|
|
case 13:
|
|
|
case 14:
|
|
|
case 15:
|
|
|
case 16:
|
|
|
case 17:
|
|
|
case 18:
|
|
|
return TMN;
|
|
|
case 19:
|
|
|
case 20:
|
|
|
case 21:
|
|
|
case 22:
|
|
|
return TMN*0.75 + NMT*0.25;
|
|
|
case 23:
|
|
|
case 24:
|
|
|
return TMN*0.50 + NMT*0.50;
|
|
|
case 25:
|
|
|
case 26:
|
|
|
return TMN*0.25 + NMT*0.75;
|
|
|
default:
|
|
|
return NMT;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
|
|
|
// calculating the coefficient for utilization of the tonality offset, depending on TMN und NMT
|
|
|
static void
|
|
|
Tonalitaetskoeffizienten ( PsyModel* m )
|
|
|
{
|
|
|
double tmp;
|
|
|
int n;
|
|
|
float bass;
|
|
|
|
|
|
bass = 0.1/8 * m->NMT;
|
|
|
if ( m->MinValChoice <= 2 && bass > 0.1 )
|
|
|
bass = 0.1f;
|
|
|
if ( m->MinValChoice <= 1 )
|
|
|
bass = 0.0f;
|
|
|
|
|
|
// alternative: calculation of the minval-values dependent on TMN and TMN
|
|
|
for ( n = 0; n < PART_LONG; n++ ) {
|
|
|
tmp = Bass ( (wl [n] + wh [n]) / 2048. * m->SampleFreq, m->TMN, m->NMT, bass );
|
|
|
MinVal [n] = POW10 ( -0.1 * tmp ); // conversion into power
|
|
|
}
|
|
|
|
|
|
// calculation of the constants for "tonality offset"
|
|
|
O_MAX = POW10 ( -0.1 * m->TMN );
|
|
|
O_MIN = POW10 ( -0.1 * m->NMT );
|
|
|
FAC1 = POW10 ( -0.1 * (m->NMT - (m->TMN - m->NMT) * 0.229) ) ;
|
|
|
FAC2 = (m->TMN - m->NMT) * (0.99011159 * 0.1);
|
|
|
}
|
|
|
|
|
|
|
|
|
// calculation of the spreading function
|
|
|
static void
|
|
|
Spread ( PsyModel* m )
|
|
|
{
|
|
|
int i;
|
|
|
int j;
|
|
|
float tmpx;
|
|
|
float tmpy;
|
|
|
float tmpz;
|
|
|
float x;
|
|
|
|
|
|
// calculation of the spreading-function for all occuring values
|
|
|
for ( i = 0; i < PART_LONG; i++ ) { // i is masking Partition, Source
|
|
|
for ( j = 0; j < PART_LONG; j++ ) { // j is masking Partition, Target
|
|
|
tmpx = LongPart2Bark (m, j) - LongPart2Bark (m, i);// Difference of the partitions in Bark
|
|
|
tmpy = tmpz = 0.; // tmpz = 0: no dip
|
|
|
|
|
|
if ( tmpx < 0 ) { // downwards (S1)
|
|
|
tmpy = -32.f * tmpx; // 32 dB per Bark, e33 (10)
|
|
|
}
|
|
|
else if ( tmpx > 0 ) { // upwards (S2)
|
|
|
#if 0
|
|
|
x = (wl[i]+wh[i])/2 * (float)(SampleFreq / 2000)/512; // center frequency in kHz ???????
|
|
|
if (i==0) x = 0.5f * (float)(SampleFreq / 2000)/512; // if first spectral line
|
|
|
#else
|
|
|
x = i ? wl[i]+wh[i] : 1;
|
|
|
x *= m->SampleFreq / 1000. / 2048; // center frequency in kHz
|
|
|
#endif
|
|
|
// dB/Bark
|
|
|
tmpy = (22.f + 0.23f / x) * tmpx; // e33 (10)
|
|
|
|
|
|
// dip (up to 6 dB)
|
|
|
tmpz = 8 * minf ( (tmpx-0.5f) * (tmpx-0.5f) - 2 * (tmpx-0.5f), 0.f );
|
|
|
}
|
|
|
|
|
|
// calculate coefficient
|
|
|
SPRD[i][j] = POW10 ( -0.1 * (tmpy+tmpz) ); // [Source] [Target]
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// Normierung e33 (10)
|
|
|
for ( i = 0; i < PART_LONG; i++ ) { // i is masked Partition
|
|
|
float norm = 0.f;
|
|
|
for ( j = 0; j < PART_LONG; j++ ) // j is masking Partition
|
|
|
norm += SPRD [j] [i];
|
|
|
for ( j = 0; j < PART_LONG; j++ ) // j is masking Partition
|
|
|
SPRD [j] [i] /= norm;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// call all initialisation procedures
|
|
|
void
|
|
|
Init_Psychoakustiktabellen ( PsyModel* m )
|
|
|
{
|
|
|
m->Max_Band = (int) ( m->BandWidth * 64. / m->SampleFreq );
|
|
|
if ( m->Max_Band < 1 ) m->Max_Band = 1;
|
|
|
if ( m->Max_Band > 31 ) m->Max_Band = 31;
|
|
|
|
|
|
Tonalitaetskoeffizienten (m);
|
|
|
Ruhehoerschwelle ( m, m->EarModelFlag, m->Ltq_offset, m->Ltq_max );
|
|
|
Loudness_Tabelle (m);
|
|
|
Spread (m);
|
|
|
}
|
|
|
|
|
|
/* end of psy_tab.c */
|
