[AT3P] Part of naive GHA processor

* extract sinusoids parameters one by one. It is simple and fast but
  probably will not work well on the real musical signals

* envelope is not passed to the output. Envelope handling must be improved

6 files changed, 908 insertions, 22 deletions

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 4242f99..2ac6049 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -103,6 +103,7 @@ set(SOURCE_ATRACDENC_IMPL
     atrac/atrac3_bitstream.cpp
     atrac/at3p/at3p.cpp
     atrac/at3p/at3p_bitstream.cpp
+    atrac/at3p/at3p_gha.cpp
     atrac/at3p/at3p_tables.cpp
     lib/mdct/mdct.cpp
 )
@@ -115,7 +116,7 @@ endif()
 add_library(oma STATIC ${SOURCE_OMA_LIB})
 add_library(bitstream STATIC ${SOURCE_BITSTREAM_LIB})
 add_library(atracdenc_impl STATIC ${SOURCE_ATRACDENC_IMPL})
-target_link_libraries(atracdenc_impl fft_impl pcm_io oma bitstream ${SNDFILE_LIBRARIES})
+target_link_libraries(atracdenc_impl fft_impl pcm_io oma bitstream ${SNDFILE_LIBRARIES} gha)
 set(SOURCE_EXE
     main.cpp
     help.cpp
diff --git a/src/atrac/at3p/at3p_bitstream_ut.cpp b/src/atrac/at3p/at3p_bitstream_ut.cpp
index 5fe6902..3129673 100644
--- a/src/atrac/at3p/at3p_bitstream_ut.cpp
+++ b/src/atrac/at3p/at3p_bitstream_ut.cpp
@@ -7,7 +7,7 @@ using namespace NAtracDEnc;
 TEST(AT3PBitstream, ToneFreqBitPack__1) {
     std::vector<TAt3PGhaData::TWaveParam> params;
 
-    params.push_back(TAt3PGhaData::TWaveParam{1, 0, 0, 0, {}});
+    params.push_back(TAt3PGhaData::TWaveParam{1, 0, 0, 0});
 
     auto r = CreateFreqBitPack(params.data(), params.size());
     EXPECT_EQ(r.UsedBits, 10);
@@ -20,9 +20,9 @@ TEST(AT3PBitstream, ToneFreqBitPack__1) {
 TEST(AT3PBitstream, ToneFreqBitPack__512_1020_1023) {
     std::vector<TAt3PGhaData::TWaveParam> params;
 
-    params.push_back(TAt3PGhaData::TWaveParam{512, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{1020, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{1023, 0, 0, 0, {}});
+    params.push_back(TAt3PGhaData::TWaveParam{512, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{1020, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{1023, 0, 0, 0});
 
     auto r = CreateFreqBitPack(params.data(), params.size());
     EXPECT_EQ(r.UsedBits, 21);
@@ -39,9 +39,9 @@ TEST(AT3PBitstream, ToneFreqBitPack__512_1020_1023) {
 TEST(AT3PBitstream, ToneFreqBitPack__1_2_3) {
     std::vector<TAt3PGhaData::TWaveParam> params;
 
-    params.push_back(TAt3PGhaData::TWaveParam{1, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{2, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{3, 0, 0, 0, {}});
+    params.push_back(TAt3PGhaData::TWaveParam{1, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{2, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{3, 0, 0, 0});
 
     auto r = CreateFreqBitPack(params.data(), params.size());
     EXPECT_EQ(r.UsedBits, 14);
@@ -58,12 +58,12 @@ TEST(AT3PBitstream, ToneFreqBitPack__1_2_3) {
 TEST(AT3PBitstream, ToneFreqBitPack__1_2_3_1020_1021_1022) {
     std::vector<TAt3PGhaData::TWaveParam> params;
 
-    params.push_back(TAt3PGhaData::TWaveParam{1, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{2, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{3, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{1020, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{1021, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{1022, 0, 0, 0, {}});
+    params.push_back(TAt3PGhaData::TWaveParam{1, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{2, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{3, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{1020, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{1021, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{1022, 0, 0, 0});
 
     auto r = CreateFreqBitPack(params.data(), params.size());
     EXPECT_EQ(r.UsedBits, 44);
@@ -86,11 +86,11 @@ TEST(AT3PBitstream, ToneFreqBitPack__1_2_3_1020_1021_1022) {
 TEST(AT3PBitstream, ToneFreqBitPack__1_2_1020_1021_1022) {
     std::vector<TAt3PGhaData::TWaveParam> params;
 
-    params.push_back(TAt3PGhaData::TWaveParam{1, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{2, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{1020, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{1021, 0, 0, 0, {}});
-    params.push_back(TAt3PGhaData::TWaveParam{1022, 0, 0, 0, {}});
+    params.push_back(TAt3PGhaData::TWaveParam{1, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{2, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{1020, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{1021, 0, 0, 0});
+    params.push_back(TAt3PGhaData::TWaveParam{1022, 0, 0, 0});
 
     auto r = CreateFreqBitPack(params.data(), params.size());
     EXPECT_EQ(r.UsedBits, 34);
diff --git a/src/atrac/at3p/at3p_gha.cpp b/src/atrac/at3p/at3p_gha.cpp
new file mode 100644
index 0000000..4dec267
--- /dev/null
+++ b/src/atrac/at3p/at3p_gha.cpp
@@ -0,0 +1,483 @@
+/*
+ * This file is part of AtracDEnc.
+ *
+ * AtracDEnc 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.
+ *
+ * AtracDEnc 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 AtracDEnc; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "at3p_gha.h"
+
+#include <assert.h>
+#include <atrac/atrac_psy_common.h>
+#include <libgha/include/libgha.h>
+
+#include <memory>
+
+#include <iostream>
+#include <map>
+#include <vector>
+
+
+using std::map;
+using std::vector;
+using std::isnan;
+using std::pair;
+using std::max;
+using std::min;
+
+namespace NAtracDEnc {
+
+namespace {
+
+uint32_t GhaFreqToIndex(float f, uint32_t sb)
+{
+    return static_cast<uint32_t>(lrintf(1024.0f * (f / M_PI)) & 1023) | (sb << 10);
+}
+
+uint32_t GhaPhaseToIndex(float p)
+{
+    return static_cast<uint32_t>(lrintf(32.0 * (p / (2 * M_PI))) & 31);
+}
+
+class TGhaProcessor : public IGhaProcessor {
+    // Number of subbands to process;
+    // No need to process all subbands.
+    static const size_t SUBBANDS = 13;
+    static const size_t CHANNEL_BUF_SZ = SUBBANDS * 128;;
+
+    using TGhaInfoMap = map<uint32_t, struct gha_info>;
+    using TWavesChannel = TAt3PGhaData::TWavesChannel;
+    using TAmpSfTab = std::array<float, 64>;
+
+    struct TChannelData {
+        const float* SrcBuf;
+        float Buf[CHANNEL_BUF_SZ];
+        pair<uint32_t, uint32_t> Envelopes[SUBBANDS] = {{0,0}};
+        uint8_t SubbandDone[SUBBANDS] = {0};
+        TGhaInfoMap GhaInfos;
+
+        void MarkSubbandDone(size_t sb) {
+            SubbandDone[sb] = 16;
+        }
+
+        bool IsSubbandDone(size_t sb) const {
+            return SubbandDone[sb] == 16;
+        }
+    };
+
+public:
+    TGhaProcessor(float* b1, float* b2)
+        : B1(b1)
+        , B2(b2)
+        , LibGhaCtx(gha_create_ctx(128))
+        , AmpSfTab(CreateAmpSfTab())
+    {
+        FillSubbandAth(&SubbandAth[0]);
+    }
+
+    const TAt3PGhaData* DoAnalize() override;
+private:
+    static void FillSubbandAth(float* out) {
+        const auto ath = CalcATH(16 * 1024, 44100);
+        #pragma GCC nounroll
+        for (size_t sb = 0; sb < SUBBANDS; sb++) {
+            float m = 999.;
+            for (size_t f = sb * 1024, i = 0; i < 1024; f++, i++) {
+                m = fmin(m, ath[f]);
+            }
+            m += 26; //Some gap to not encode too much
+            out[sb] = pow(10, 0.1 * (m + 90)); //adjust to 0db level = 32768, convert to power
+        }
+    }
+
+    static TAmpSfTab CreateAmpSfTab() {
+        TAmpSfTab AmpSfTab;
+        for (int i = 0; i < (int)AmpSfTab.size(); i++) {
+            AmpSfTab[i] = exp2f((i - 3) / 4.0f);
+        }
+        return AmpSfTab;
+    }
+
+    uint32_t FillFolowerRes(const TGhaInfoMap& fGhaInfos, TGhaInfoMap::const_iterator& it, uint32_t leaderSb);
+
+    uint32_t AmplitudeToSf(float amp);
+
+    bool DoRound(TChannelData& data, size_t& totalTones) const;
+    bool PsyPreCheck(size_t sb, const struct gha_info& gha) const;
+    pair<uint32_t, uint32_t> FindPos(const float* src, const float* analized, float magn, uint32_t freqIdx) const;
+    void FillResultBuf(const vector<TChannelData>& data);
+
+    bool IsStereo() const {
+        return (bool)B2;
+    }
+
+    float* const B1;
+    float* const B2;
+    gha_ctx_t LibGhaCtx;
+    const TAmpSfTab AmpSfTab;
+    float SubbandAth[SUBBANDS];
+    TAt3PGhaData ResultBuf;
+};
+
+const TAt3PGhaData* TGhaProcessor::DoAnalize()
+{
+    vector<TChannelData> data((size_t)IsStereo() + 1);
+    bool progress[2] = {false};
+
+    for (size_t ch = 0; ch < data.size(); ch++) {
+        const float* b = (ch == 0) ? B1 : B2;
+        data[ch].SrcBuf = b;
+        memcpy(&data[ch].Buf[0], b, sizeof(float) * CHANNEL_BUF_SZ);
+    }
+
+    size_t totalTones = 0;
+    do {
+        for (size_t ch = 0; ch < data.size(); ch++) {
+            progress[ch] = DoRound(data[ch], totalTones);
+        }
+    } while ((progress[0] || progress[1]) && totalTones < 48);
+
+    if (totalTones == 0) {
+        return nullptr;
+    }
+
+    FillResultBuf(data);
+
+    return &ResultBuf;
+}
+
+bool TGhaProcessor::DoRound(TChannelData& data, size_t& totalTones) const
+{
+    bool progress = false;
+    for (size_t sb = 0; sb < SUBBANDS; sb++) {
+        if (data.IsSubbandDone(sb)) {
+            continue;
+        }
+
+        if (totalTones > 48) {
+            return false;
+        }
+
+        float* b = &data.Buf[sb * 128];
+        struct gha_info res;
+
+        gha_analyze_one(b, &res, LibGhaCtx);
+
+#if 0
+        {
+            const float* srcB = data.SrcBuf + (sb * 128);
+            auto it = data.GhaInfos.lower_bound(sb << 10);
+            vector<gha_info> tmp;
+            for(; it != data.GhaInfos.end() && it->first < (sb + 1) << 10; it++) {
+                tmp.push_back(it->second);
+            }
+            if (tmp.size() > 1) {
+                for (const auto& x : tmp) {
+                    std::cerr << "to adjust: " << x.frequency << " " << x.magnitude << std::endl;
+                }
+                int ar = gha_adjust_info(srcB, tmp.data(), tmp.size(), LibGhaCtx);
+                for (const auto& x : tmp) {
+                    std::cerr << "after adjust: " << x.frequency << " " << x.magnitude << " ar: " << ar << std::endl;
+                }
+
+            }
+        }
+#endif
+        auto freqIndex = GhaFreqToIndex(res.frequency, sb);
+        //std::cerr << "sb: " << sb << " findex: " << freqIndex << " magn " << res.magnitude <<  std::endl;
+        if (PsyPreCheck(sb, res) == false) {
+            data.MarkSubbandDone(sb);
+            //std::cerr << "PreCheck failed" << std::endl;
+        } else {
+            const float* tone = gha_get_analyzed(LibGhaCtx);
+            auto pos = FindPos(b, tone, res.magnitude, freqIndex);
+            if (pos.second == 0) {
+                data.MarkSubbandDone(sb);
+            } else {
+                auto& envelope = data.Envelopes[sb];
+                // Envelope for whole subband not initialized
+                if (envelope.second == 0) {
+                    assert(data.SubbandDone[sb] == 0);
+                    envelope = pos;
+                    bool ins = data.GhaInfos.insert({freqIndex, res}).second;
+                    assert(ins);
+                } else {
+                    //TODO:
+                    // Nth sine for subband was out of existed envelope. This case requires more
+                    // complicated analisys - just skip it for a while
+                    if ((pos.first < envelope.first && pos.second <= envelope.first) ||
+                        (pos.first >= envelope.second)) {
+                        data.MarkSubbandDone(sb);
+                        continue;
+                    }
+
+                    const auto it = data.GhaInfos.lower_bound(freqIndex);
+
+                    const size_t minFreqDistanse = 40; // Now we unable to handle tones with close frequency
+                    if (it != data.GhaInfos.end()) {
+                        if (it->first == freqIndex) {
+                            //std::cerr << "already added: " << freqIndex << std::endl;
+                            data.MarkSubbandDone(sb);
+                            continue;
+                        }
+                        //std::cerr << "right: " << it->first << " " << freqIndex << std::endl;
+                        if (it->first - freqIndex < minFreqDistanse) {
+                            data.MarkSubbandDone(sb);
+                            continue;
+                        }
+                    }
+                    if (it != data.GhaInfos.begin()) {
+                        auto prev = it;
+                        prev--;
+                        //std::cerr << "left: " << prev->first << " " << freqIndex << std::endl;
+                        if (freqIndex - prev->first < minFreqDistanse) {
+                            data.MarkSubbandDone(sb);
+                            continue;
+                        }
+                    }
+                    data.GhaInfos.insert(it, {freqIndex, res});
+                    envelope.first = max(envelope.first, pos.first);
+                    envelope.second = min(envelope.second, pos.second);
+                }
+                // Here we have updated envelope, we are ready to add extracted tone
+
+                for (size_t j = envelope.first; j < envelope.second; j++) {
+                    b[j] -= tone[j] * res.magnitude;
+                }
+
+                data.SubbandDone[sb]++;
+                totalTones++;
+                progress = true;
+
+                //std::cerr << "envelop " << envelope.first << " " << envelope.second << std::endl;
+            }
+        }
+
+    }
+    return progress;
+}
+
+bool TGhaProcessor::PsyPreCheck(size_t sb, const struct gha_info& gha) const
+{
+    if (isnan(gha.magnitude)) {
+        return false;
+    }
+
+    //std::cerr << "sb: " << sb << " " << gha.magnitude << " ath: " << SubbandAth[sb] << std::endl;
+    // TODO: improve it
+    // Just to start. Actualy we need to consider spectral leakage during MDCT
+    if ((gha.magnitude * gha.magnitude) > SubbandAth[sb]) {
+        return true;
+    } else {
+        return false;
+    }
+    return true;
+}
+
+pair<uint32_t, uint32_t> TGhaProcessor::FindPos(const float* src, const float* tone, float magn, uint32_t freqIdx) const
+{
+    freqIdx &= 1023; // clear subband part
+    size_t windowSz = 64;
+    // TODO: place to make experiments
+    // we should compare rms on some window. Optimal window sile depends on frequency so right now
+    // here just some heuristic to match window size to freq index
+
+    // Ideas:
+    // - generate FIR filter for each frequncy index and use it fo filter source signal before compare level
+    // - may be just perform search all sinusoid components, run N dimension GHA optimization, and then find positions
+    if (freqIdx > 256) {
+        windowSz = 8;
+    } else if (freqIdx > 128) {
+        windowSz = 16;
+    } else if (freqIdx > 64) {
+        windowSz = 32;
+    } else {
+        windowSz = 64;
+    }
+
+    uint32_t start = 0;
+    uint32_t count = 0;
+    uint32_t len = 0;
+    bool found = false;
+
+    for (size_t i = 0; i < 128; i += windowSz) {
+        float rmsIn = 0.0;
+        float rmsOut = 0.0;
+        for (size_t j = 0; j < windowSz; j++) {
+            rmsIn += src[i + j] * src[i + j];
+            float r = src[i + j] - tone[i + j] * magn;
+            rmsOut += r * r;
+        }
+        rmsIn /= windowSz;
+        rmsOut /= windowSz;
+        rmsIn = sqrt(rmsIn);
+        rmsOut = sqrt(rmsOut);
+
+        //std::cerr << i << " rms : " << rmsIn << " " << rmsOut  << "\t\t\t" << ((rmsOut < rmsIn) ? "+" : "-")  << std::endl;
+        if (rmsIn / rmsOut < 1) {
+            count = 0;
+            found = false;
+        } else {
+            count++;
+            if (count > len) {
+                len = count;
+                if (!found) {
+                    start = i;
+                    found = true;
+                }
+            }
+        }
+
+    }
+
+    //std::cerr << "start pos: " << start << " len: " << len * windowSz << " end: " << start + len * windowSz << std::endl;
+
+    auto end = start + len * windowSz;
+
+    //if (end - start < 64) {
+    //    return {0, 0};
+    //}
+
+    return {start, end};
+}
+
+void TGhaProcessor::FillResultBuf(const vector<TChannelData>& data)
+{
+    uint32_t usedContiguousSb[2] = {0, 0};
+    uint32_t numTones[2] = {0, 0};
+    //TODO: Calc used SB just during DoRound routine
+    for (size_t ch = 0; ch < data.size(); ch++) {
+        int cur = -1;
+        for (const auto& info : data[ch].GhaInfos) {
+            int sb = info.first >> 10;
+            if (sb == cur + 1) {
+                usedContiguousSb[ch]++;
+                numTones[ch]++;
+                cur = sb;
+            } else if (sb == cur) {
+                numTones[ch]++;
+                continue;
+            } else {
+                break;
+            }
+        }
+        //std::cerr << "ch: " << ch << " usedSb: " << usedContiguousSb[ch] << " numTones: " << numTones[ch] << std::endl;
+    }
+
+    bool leader = usedContiguousSb[1] > usedContiguousSb[0];
+
+    ResultBuf.NumToneBands = usedContiguousSb[leader];
+
+    TGhaInfoMap::const_iterator folowerIt;
+    if (data.size() == 2) {
+        TWavesChannel& fWaves = ResultBuf.Waves[1];
+        fWaves.WaveParams.clear();
+        fWaves.WaveSbInfos.clear();
+        // Yes, see bitstream code
+        fWaves.WaveSbInfos.resize(usedContiguousSb[leader]);
+        folowerIt = data[!leader].GhaInfos.begin();
+
+        for (size_t i = 0; i < usedContiguousSb[leader]; i++) {
+            ResultBuf.SecondChBands[i] = false;
+        }
+    }
+
+    const auto& ghaInfos = data[leader].GhaInfos;
+    TWavesChannel& waves = ResultBuf.Waves[0];
+    waves.WaveParams.clear();
+    waves.WaveSbInfos.clear();
+    waves.WaveSbInfos.resize(usedContiguousSb[leader]);
+    auto it = ghaInfos.begin();
+    uint32_t nextSb = 0;
+    uint32_t index = 0;
+
+    uint32_t nextFolowerSb = 1;
+    while (it != ghaInfos.end()) {
+        const auto sb = ((it->first) >> 10);
+        if (sb >= usedContiguousSb[leader]) {
+            break;
+        }
+
+        const auto freqIndex = it->first & 1023;
+        const auto phaseIndex = GhaPhaseToIndex(it->second.phase);
+        const auto ampSf = AmplitudeToSf(it->second.magnitude);
+
+        waves.WaveSbInfos[sb].WaveNums++;
+        if (sb != nextSb) {
+            // cur subband done
+            // update index sb -> wave position index
+            waves.WaveSbInfos[sb].WaveIndex = index;
+
+            // process folower if present
+            if (data.size() == 2) {
+                nextFolowerSb = FillFolowerRes(data[!leader].GhaInfos, folowerIt, nextSb);
+            }
+
+            nextSb = sb;
+        }
+        waves.WaveParams.push_back(TAt3PGhaData::TWaveParam{freqIndex, ampSf, 1, phaseIndex});
+        it++;
+        index++;
+    }
+
+    if (data.size() == 2 && nextFolowerSb <= usedContiguousSb[leader]) {
+        FillFolowerRes(data[!leader].GhaInfos, folowerIt, nextSb);
+    }
+}
+
+uint32_t TGhaProcessor::FillFolowerRes(const TGhaInfoMap& fGhaInfos, TGhaInfoMap::const_iterator& it, const uint32_t leaderSb)
+{
+    TWavesChannel& waves = ResultBuf.Waves[1];
+    uint32_t fSb;
+    if (it != fGhaInfos.end()) {
+        fSb = ((it->first) >> 10);
+        ResultBuf.SecondChBands[fSb] = true;
+        waves.WaveSbInfos[fSb].WaveIndex = waves.WaveParams.size();
+    }
+
+    while (it != fGhaInfos.end()) {
+        fSb = ((it->first) >> 10);
+        if (fSb > leaderSb) {
+            return fSb;
+        }
+
+        const auto freqIndex = it->first & 1023;
+        const auto phaseIndex = GhaPhaseToIndex(it->second.phase);
+        const auto ampSf = AmplitudeToSf(it->second.magnitude);
+
+        waves.WaveSbInfos[fSb].WaveNums++;
+        waves.WaveParams.push_back(TAt3PGhaData::TWaveParam{freqIndex, ampSf, 1, phaseIndex});
+
+        it++;
+    }
+    return 0;
+}
+
+uint32_t TGhaProcessor::AmplitudeToSf(float amp)
+{
+    auto it = std::upper_bound(AmpSfTab.begin(), AmpSfTab.end(), amp);
+    if (it != AmpSfTab.begin()) {
+        it--;
+    }
+    return it - AmpSfTab.begin();
+}
+
+} // namespace
+
+std::unique_ptr<IGhaProcessor> MakeGhaProcessor0(float* b1, float* b2)
+{
+    return std::unique_ptr<TGhaProcessor>(new TGhaProcessor(b1, b2));
+}
+
+} // namespace NAtracDEnc
diff --git a/src/atrac/at3p/at3p_gha.h b/src/atrac/at3p/at3p_gha.h
index 4625ac3..951ece7 100644
--- a/src/atrac/at3p/at3p_gha.h
+++ b/src/atrac/at3p/at3p_gha.h
@@ -1,7 +1,24 @@
+/*
+ * This file is part of AtracDEnc.
+ *
+ * AtracDEnc 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.
+ *
+ * AtracDEnc 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 AtracDEnc; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
 #pragma once
 
 #include <config.h>
-#include <libgha/include/libgha.h>
 
 #include <vector>
 
@@ -14,7 +31,6 @@ struct TAt3PGhaData {
         uint32_t AmpSf;
         uint32_t AmpIndex;
         uint32_t PhaseIndex;
-        struct gha_info orig_gha_info;
     };
     struct TWaveSbInfo {
         size_t WaveIndex;
@@ -27,7 +43,7 @@ struct TAt3PGhaData {
     std::array<TWavesChannel, 2> Waves;
 
     uint8_t NumToneBands;
-    uint8_t SecondChBands[16];
+    bool SecondChBands[16];
 
     size_t GetNumWaves(size_t ch, size_t sb) const {
         return Waves[ch].WaveSbInfos.at(sb).WaveNums;
@@ -39,5 +55,13 @@ struct TAt3PGhaData {
     }
 };
 
+class IGhaProcessor {
+public:
+    virtual ~IGhaProcessor() {}
+    virtual const TAt3PGhaData* DoAnalize() = 0;
+};
+
+std::unique_ptr<IGhaProcessor> MakeGhaProcessor0(float* b1, float* b2);
+
 } // namespace NAtracDEnc
 
diff --git a/src/atrac/at3p/at3p_gha_ut.cpp b/src/atrac/at3p/at3p_gha_ut.cpp
new file mode 100644
index 0000000..718bb04
--- /dev/null
+++ b/src/atrac/at3p/at3p_gha_ut.cpp
@@ -0,0 +1,359 @@
+#include "at3p_gha.h"
+#include <gtest/gtest.h>
+#include <cmath>
+#include <vector>
+
+using std::vector;
+
+using namespace NAtracDEnc;
+
+struct TTestParam {
+    float freq;
+    uint16_t amplitude;
+    uint16_t start;
+    uint16_t end;
+};
+
+static void __attribute__ ((noinline)) Gen(const TTestParam& p, vector<float>& out)
+{
+    float freq = p.freq / (44100.0 / 16.0);
+    float a = p.amplitude;
+    int end = p.end;
+    for (int i = p.start; i < end; i++) {
+        out[i] += sin(freq * (float)(i % 128) * 2 * M_PI) * a;
+    }
+}
+
+static const TAt3PGhaData __attribute__ ((noinline)) GenAndRunGha(vector<TTestParam> p1, vector<TTestParam> p2)
+{
+    vector<float> buf1(2048);
+
+    for (const auto& p : p1) {
+        Gen(p, buf1);
+    }
+
+    vector<float> buf2;
+
+    if (!p2.empty()) {
+        buf2.resize(2048);
+
+        for (const auto& p : p2) {
+            Gen(p, buf2);
+        }
+    }
+
+    auto processor = MakeGhaProcessor0(buf1.data(), buf2.empty() ? nullptr : buf2.data());
+
+
+    return *(processor->DoAnalize());
+}
+
+// Single channel simple cases
+
+TEST(AT3PGHA, 689hz0625__full_frame_mono) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 0, 128}}, {});
+    EXPECT_EQ(res.NumToneBands, 1);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 1);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first->FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 0).first->AmpSf, 63);
+}
+
+
+TEST(AT3PGHA, 689hz0625__partial_frame_mono) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 32, 128}}, {});
+    EXPECT_EQ(res.NumToneBands, 1);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 1);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first->FreqIndex, 512);
+}
+
+TEST(AT3PGHA, 689hz0625_900hz__full_frame_mono) {
+    auto res = GenAndRunGha({{689.0625f, 16384, 0, 128}, {900.0, 8192, 0, 128}}, {});
+    EXPECT_EQ(res.NumToneBands, 1);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 2);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 2);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 2);
+    EXPECT_EQ(res.GetWaves(0, 0).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 0).first[1].FreqIndex, 669);
+}
+
+TEST(AT3PGHA, 400hz_800hz__full_frame_mono) {
+    auto res = GenAndRunGha({{400.0, 16384, 0, 128}, {800.0, 4096, 0, 128}}, {});
+    EXPECT_EQ(res.NumToneBands, 1);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 2);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 2);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 2);
+    EXPECT_EQ(res.GetWaves(0, 0).first[0].FreqIndex, 297);
+    EXPECT_EQ(res.GetWaves(0, 0).first[1].FreqIndex, 594);
+}
+
+TEST(AT3PGHA, 689hz0625_2067hz1875__full_frame_mono) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 0, 128}, {689.0625f, 16384, 128, 256}}, {});
+    EXPECT_EQ(res.NumToneBands, 2);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 2);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 2);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 1), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 1).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 1).first[0].FreqIndex, 512);
+}
+
+TEST(AT3PGHA, 689hz0625_4823hz4375__full_frame_mono) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 0, 128}, {689.0625f, 16384, 256, 384}}, {});
+    EXPECT_EQ(res.NumToneBands, 1);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 1);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first->FreqIndex, 512);
+}
+
+// Two channels simple cases
+
+TEST(AT3PGHA, 689hz0625__full_frame_stereo) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 0, 128}}, {{689.0625f, 32768, 0, 128}});
+    EXPECT_EQ(res.NumToneBands, 1);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 1);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first->FreqIndex, 512);
+
+    EXPECT_EQ(res.Waves[1].WaveParams.size(), 1);
+    EXPECT_EQ(res.Waves[1].WaveSbInfos.size(), 1);
+    EXPECT_EQ(res.GetNumWaves(1, 0), 1);
+    EXPECT_EQ(res.GetWaves(1, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(1, 0).first->FreqIndex, 512);
+}
+
+TEST(AT3PGHA, 689hz0625__full_frame_stereo_multiple_second) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 0, 128}}, {{689.0625f, 16384, 0, 128}, {900.0, 8192, 0, 128}});
+    EXPECT_EQ(res.NumToneBands, 1);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 1);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first->FreqIndex, 512);
+
+    EXPECT_EQ(res.Waves[1].WaveParams.size(), 2);
+    EXPECT_EQ(res.Waves[1].WaveSbInfos.size(), 1);
+    EXPECT_EQ(res.GetNumWaves(1, 0), 2);
+    EXPECT_EQ(res.GetWaves(1, 0).second, 2);
+    EXPECT_EQ(res.GetWaves(1, 0).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(1, 0).first[1].FreqIndex, 669);
+}
+
+TEST(AT3PGHA, 689hz0625_2067hz1875__full_frame_stereo_second_is_leader) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 0, 128}},
+                            {{689.0625f, 32768, 0, 128}, {689.0625f, 16384, 128, 256}});
+    EXPECT_EQ(res.NumToneBands, 2);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 2);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 2);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 1), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 1).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 1).first[0].FreqIndex, 512);
+
+    EXPECT_EQ(res.Waves[1].WaveParams.size(), 1);
+    EXPECT_EQ(res.Waves[1].WaveSbInfos.size(), 2);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetWaves(1, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(1, 0).first[0].FreqIndex, 512);
+}
+
+TEST(AT3PGHA, 689hz0625_2067hz1875_3445hz3125__full_frame_stereo_folower_has_0_2) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 0, 128}, {689.0625f, 32768, 128, 256}, {689.0625f, 16384, 256, 384}},
+                            {{689.0625f, 32768, 0, 128}, {689.0625f, 16384, 256, 384}});
+    EXPECT_EQ(res.NumToneBands, 3);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 3);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 3);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 1), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 2), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 1).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 1).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 2).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 2).first[0].FreqIndex, 512);
+
+    EXPECT_EQ(res.SecondChBands[0], true);
+    EXPECT_EQ(res.SecondChBands[1], false);
+    EXPECT_EQ(res.SecondChBands[2], true);
+    EXPECT_EQ(res.Waves[1].WaveParams.size(), 2);
+    EXPECT_EQ(res.Waves[1].WaveSbInfos.size(), 3);
+    EXPECT_EQ(res.GetNumWaves(1, 0), 1);
+    EXPECT_EQ(res.GetNumWaves(1, 2), 1);
+    EXPECT_EQ(res.GetWaves(1, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(1, 0).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(1, 2).second, 1);
+    EXPECT_EQ(res.GetWaves(1, 2).first[0].FreqIndex, 512);
+}
+
+TEST(AT3PGHA, 689hz0625_2067hz1875_3445hz3125__full_frame_stereo_folower_has_2) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 0, 128}, {689.0625f, 32768, 128, 256}, {689.0625f, 16384, 256, 384}},
+                            {{689.0625f, 16384, 256, 384}});
+    EXPECT_EQ(res.NumToneBands, 3);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 3);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 3);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 1), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 2), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 1).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 1).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 2).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 2).first[0].FreqIndex, 512);
+
+    EXPECT_EQ(res.SecondChBands[0], false);
+    EXPECT_EQ(res.SecondChBands[1], false);
+    EXPECT_EQ(res.SecondChBands[2], true);
+    EXPECT_EQ(res.Waves[1].WaveParams.size(), 1);
+    EXPECT_EQ(res.Waves[1].WaveSbInfos.size(), 3);
+    EXPECT_EQ(res.GetNumWaves(1, 2), 1);
+    EXPECT_EQ(res.GetWaves(1, 2).second, 1);
+    EXPECT_EQ(res.GetWaves(1, 2).first[0].FreqIndex, 512);
+}
+
+TEST(AT3PGHA, 689hz0625_2067hz1875_3445hz3125__full_frame_stereo_folower_has_1) {
+    auto res = GenAndRunGha({{689.0625f, 32768, 0, 128}, {689.0625f, 32768, 128, 256}, {689.0625f, 16384, 256, 384}},
+                            {{689.0625f, 16384, 128, 256}});
+    EXPECT_EQ(res.NumToneBands, 3);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 3);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 3);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 1), 1);
+    EXPECT_EQ(res.GetNumWaves(0, 2), 1);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 0).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 1).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 1).first[0].FreqIndex, 512);
+    EXPECT_EQ(res.GetWaves(0, 2).second, 1);
+    EXPECT_EQ(res.GetWaves(0, 2).first[0].FreqIndex, 512);
+
+    EXPECT_EQ(res.SecondChBands[0], false);
+    EXPECT_EQ(res.SecondChBands[1], true);
+    EXPECT_EQ(res.SecondChBands[2], false);
+    EXPECT_EQ(res.Waves[1].WaveParams.size(), 1);
+    EXPECT_EQ(res.Waves[1].WaveSbInfos.size(), 3);
+    EXPECT_EQ(res.GetNumWaves(1, 1), 1);
+    EXPECT_EQ(res.GetWaves(1, 1).second, 1);
+    EXPECT_EQ(res.GetWaves(1, 1).first[0].FreqIndex, 512);
+}
+
+TEST(AT3PGHA, max_tones_multiple_bands_full_frame_stereo) {
+    auto res = GenAndRunGha({
+                                {60.0f, 8192, 0,   128}, {120.0f, 8192, 0,   128}, {180.0f, 4096, 0,   128}, {240.0f, 2048, 0,   128},
+                                {60.0f, 8192, 128, 256}, {120.0f, 8192, 128, 256}, {180.0f, 4096, 128, 256}, {240.0f, 2048, 128, 256},
+                                {60.0f, 8192, 256, 384}, {120.0f, 8192, 256, 384}, {180.0f, 4096, 256, 384}, {240.0f, 2048, 256, 384},
+                                {60.0f, 8192, 384, 512}, {120.0f, 8192, 384, 512}, {180.0f, 4096, 384, 512}, {240.0f, 2048, 384, 512},
+                                {60.0f, 8192, 512, 640}, {120.0f, 8192, 512, 640}, {180.0f, 4096, 512, 640}, {240.0f, 2048, 512, 640},
+                                {60.0f, 8192, 640, 768}, {120.0f, 8192, 640, 768}, {180.0f, 4096, 640, 768}, {240.0f, 2048, 640, 768},
+                                {60.0f, 8192, 768, 896}, {120.0f, 8192, 768, 896}, {180.0f, 4096, 768, 896}, {240.0f, 2048, 768, 896},
+                            }, {
+                                {60.0f, 8192, 0,   128}, {120.0f, 8192, 0,   128}, {180.0f, 4096, 0,   128}, {240.0f, 2048, 0,   128},
+                                {60.0f, 8192, 128, 256}, {120.0f, 8192, 128, 256}, {180.0f, 4096, 128, 256}, {240.0f, 2048, 128, 256},
+                                {60.0f, 8192, 256, 384}, {120.0f, 8192, 256, 384}, {180.0f, 4096, 256, 384}, {240.0f, 2048, 256, 384},
+                                {60.0f, 8192, 384, 512}, {120.0f, 8192, 384, 512}, {180.0f, 4096, 384, 512}, {240.0f, 2048, 384, 512},
+                                {60.0f, 8192, 512, 640}, {120.0f, 8192, 512, 640}, {180.0f, 4096, 512, 640}, {240.0f, 2048, 512, 640},
+                                {60.0f, 8192, 640, 768}, {120.0f, 8192, 640, 768}, {180.0f, 4096, 640, 768}, {240.0f, 2048, 640, 768},
+                                {60.0f, 8192, 768, 896}, {120.0f, 8192, 768, 896}, {180.0f, 4096, 768, 896}, {240.0f, 2048, 768, 896},
+                            });
+    EXPECT_EQ(res.NumToneBands, 7);
+    EXPECT_EQ(res.Waves[0].WaveParams.size(), 28);
+    EXPECT_EQ(res.Waves[0].WaveSbInfos.size(), 7);
+    EXPECT_EQ(res.GetNumWaves(0, 0), 4);
+    EXPECT_EQ(res.GetNumWaves(0, 1), 4);
+    EXPECT_EQ(res.GetNumWaves(0, 2), 4);
+    EXPECT_EQ(res.GetNumWaves(0, 3), 4);
+    EXPECT_EQ(res.GetNumWaves(0, 4), 4);
+    EXPECT_EQ(res.GetNumWaves(0, 5), 4);
+    EXPECT_EQ(res.GetNumWaves(0, 6), 4);
+    EXPECT_EQ(res.GetWaves(0, 0).second, 4);
+    EXPECT_EQ(res.GetWaves(0, 0).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(0, 0).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(0, 0).first[2].FreqIndex, 134);
+    EXPECT_EQ(res.GetWaves(0, 0).first[3].FreqIndex, 178);
+
+    EXPECT_EQ(res.GetWaves(0, 1).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(0, 1).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(0, 1).first[2].FreqIndex, 134);
+    EXPECT_EQ(res.GetWaves(0, 1).first[3].FreqIndex, 178);
+
+    EXPECT_EQ(res.GetWaves(0, 2).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(0, 2).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(0, 2).first[2].FreqIndex, 134);
+    EXPECT_EQ(res.GetWaves(0, 2).first[3].FreqIndex, 178);
+
+    EXPECT_EQ(res.GetWaves(0, 3).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(0, 3).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(0, 3).first[2].FreqIndex, 134);
+    EXPECT_EQ(res.GetWaves(0, 3).first[3].FreqIndex, 178);
+
+    EXPECT_EQ(res.GetWaves(0, 4).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(0, 4).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(0, 4).first[2].FreqIndex, 134);
+    EXPECT_EQ(res.GetWaves(0, 4).first[3].FreqIndex, 178);
+
+    EXPECT_EQ(res.GetWaves(0, 5).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(0, 5).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(0, 5).first[2].FreqIndex, 134);
+    EXPECT_EQ(res.GetWaves(0, 5).first[3].FreqIndex, 178);
+
+    EXPECT_EQ(res.GetWaves(0, 6).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(0, 6).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(0, 6).first[2].FreqIndex, 134);
+    EXPECT_EQ(res.GetWaves(0, 6).first[3].FreqIndex, 178);
+
+    EXPECT_EQ(res.Waves[1].WaveParams.size(), 21);
+    EXPECT_EQ(res.Waves[1].WaveSbInfos.size(), 7);
+    EXPECT_EQ(res.GetNumWaves(1, 0), 3);
+    EXPECT_EQ(res.GetNumWaves(1, 1), 3);
+    EXPECT_EQ(res.GetNumWaves(1, 2), 3);
+    EXPECT_EQ(res.GetNumWaves(1, 3), 3);
+    EXPECT_EQ(res.GetNumWaves(1, 4), 3);
+    EXPECT_EQ(res.GetNumWaves(1, 5), 3);
+    EXPECT_EQ(res.GetNumWaves(1, 6), 3);
+    EXPECT_EQ(res.GetWaves(1, 0).second, 3);
+    EXPECT_EQ(res.GetWaves(1, 0).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(1, 0).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(1, 0).first[2].FreqIndex, 134);
+
+    EXPECT_EQ(res.GetWaves(1, 1).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(1, 1).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(1, 1).first[2].FreqIndex, 134);
+
+    EXPECT_EQ(res.GetWaves(1, 2).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(1, 2).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(1, 2).first[2].FreqIndex, 134);
+
+    EXPECT_EQ(res.GetWaves(1, 3).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(1, 3).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(1, 3).first[2].FreqIndex, 134);
+
+    EXPECT_EQ(res.GetWaves(1, 4).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(1, 4).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(1, 4).first[2].FreqIndex, 134);
+
+    EXPECT_EQ(res.GetWaves(1, 5).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(1, 5).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(1, 5).first[2].FreqIndex, 134);
+
+    EXPECT_EQ(res.GetWaves(1, 6).first[0].FreqIndex, 45);
+    EXPECT_EQ(res.GetWaves(1, 6).first[1].FreqIndex, 89);
+    EXPECT_EQ(res.GetWaves(1, 6).first[2].FreqIndex, 134);
+}
+
+
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 66218d8..f6d1827 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -71,7 +71,26 @@ target_link_libraries(at3plus_bitstream_ut
     GTest::gtest_main
 )
 
+###
+
+set(at3plus_gha_ut
+    ${CMAKE_SOURCE_DIR}/src/atrac/at3p/at3p_gha_ut.cpp
+)
+
+add_executable(at3plus_gha_ut ${at3plus_gha_ut})
+
+target_link_libraries(at3plus_gha_ut
+    m
+    fft_impl
+    atracdenc_impl
+    GTest::gtest_main
+)
+
+###
+
+
 enable_testing()
 add_test(ut atracdenc_ut)
 add_test(at3plus_pqf_ut at3plus_pqf_ut)
 add_test(at3plus_bitstream_ut at3plus_bitstream_ut)
+add_test(at3plus_gha_ut at3plus_gha_ut)