fast MIDCT implementation added

6 files changed, 204 insertions, 4 deletions

diff --git a/src/atracdenc.cpp b/src/atracdenc.cpp
index c2ecad7..f9afaf5 100644
--- a/src/atracdenc.cpp
+++ b/src/atracdenc.cpp
@@ -46,7 +46,7 @@ vector<double> midct(double* x, int N) {
         for (int k = 0; k < N; k++) {
             sum += (x[k] * cos((M_PI/N) * ((double)n + 0.5 + N/2) * ((double)k + 0.5)));
         }
-        sum = sum / N;
+
         res.push_back(sum);
     }
     return res;
@@ -101,7 +101,6 @@ void TAtrac1Processor::IMdct(double Specs[512], const TBlockSize& mode, double*
         uint32_t start = 0;
 
         double* dstBuf = (band == 0) ? low : (band == 1) ? mid : hi;
-        const double gainMul = (band == 2) ? 256.0 : 128.0;
 
 
         vector<double> invBuf;
@@ -118,10 +117,10 @@ void TAtrac1Processor::IMdct(double Specs[512], const TBlockSize& mode, double*
                     Specs[pos + blockSize - 1 -j] = tmp;
                 }
             }
-            vector<double> inv = midct(&Specs[pos], blockSize);
+            vector<double> inv = (blockSize == 32) ? midct(&Specs[pos], blockSize) : (blockSize == 128) ? Midct256(&Specs[pos]) : Midct512(&Specs[pos]);
             for (int i = 0; i < (inv.size()/2); i++) {
 
-                invBuf[start+i] = ((blockSize == 32) ? 32.0 : gainMul) * inv[i + inv.size()/4];
+                invBuf[start+i] = inv[i + inv.size()/4];
             }
 
             vector_fmul_window(dstBuf + start, prevBuf, &invBuf[start], &TAtrac1Data::SineWindow[0], 16);
diff --git a/src/atracdenc.h b/src/atracdenc.h
index a0ab786..c1bce48 100644
--- a/src/atracdenc.h
+++ b/src/atracdenc.h
@@ -25,6 +25,9 @@ class TAtrac1Processor : public TAtrac1Data {
     Atrac1SplitFilterBank<double> SplitFilterBank[2];
     NMDCT::TMDCT<512> Mdct512;
     NMDCT::TMDCT<256> Mdct256;
+    NMDCT::TMIDCT<512> Midct512;
+    NMDCT::TMIDCT<256> Midct256;
+ 
     void IMdct(double specs[512], const TBlockSize& mode, double* low, double* mid, double* hi);
     void Mdct(double specs[512], double* low, double* mid, double* hi);
     NAtrac1::TScaler Scaler;
diff --git a/src/mdct/mdct.h b/src/mdct/mdct.h
index 0637651..cbb842f 100644
--- a/src/mdct/mdct.h
+++ b/src/mdct/mdct.h
@@ -29,6 +29,19 @@ public:
         mdct(&Ctx, &Buf[0], in);
         return Buf;
     }
+};
 
+template<int N>
+class TMIDCT : public TMDCTBase {
+    std::vector<double> Buf;
+public:
+    TMIDCT(float scale = 1.0)
+        : Buf(N)
+        , TMDCTBase(N, scale)
+    {}
+    const std::vector<double>& operator()(double* in) {
+        midct(&Ctx, &Buf[0], in);
+        return Buf;
+    }
 };
 } //namespace NMDCT
diff --git a/src/mdct/mdct_ut.cpp b/src/mdct/mdct_ut.cpp
index a63544e..e81bea1 100644
--- a/src/mdct/mdct_ut.cpp
+++ b/src/mdct/mdct_ut.cpp
@@ -2,6 +2,7 @@
 #include <gtest/gtest.h>
 
 #include <vector>
+#include <cstdlib>
 
 using std::vector;
 using namespace NMDCT;
@@ -18,6 +19,19 @@ static vector<double> mdct(double* x, int N) {
     return res;
 }  
 
+static vector<double> midct(double* x, int N) {
+    vector<double> res;
+    for (int n = 0; n < 2 * N; n++) {
+        double sum = 0;
+        for (int k = 0; k < N; k++)
+            sum += (x[k] * cos((M_PI/N) * ((double)n + 0.5 + N/2) * ((double)k + 0.5)));
+
+        res.push_back(sum);
+    }
+    return res;
+}
+
+
 TEST(TBitStream, MDCT64) {
     const int N = 64;
     TMDCT<N> transform(N);
@@ -62,3 +76,80 @@ TEST(TBitStream, MDCT256) {
         EXPECT_NEAR(res1[i], res2[i], 0.00000001);
     }
 }
+
+TEST(TBitStream, MDCT256_RAND) {
+    const int N = 256;
+    TMDCT<N> transform(N);
+    vector<double> src(N);
+    for (int i = 0; i < N; i++) {
+        src[i] = rand();
+    }
+    const vector<double> res1 = mdct(&src[0], N/2);
+    const vector<double> res2 = transform(&src[0]);
+    EXPECT_EQ(res1.size(), res2.size());
+    for (int i = 0; i < res1.size(); i++) {
+        EXPECT_NEAR(res1[i], res2[i], 0.01);
+    }
+}
+
+
+TEST(TBitStream, MIDCT64) {
+    const int N = 64;
+    TMIDCT<N> transform(1);
+    vector<double> src(N);
+    for (int i = 0; i < N/2; i++) {
+        src[i] = i;
+    }
+    const vector<double> res1 = midct(&src[0], N/2);
+    const vector<double> res2 = transform(&src[0]);
+    EXPECT_EQ(res1.size(), res2.size());
+    for (int i = 0; i < N; i++) {
+        EXPECT_NEAR(res1[i], res2[i], 0.0000000001);
+    }
+}
+
+TEST(TBitStream, MIDCT128) {
+    const int N = 128;
+    TMIDCT<N> transform(1);
+    vector<double> src(N);
+    for (int i = 0; i < N/2; i++) {
+        src[i] = i;
+    }
+    const vector<double> res1 = midct(&src[0], N/2);
+    const vector<double> res2 = transform(&src[0]);
+    EXPECT_EQ(res1.size(), res2.size());
+    for (int i = 0; i < N; i++) {
+        EXPECT_NEAR(res1[i], res2[i], 0.0000000001);
+    }
+}
+
+TEST(TBitStream, MIDCT256) {
+    const int N = 256;
+    TMIDCT<N> transform(1);
+    vector<double> src(N);
+    for (int i = 0; i < N/2; i++) {
+        src[i] = i;
+    }
+    const vector<double> res1 = midct(&src[0], N/2);
+    const vector<double> res2 = transform(&src[0]);
+    EXPECT_EQ(res1.size(), res2.size());
+    for (int i = 0; i < N; i++) {
+        EXPECT_NEAR(res1[i], res2[i], 0.000000001);
+    }
+}
+
+TEST(TBitStream, MIDCT256_RAND) {
+    const int N = 256;
+    TMIDCT<N> transform(1);
+    vector<double> src(N);
+    for (int i = 0; i < N/2; i++) {
+        src[i] = rand();
+    }
+    const vector<double> res1 = midct(&src[0], N/2);
+    const vector<double> res2 = transform(&src[0]);
+    EXPECT_EQ(res1.size(), res2.size());
+    for (int i = 0; i < N; i++) {
+        EXPECT_NEAR(res1[i], res2[i], 0.01);
+    }
+}
+
diff --git a/src/mdct/vorbis_impl/mdct.c b/src/mdct/vorbis_impl/mdct.c
index b297089..baa7b43 100644
--- a/src/mdct/vorbis_impl/mdct.c
+++ b/src/mdct/vorbis_impl/mdct.c
@@ -448,3 +448,96 @@ mdct(MDCTContext *mdct, FLOAT *out, FLOAT *in)
     }
 }
 
+void
+midct(MDCTContext *mdct, FLOAT *out, FLOAT *in)
+{
+    int n = mdct->n;
+    int n2 = n>>1;
+    int n4 = n>>2;
+
+    FLOAT *iX = in+n2-7;
+    FLOAT *oX = out+n2+n4;
+    FLOAT *T  = mdct->trig+n4;
+    do {
+        oX -= 4;
+        oX[0] = (-iX[2] * T[3] - iX[0]  * T[2]);
+        oX[1] = (iX[0] * T[3] - iX[2]  * T[2]);
+        oX[2] = (-iX[6] * T[1] - iX[4]  * T[0]);
+        oX[3] = (iX[4] * T[1] - iX[6]  * T[0]);
+        iX -= 8;
+        T += 4;
+    } while (iX >= in);
+
+    iX = in+n2-8;
+    oX = out+n2+n4;
+    T = mdct->trig+n4;
+
+    do {
+        T -= 4;
+        oX[0] = (iX[4] * T[3] + iX[6] * T[2]);
+        oX[1] = (iX[4] * T[2] - iX[6] * T[3]);
+        oX[2] = (iX[0] * T[1] + iX[2] * T[0]);
+        oX[3] = (iX[0] * T[0] - iX[2] * T[1]);
+        iX -= 8;
+        oX += 4;
+    } while (iX >= in);
+
+    mdct_butterflies(mdct, out+n2, n2);
+    mdct_bitreverse(mdct, out);
+
+    /* roatate + window */
+
+    {
+        FLOAT *oX1= out + n2 + n4;
+        FLOAT *oX2= out + n2 + n4;
+        FLOAT *iX = out;
+        T = mdct->trig + n2;
+
+        do {
+            oX1-=4;
+            oX1[3] = (iX[0] * T[1] - iX[1] * T[0]);
+            oX2[0] = -(iX[0] * T[0] + iX[1] * T[1]);
+
+            oX1[2] = (iX[2] * T[3] - iX[3] * T[2]);
+            oX2[1] = -(iX[2] * T[2] + iX[3] * T[3]);
+
+            oX1[1] = (iX[4] * T[5] - iX[5] * T[4]);
+            oX2[2] = -(iX[4] * T[4] + iX[5] * T[5]);
+
+            oX1[0] = (iX[6] * T[7] - iX[7] * T[6]);
+            oX2[3] = -(iX[6] * T[6] + iX[7] * T[7]);
+
+            oX2 += 4;
+            iX += 8;
+            T += 8;
+        } while (iX < oX1);
+
+        iX = out + n2 + n4;
+        oX1 = out + n4;
+        oX2 = oX1;
+
+        do{
+            oX1-=4;
+            iX-=4;
+            oX2[0] = -(oX1[3] = iX[3]);
+            oX2[1] = -(oX1[2] = iX[2]);
+            oX2[2] = -(oX1[1] = iX[1]);
+            oX2[3] = -(oX1[0] = iX[0]);
+            oX2+=4;
+        } while (oX2 < iX);
+
+        iX=out+n2+n4;
+        oX1=out+n2+n4;
+        oX2=out+n2;
+
+        do{
+            oX1-=4;
+            oX1[0]= iX[3];
+            oX1[1]= iX[2];
+            oX1[2]= iX[1];
+            oX1[3]= iX[0];
+            iX+=4;
+        } while (oX1 > oX2);
+    }
+
+}
diff --git a/src/mdct/vorbis_impl/mdct.h b/src/mdct/vorbis_impl/mdct.h
index c80a351..de415af 100644
--- a/src/mdct/vorbis_impl/mdct.h
+++ b/src/mdct/vorbis_impl/mdct.h
@@ -87,6 +87,7 @@ extern "C" {
 void mdct_ctx_init(MDCTContext *mdct, int n, FLOAT scale);
 void mdct_ctx_close(MDCTContext *mdct);
 void mdct(MDCTContext *mdct, FLOAT *out, FLOAT *in);
+void midct(MDCTContext *mdct, FLOAT *out, FLOAT *in);
 #ifdef __cplusplus
 }