[AT3P] Use fast DCT-IV calculation for PQF

Author: dcherednik

src/atrac/atrac3plus_pqf/atrac3plus_pqf.c

@@ -30,6 +30,8 @@
 #include "atrac3plus_pqf.h"
 #include "atrac3plus_pqf_data.h"
 
+#include "../../mdct/dct.h"
+
 /*
  * Number of subbands to split input signal
  */
@@ -49,18 +51,9 @@
 
 static float fir[PROTO_SZ];
 
-static void dct4(float* out, const float* x, int N, float scale) {
-    for (int k = 0; k < N; k++) {
-        double sum = 0;
-        for (int n = 0; n < N; n++) {
-            sum += ((double)x[n] * cosl((M_PI/(double)N) * ((double)n + 0.5) * ((double)k + 0.5)));
-        }
-        out[N - 1 - k] = sum * scale;
-    }
-}
-
 struct at3plus_pqf_a_ctx {
     float buf[FRAME_SZ + OVERLAP_SZ];
+    atde_dct_ctx_t dct_ctx;
 };
 
 static void init(void)
@@ -95,7 +88,7 @@ static void vectoring(const float* const x, float* y)
     }
 }
 
-static void matrixing(const float* y, float* samples )
+static void matrixing(atde_dct_ctx_t ctx, const float* y, float* samples )
 {
     float yy[SUBBANDS_NUM];
     float res[SUBBANDS_NUM];
@@ -105,7 +98,7 @@ static void matrixing(const float* y, float* samples )
         yy[i + 8] = y[i + 16] + y[31 - i];
     }
 
-    dct4(res, yy, SUBBANDS_NUM, 128.0 * 512.0);
+    atde_do_dct4_16(ctx, yy, res);
 
     for (int i = 0; i < SUBBANDS_NUM; i++) {
         samples[i * SUBBAND_SIZE] = res[SUBBANDS_NUM - 1 - i];
@@ -120,13 +113,17 @@ at3plus_pqf_a_ctx_t at3plus_pqf_create_a_ctx()
         ctx->buf[i] = 0.0;
     }
 
+    ctx->dct_ctx = atde_create_dct4_16(128 * 512.0);
+
     init();
 
     return ctx;
 }
 
 void at3plus_pqf_free_a_ctx(at3plus_pqf_a_ctx_t ctx)
 {
+    atde_free_dct_ctx(ctx->dct_ctx);
+
     free(ctx);
 }
 
@@ -142,7 +139,7 @@ void at3plus_pqf_do_analyse(at3plus_pqf_a_ctx_t ctx, const float* in, float* out
 
     for (int i = 0; i < SUBBAND_SIZE; i++) {
         vectoring(x, y);
-        matrixing (y, &out[i]);
+        matrixing (ctx->dct_ctx, y, &out[i]);
         x += SUBBANDS_NUM;
     }
 

src/mdct/mdct.cpp

@@ -17,6 +17,7 @@
  */
 
 #include "mdct.h"
+#include "dct.h"
 #include <iostream>
 
 namespace NMDCT {
@@ -51,3 +52,31 @@ TMDCTBase::~TMDCTBase()
 }
 
 } // namespace NMDCT
+
+struct atde_dct_ctx {
+    atde_dct_ctx(float scale)
+        : mdct(scale)
+    {}
+    NMDCT::TMIDCT<32, float> mdct;
+};
+
+atde_dct_ctx_t atde_create_dct4_16(float scale)
+{
+    return new atde_dct_ctx(32.0 * scale);
+}
+
+void atde_free_dct_ctx(atde_dct_ctx_t ctx)
+{
+    delete ctx;
+}
+
+void atde_do_dct4_16(atde_dct_ctx_t ctx, const float* in, float* out)
+{
+    //TODO: rewrire more optimal
+    const auto& x = ctx->mdct(in);
+
+    for (int i = 0; i < 16; i++) {
+        out[i] = x[i + 8] * -1.0;
+    }
+}
+

src/mdct/mdct.h

@@ -39,16 +39,16 @@ class TMDCTBase {
 };
 
 
-template<size_t TN>
+template<size_t TN, typename TIO = TFloat>
 class TMDCT : public TMDCTBase {
-    std::vector<TFloat> Buf;
+    std::vector<TIO> Buf;
 public:
     TMDCT(float scale = 1.0)
         : TMDCTBase(TN, scale)
         , Buf(TN/2)
     {
     }
-    const std::vector<TFloat>& operator()(TFloat* in) {
+    const std::vector<TIO>& operator()(const TIO* in) {
 
         const size_t n2 = N >> 1;
         const size_t n4 = N >> 2;
@@ -104,15 +104,15 @@ class TMDCT : public TMDCTBase {
     }
 };
 
-template<size_t TN>
+template<size_t TN, typename TIO = TFloat>
 class TMIDCT : public TMDCTBase {
-    std::vector<TFloat> Buf;
+    std::vector<TIO> Buf;
 public:
     TMIDCT(float scale = TN)
         : TMDCTBase(TN, scale/2)
         , Buf(TN)
     {}
-    const std::vector<TFloat>& operator()(TFloat* in) {
+    const std::vector<TIO>& operator()(const TIO* in) {
 
         const size_t n2 = N >> 1;
         const size_t n4 = N >> 2;

test/CMakeLists.txt

@@ -37,12 +37,14 @@ set(at3plus_pqf_ut
     ../src/atrac/atrac3plus_pqf/ut/ipqf_ut.cpp
     ../src/atrac/atrac3plus_pqf/ut/atrac3plusdsp.c
     ../src/atrac/atrac3plus_pqf/atrac3plus_pqf.c
+    ../src/mdct/mdct.cpp
 )
 
 add_executable(at3plus_pqf_ut ${at3plus_pqf_ut})
 
 target_link_libraries(at3plus_pqf_ut
     m
+    fft_impl
     GTest::gtest_main
 )