Bring some ATH related code from musepack project

Motivations:
- We need some psychoacoustic model to apply it
in the AT3P during GHA. The idea is try to use ATH
along with some tonality measurement to determine
whether GHA gives profit or just penalty due to
non effitient GHA bit encoding.

- We can try to improve AT1 and AT3 quality applying ATH

1 files changed, 96 insertions, 2 deletions

diff --git a/src/atrac/atrac_psy_common.cpp b/src/atrac/atrac_psy_common.cpp
index 257b8f3..089bf47 100644
--- a/src/atrac/atrac_psy_common.cpp
+++ b/src/atrac/atrac_psy_common.cpp
@@ -18,10 +18,88 @@
 
 #include "atrac_psy_common.h"
 
+
+////////////////////////////////////////////////////////////////////////////////
+namespace {
+
+/*
+ * Borrowed from Musepack.
+ * https://www.musepack.net/
+ */
+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));
+}
+
+} // namespace
+////////////////////////////////////////////////////////////////////////////////
+
 namespace NAtracDEnc {
 
+using std::vector;
+
 //returns 1 for tone-like, 0 - noise-like
-TFloat AnalizeScaleFactorSpread(const std::vector<TScaledBlock>& scaledBlocks) {
+TFloat AnalizeScaleFactorSpread(const vector<TScaledBlock>& scaledBlocks)
+{
     TFloat s = 0.0;
     for (size_t i = 0; i < scaledBlocks.size(); ++i) {
         s += scaledBlocks[i].ScaleFactorIndex;
@@ -41,4 +119,20 @@ TFloat AnalizeScaleFactorSpread(const std::vector<TScaledBlock>& scaledBlocks) {
     return sigma/14.0;
 }
 
-} //namespace NAtracDEnc
+vector<float> CalcATH(int len, int sampleRate)
+{
+    vector<float> res(len);
+    float mf = (float)sampleRate / 2000.0;
+    for (size_t i = 0; i < res.size(); i++) {
+        const float f = (float)(i+1) * mf / len;   // Frequency in kHz
+        float trh = ATHformula_Frank(1.e3 * f) - 100;
+
+	// tmp -= f * f * (int)(EarModelFlag % 100 - 50) * 0.0015;  // 00: +30 dB, 100: -30 dB  @20 kHz
+        trh -= f * f * 0.015;
+
+	res[i] = trh;
+    }
+    return res;
+}
+
+} // namespace NAtracDEnc