[AT3P] PQF implementation

M=16 perfect reconstruction polyphase quadrature analysis filter bank

3 files changed, 224 insertions, 109 deletions

diff --git a/src/atrac/atrac3plus_pqf/ut/atrac3plusdsp.c b/src/atrac/atrac3plus_pqf/ut/atrac3plusdsp.c
index 940cf05..85f712d 100644
--- a/src/atrac/atrac3plus_pqf/ut/atrac3plusdsp.c
+++ b/src/atrac/atrac3plus_pqf/ut/atrac3plusdsp.c
@@ -35,116 +35,14 @@
 
 #include "atrac3plusdsp.h"
 
+#include "../atrac3plus_pqf_data.h"
+
 /* lookup table for fast modulo 23 op required for cyclic buffers of the IPQF */
 static const int mod23_lut[26] = {
     23,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11,
     12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0
 };
 
-/* First half of the 384-tap IPQF filtering coefficients. */
-static const float ipqf_coeffs1[ATRAC3P_PQF_FIR_LEN][16] = {
-    { -5.8336207e-7,    -8.0604229e-7,    -4.2005411e-7,    -4.4400572e-8,
-       3.226247e-8,      3.530856e-8,      1.2660377e-8,     0.000010516783,
-      -0.000011838618,   6.005389e-7,      0.0000014333754,  0.0000023108685,
-       0.0000032569742,  0.0000046192422,  0.0000063894258,  0.0000070302972 },
-    { -0.0000091622824, -0.000010502935,  -0.0000079212787, -0.0000041712024,
-      -0.0000026336629, -0.0000015432918, -5.7168614e-7,     0.0000018111954,
-       0.000023530851,   0.00002780562,    0.000032302323,   0.000036968919,
-       0.000041575615,   0.000045337845,   0.000046043948,   0.000048585582 },
-    { -0.000064464548,  -0.000068306952,  -0.000073081472,  -0.00007612785,
-      -0.000074850752,  -0.000070208509,  -0.000062285151,  -0.000058270442,
-      -0.000056296329,  -0.000049888811,  -0.000035615325,  -0.000018532943,
-       0.0000016657353,  0.00002610587,    0.000053397067,   0.00008079566 },
-    { -0.00054488552,   -0.00052537228,   -0.00049731287,   -0.00045778,
-      -0.00040612387,   -0.00034301577,   -0.00026866337,   -0.00018248901,
-      -0.000084307925,   0.000025081157,   0.00014135583,    0.00026649953,
-       0.00039945057,    0.00053928449,    0.00068422867,    0.00083093712 },
-    { -0.0014771431,    -0.001283227,     -0.0010566821,    -0.00079780724,
-      -0.00050782406,   -0.00018855913,    0.00015771533,    0.00052769453,
-       0.00091862219,    0.001326357,      0.0017469483,     0.0021754825,
-       0.0026067684,     0.0030352892,     0.0034549395,     0.0038591374 },
-    { -0.0022995141,    -0.001443546,     -0.00049266568,    0.00055068987,
-       0.001682895,      0.0028992873,     0.0041943151,     0.0055614738,
-       0.0069935122,     0.0084823566,     0.010018963,      0.011593862,
-       0.013196872,      0.014817309,      0.016444042,      0.018065533 },
-    { -0.034426283,     -0.034281436,     -0.033992987,     -0.033563249,
-      -0.032995768,     -0.032295227,     -0.031467363,     -0.030518902,
-      -0.02945766,      -0.028291954,     -0.027031265,     -0.025685543,
-      -0.024265358,     -0.022781773,     -0.021246184,     -0.019670162 },
-    { -0.0030586775,    -0.0037203205,    -0.0042847847,    -0.0047529764,
-      -0.0051268316,    -0.0054091476,    -0.0056034233,    -0.005714261,
-      -0.0057445862,    -0.0057025906,    -0.0055920109,    -0.0054194843,
-      -0.0051914565,    -0.0049146507,    -0.0045959447,    -0.0042418269 },
-    { -0.0016376863,    -0.0017651899,    -0.0018608454,    -0.0019252141,
-      -0.0019593791,    -0.0019653172,    -0.0019450618,    -0.0018990048,
-      -0.00183808,      -0.0017501717,    -0.0016481078,    -0.0015320742,
-      -0.0014046903,    -0.0012685474,    -0.001125814,     -0.00097943726 },
-    { -0.00055432378,   -0.00055472925,   -0.00054783461,   -0.00053276919,
-      -0.00051135791,   -0.00048466062,   -0.00045358928,   -0.00042499689,
-      -0.00036942671,   -0.0003392619,    -0.00030001783,   -0.00025986304,
-      -0.0002197204,    -0.00018116167,   -0.00014691355,   -0.00011279432 },
-    { -0.000064147389,  -0.00006174868,   -0.000054267788,  -0.000047133824,
-      -0.000042927582,  -0.000039477309,  -0.000036340745,  -0.000029687517,
-      -0.000049787737,  -0.000041577889,  -0.000033864744,  -0.000026534748,
-      -0.000019841305,  -0.000014789486,  -0.000013131184,  -0.0000099198869 },
-    { -0.0000062990207, -0.0000072701259, -0.000011984052,  -0.000017348082,
-      -0.000019907106,  -0.000021348773,  -0.000021961965,  -0.000012203576,
-      -0.000010840992,   4.6299544e-7,     5.2588763e-7,     2.7792686e-7,
-      -2.3649704e-7,    -0.0000010897784, -9.171448e-7,     -5.22682e-7 }
-};
-
-/* Second half of the 384-tap IPQF filtering coefficients. */
-static const float ipqf_coeffs2[ATRAC3P_PQF_FIR_LEN][16] = {
-    {  5.22682e-7,       9.171448e-7,      0.0000010897784,  2.3649704e-7,
-      -2.7792686e-7,    -5.2588763e-7,    -4.6299544e-7,     0.000010840992,
-      -0.000012203576,  -0.000021961965,  -0.000021348773,  -0.000019907106,
-      -0.000017348082,  -0.000011984052,  -0.0000072701259, -0.0000062990207 },
-    {  0.0000099198869,  0.000013131184,   0.000014789486,   0.000019841305,
-       0.000026534748,   0.000033864744,   0.000041577889,   0.000049787737,
-      -0.000029687517,  -0.000036340745,  -0.000039477309,  -0.000042927582,
-      -0.000047133824,  -0.000054267788,  -0.00006174868,   -0.000064147389 },
-    {  0.00011279432,    0.00014691355,    0.00018116167,    0.0002197204,
-       0.00025986304,    0.00030001783,    0.0003392619,     0.00036942671,
-      -0.00042499689,   -0.00045358928,   -0.00048466062,   -0.00051135791,
-      -0.00053276919,   -0.00054783461,   -0.00055472925,   -0.00055432378 },
-    {  0.00097943726,    0.001125814,      0.0012685474,     0.0014046903,
-       0.0015320742,     0.0016481078,     0.0017501717,     0.00183808,
-      -0.0018990048,    -0.0019450618,    -0.0019653172,    -0.0019593791,
-      -0.0019252141,    -0.0018608454,    -0.0017651899,    -0.0016376863 },
-    {  0.0042418269,     0.0045959447,     0.0049146507,     0.0051914565,
-       0.0054194843,     0.0055920109,     0.0057025906,     0.0057445862,
-      -0.005714261,     -0.0056034233,    -0.0054091476,    -0.0051268316,
-      -0.0047529764,    -0.0042847847,    -0.0037203205,    -0.0030586775 },
-    {  0.019670162,      0.021246184,      0.022781773,      0.024265358,
-       0.025685543,      0.027031265,      0.028291954,      0.02945766,
-      -0.030518902,     -0.031467363,     -0.032295227,     -0.032995768,
-      -0.033563249,     -0.033992987,     -0.034281436,     -0.034426283 },
-    { -0.018065533,     -0.016444042,     -0.014817309,     -0.013196872,
-      -0.011593862,     -0.010018963,     -0.0084823566,    -0.0069935122,
-       0.0055614738,     0.0041943151,     0.0028992873,     0.001682895,
-       0.00055068987,   -0.00049266568,   -0.001443546,     -0.0022995141 },
-    { -0.0038591374,    -0.0034549395,    -0.0030352892,    -0.0026067684,
-      -0.0021754825,    -0.0017469483,    -0.001326357,     -0.00091862219,
-       0.00052769453,    0.00015771533,   -0.00018855913,   -0.00050782406,
-      -0.00079780724,   -0.0010566821,    -0.001283227,     -0.0014771431 },
-    { -0.00083093712,   -0.00068422867,   -0.00053928449,   -0.00039945057,
-      -0.00026649953,   -0.00014135583,   -0.000025081157,   0.000084307925,
-      -0.00018248901,   -0.00026866337,   -0.00034301577,   -0.00040612387,
-      -0.00045778,      -0.00049731287,   -0.00052537228,   -0.00054488552 },
-    { -0.00008079566,   -0.000053397067,  -0.00002610587,   -0.0000016657353,
-       0.000018532943,   0.000035615325,   0.000049888811,   0.000056296329,
-      -0.000058270442,  -0.000062285151,  -0.000070208509,  -0.000074850752,
-      -0.00007612785,   -0.000073081472,  -0.000068306952,  -0.000064464548 },
-    { -0.000048585582,  -0.000046043948,  -0.000045337845,  -0.000041575615,
-      -0.000036968919,  -0.000032302323,  -0.00002780562,   -0.000023530851,
-       0.0000018111954, -5.7168614e-7,    -0.0000015432918, -0.0000026336629,
-      -0.0000041712024, -0.0000079212787, -0.000010502935,  -0.0000091622824 },
-    { -0.0000070302972, -0.0000063894258, -0.0000046192422, -0.0000032569742,
-      -0.0000023108685, -0.0000014333754, -6.005389e-7,      0.000011838618,
-       0.000010516783,   1.2660377e-8,     3.530856e-8,      3.226247e-8,
-      -4.4400572e-8,    -4.2005411e-7,    -8.0604229e-7,    -5.8336207e-7 }
-};
-
 // Just for test
 static void dct4(float* out, const float* x, int N, float scale) {
     for (int k = 0; k < N; k++) {
@@ -182,10 +80,10 @@ void ff_atrac3p_ipqf(Atrac3pIPQFChannelCtx *hist, const float *in, float *out)
 
         for (t = 0; t < ATRAC3P_PQF_FIR_LEN; t++) {
             for (i = 0; i < 8; i++) {
-                out[s * 16 + i + 0] += hist->buf1[pos_now][i]  * ipqf_coeffs1[t][i] +
-                                       hist->buf2[pos_next][i] * ipqf_coeffs2[t][i];
-                out[s * 16 + i + 8] += hist->buf1[pos_now][7 - i]  * ipqf_coeffs1[t][i + 8] +
-                                       hist->buf2[pos_next][7 - i] * ipqf_coeffs2[t][i + 8];
+                out[s * 16 + i + 0] += hist->buf1[pos_now][i]  * ff_ipqf_coeffs1[t][i] +
+                                       hist->buf2[pos_next][i] * ff_ipqf_coeffs2[t][i];
+                out[s * 16 + i + 8] += hist->buf1[pos_now][7 - i]  * ff_ipqf_coeffs1[t][i + 8] +
+                                       hist->buf2[pos_next][7 - i] * ff_ipqf_coeffs2[t][i + 8];
             }
             pos_now  = mod23_lut[pos_next + 2]; // pos_now  = (pos_now  + 2) % 23;
             pos_next = mod23_lut[pos_now + 2];  // pos_next = (pos_next + 2) % 23;
diff --git a/src/atrac/atrac3plus_pqf/ut/atrac3plusdsp.h b/src/atrac/atrac3plus_pqf/ut/atrac3plusdsp.h
index af45d92..5418f0a 100644
--- a/src/atrac/atrac3plus_pqf/ut/atrac3plusdsp.h
+++ b/src/atrac/atrac3plus_pqf/ut/atrac3plusdsp.h
@@ -36,6 +36,8 @@ extern "C"
 {
 #endif /* __cplusplus */
 
+extern const float ipqf_coeffs1[ATRAC3P_PQF_FIR_LEN][16];
+extern const float ipqf_coeffs2[ATRAC3P_PQF_FIR_LEN][16];
 void ff_atrac3p_ipqf(Atrac3pIPQFChannelCtx *hist, const float *in, float *out);
 
 #ifdef __cplusplus
diff --git a/src/atrac/atrac3plus_pqf/ut/ipqf_ut.cpp b/src/atrac/atrac3plus_pqf/ut/ipqf_ut.cpp
index ababfcd..c7eec07 100644
--- a/src/atrac/atrac3plus_pqf/ut/ipqf_ut.cpp
+++ b/src/atrac/atrac3plus_pqf/ut/ipqf_ut.cpp
@@ -17,6 +17,8 @@
  */
 
 #include "atrac3plusdsp.h"
+#include "../atrac3plus_pqf.h"
+#include "../atrac3plus_pqf_data.h"
 
 #include <gtest/gtest.h>
 
@@ -25,6 +27,8 @@
 #include <string.h>
 #include <math.h>
 
+#include <cstdlib>
+
 #define SAMPLES 8192
 
 static int read_file(FILE* f, float buf[SAMPLES]) {
@@ -36,6 +40,14 @@ static int read_file(FILE* f, float buf[SAMPLES]) {
     return 0;
 }
 
+static void create_chirp(int sz, float* buf) {
+    int i = 0;
+    for(i = 0; i < sz; i++) {
+        float t = i;
+        buf[i] = sinf((t + t * t * 0.5 / 2.0) * 2.0 * M_PI/(float)sz);
+    }
+}
+
 TEST(ipqf, CheckOnRefData) {
     FILE* mr_f = fopen("test_data/ipqftest_pcm_mr.dat", "r");
     if (!mr_f) {
@@ -74,7 +86,210 @@ TEST(ipqf, CheckOnRefData) {
 
     const static float err = 1.0 / (float)(1<<26);
     for (int i = 0; i < SAMPLES; i++) {
-        //fprintf(stdout, "%f %f\n", tmp[i], ref_data[i]);
         EXPECT_NEAR(tmp[i], ref_data[i], err);
     }
 }
+
+TEST(ipqf, CmpEnergy) {
+   double e1 = 0.0;
+   double e2 = 0.0;
+   double e = 0.0;
+   const static double err = 1.0 / (double)(1ull<<32);
+   for (int i = 0; i < ATRAC3P_SUBBANDS; i++) {
+       e1 = 0.0;
+       e2 = 0.0;
+       for (int j = 0; j < ATRAC3P_PQF_FIR_LEN; j++) {
+           e1 += ff_ipqf_coeffs1[j][i] * ff_ipqf_coeffs1[j][i];
+           e2 += ff_ipqf_coeffs2[j][i] * ff_ipqf_coeffs2[j][i];
+       }
+       if (i) {
+           EXPECT_NEAR(e, e1 + e2, err);
+       }
+       e = e1 + e2;
+   }
+}
+
+TEST(pqf, DC_Short) {
+    int i = 0;
+    float x[2048] = {0};
+    float subbands[2048] = {0};
+    for (i = 0; i < 2048; i++)
+        x[i] = 1;
+
+    at3plus_pqf_a_ctx_t actx = at3plus_pqf_create_a_ctx();
+
+    at3plus_pqf_do_analyse(actx, x, subbands);
+
+    float tmp[2048] = {0};
+
+    Atrac3pIPQFChannelCtx sctx;
+    memset(&sctx, 0, sizeof(sctx));
+
+
+    ff_atrac3p_ipqf(&sctx, &subbands[0], &tmp[0]);
+
+    const static float err = 1.0 / (float)(1<<21);
+
+    for (int i = 368; i < 2048; i++) {
+        EXPECT_NEAR(tmp[i], x[i], err);
+    }
+}
+
+TEST(pqf, DC_Long) {
+    int i = 0;
+    float x[4096] = {0};
+    float subbands[4096] = {0};
+    for (i = 0; i < 4096; i++)
+        x[i] = 1;
+
+    at3plus_pqf_a_ctx_t actx = at3plus_pqf_create_a_ctx();
+
+    at3plus_pqf_do_analyse(actx, x, subbands);
+    at3plus_pqf_do_analyse(actx, x + 2048, subbands + 2048);
+
+    float tmp[4096] = {0};
+
+    Atrac3pIPQFChannelCtx sctx;
+    memset(&sctx, 0, sizeof(sctx));
+
+    ff_atrac3p_ipqf(&sctx, &subbands[0], &tmp[0]);
+    ff_atrac3p_ipqf(&sctx, &subbands[2048], &tmp[2048]);
+
+    const static float err = 1.0 / (float)(1<<21);
+
+    for (int i = 368; i < 4096; i++) {
+        EXPECT_NEAR(tmp[i], x[i-368], err);
+    }
+}
+
+TEST(pqf, Seq_Short) {
+    int i = 0;
+    float x[2048] = {0};
+    float subbands[2048] = {0};
+
+    for (i = 0; i < 2048; i++)
+        x[i] = i;
+
+    at3plus_pqf_a_ctx_t actx = at3plus_pqf_create_a_ctx();
+
+    at3plus_pqf_do_analyse(actx, x, subbands);
+
+    float tmp[2048] = {0};
+
+    Atrac3pIPQFChannelCtx sctx;
+    memset(&sctx, 0, sizeof(sctx));
+
+
+    ff_atrac3p_ipqf(&sctx, &subbands[0], &tmp[0]);
+
+    const static float err = 2048.0 / (float)(1<<22);
+
+    for (int i = 368; i < 2048; i++) {
+        EXPECT_NEAR(tmp[i], x[i - 368], err);
+    }
+}
+
+TEST(pqf, Seq_Long) {
+    int i = 0;
+    float x[4096] = {0};
+    float subbands[4096] = {0};
+    for (i = 0; i < 4096; i++)
+        x[i] = i;
+
+    at3plus_pqf_a_ctx_t actx = at3plus_pqf_create_a_ctx();
+
+    at3plus_pqf_do_analyse(actx, x, subbands);
+    at3plus_pqf_do_analyse(actx, x + 2048, subbands + 2048);
+
+    float tmp[4096] = {0};
+
+    Atrac3pIPQFChannelCtx sctx;
+    memset(&sctx, 0, sizeof(sctx));
+
+    ff_atrac3p_ipqf(&sctx, &subbands[0], &tmp[0]);
+    ff_atrac3p_ipqf(&sctx, &subbands[2048], &tmp[2048]);
+
+    const static float err = 4096.0 / (float)(1<<21);
+    for (int i = 368; i < 4096; i++) {
+        EXPECT_NEAR(tmp[i], x[i-368], err);
+    }
+}
+
+TEST(pqf, Chirp_Short) {
+    int i = 0;
+    float x[2048] = {0};
+    float subbands[2048] = {0};
+
+    create_chirp(2048, x);
+
+    at3plus_pqf_a_ctx_t actx = at3plus_pqf_create_a_ctx();
+
+    at3plus_pqf_do_analyse(actx, x, subbands);
+
+    float tmp[2048] = {0};
+
+    Atrac3pIPQFChannelCtx sctx;
+    memset(&sctx, 0, sizeof(sctx));
+
+    ff_atrac3p_ipqf(&sctx, &subbands[0], &tmp[0]);
+
+    const static float err = 1.0 / (float)(1<<22);
+
+    for (int i = 368; i < 2048; i++) {
+        EXPECT_NEAR(tmp[i], x[i - 368], err);
+    }
+}
+
+TEST(pqf, Chirp_Long) {
+    int i = 0;
+    float x[4096] = {0};
+    float subbands[4096] = {0};
+
+    create_chirp(4096, x);
+
+    at3plus_pqf_a_ctx_t actx = at3plus_pqf_create_a_ctx();
+
+    at3plus_pqf_do_analyse(actx, x, subbands);
+    at3plus_pqf_do_analyse(actx, x + 2048, subbands + 2048);
+
+    float tmp[4096] = {0};
+
+    Atrac3pIPQFChannelCtx sctx;
+    memset(&sctx, 0, sizeof(sctx));
+
+    ff_atrac3p_ipqf(&sctx, &subbands[0], &tmp[0]);
+    ff_atrac3p_ipqf(&sctx, &subbands[2048], &tmp[2048]);
+
+    const static float err = 4096.0 / (float)(1<<21);
+    for (int i = 368; i < 4096; i++) {
+        EXPECT_NEAR(tmp[i], x[i-368], err);
+    }
+}
+
+TEST(pqf, Noise_Long) {
+    int i = 0;
+    float x[4096] = {0};
+    float subbands[4096] = {0};
+    for (i = 0; i < 4096; i++)
+        x[i] = (float)rand() / (float)RAND_MAX - 0.5;
+
+    at3plus_pqf_a_ctx_t actx = at3plus_pqf_create_a_ctx();
+
+    at3plus_pqf_do_analyse(actx, x, subbands);
+    at3plus_pqf_do_analyse(actx, x + 2048, subbands + 2048);
+
+    float tmp[4096] = {0};
+
+    Atrac3pIPQFChannelCtx sctx;
+    memset(&sctx, 0, sizeof(sctx));
+
+    ff_atrac3p_ipqf(&sctx, &subbands[0], &tmp[0]);
+    ff_atrac3p_ipqf(&sctx, &subbands[2048], &tmp[2048]);
+
+    const static float err = 1.0 / (float)(1<<21);
+    for (int i = 368; i < 4096; i++) {
+        EXPECT_NEAR(tmp[i], x[i-368], err);
+    }
+}
+
+