ggwave : more resampling fixes

This time the sound cracking should be fixed for real.

Also adding option to generate noise in the cpp tests

include/ggwave/ggwave.h

@@ -294,9 +294,17 @@ public:
     bool init(int dataSize, const char * dataBuffer, const TxProtocol & txProtocol, const int volume = kDefaultVolume);
 
     // expected waveform size of the encoded Tx data in bytes
+    //
+    //   When the output sampling rate is not equal to kBaseSampleRate the result of this method is overestimation of
+    //   the actual number of bytes that would be produced
+    //
     uint32_t encodeSize_bytes() const;
 
     // expected waveform size of the encoded Tx data in samples
+    //
+    //   When the output sampling rate is not equal to kBaseSampleRate the result of this method is overestimation of
+    //   the actual number of samples that would be produced
+    //
     uint32_t encodeSize_samples() const;
 
     // encode Tx data into an audio waveform

src/ggwave.cpp

@@ -473,8 +473,9 @@ uint32_t GGWave::encodeSize_samples() const {
     float factor = 1.0f;
     int samplesPerFrameOut = m_samplesPerFrame;
     if (m_sampleRateOut != kBaseSampleRate) {
-        factor = kBaseSampleRate/m_sampleRateOut;
+        factor = float(kBaseSampleRate)/m_sampleRateOut;
-        samplesPerFrameOut = m_impl->resampler.resample(factor, m_samplesPerFrame, m_outputBlock.data(), nullptr);
+        // note : +1 extra sample in order to overestimate the buffer size
+        samplesPerFrameOut = m_impl->resampler.resample(factor, m_samplesPerFrame, m_outputBlock.data(), nullptr) + 1;
     }
     int nECCBytesPerTx = getECCBytesForLength(m_txDataLength);
     int sendDataLength = m_txDataLength + m_encodedDataOffset;
@@ -489,6 +490,8 @@ uint32_t GGWave::encodeSize_samples() const {
 bool GGWave::encode(const CBWaveformOut & cbWaveformOut) {
     int frameId = 0;
 
+    m_impl->resampler.reset();
+
     std::vector<double> phaseOffsets(kMaxDataBits);
 
     for (int k = 0; k < (int) phaseOffsets.size(); ++k) {
@@ -540,10 +543,7 @@ bool GGWave::encode(const CBWaveformOut & cbWaveformOut) {
     rsData.Encode(m_txData.data() + 1, m_txDataEncoded.data() + m_encodedDataOffset);
 
     float factor = float(kBaseSampleRate)/m_sampleRateOut;
-    int samplesPerFrameOut = m_samplesPerFrame;
+    uint32_t offset = 0;
-    if (m_sampleRateOut != kBaseSampleRate) {
-        samplesPerFrameOut = m_impl->resampler.resample(factor, m_samplesPerFrame, m_outputBlock.data(), m_outputBlockResampled.data());
-    }
 
     while (m_hasNewTxData) {
         std::fill(m_outputBlock.begin(), m_outputBlock.end(), 0.0f);
@@ -610,14 +610,13 @@ bool GGWave::encode(const CBWaveformOut & cbWaveformOut) {
             m_outputBlock[i] *= scale;
         }
 
-        if (samplesPerFrameOut != m_samplesPerFrame) {
+        int samplesPerFrameOut = m_samplesPerFrame;
-            m_impl->resampler.resample(factor, m_samplesPerFrame, m_outputBlock.data(), m_outputBlockResampled.data());
+        if (m_sampleRateOut != kBaseSampleRate) {
+            samplesPerFrameOut = m_impl->resampler.resample(factor, m_samplesPerFrame, m_outputBlock.data(), m_outputBlockResampled.data());
         } else {
             m_outputBlockResampled = m_outputBlock;
         }
 
-        uint32_t offset = frameId*samplesPerFrameOut;
-
         // default output is in 16-bit signed int so we always compute it
         for (int i = 0; i < samplesPerFrameOut; ++i) {
             m_outputBlockI16[offset + i] = 32768*m_outputBlockResampled[i];
@@ -665,25 +664,26 @@ bool GGWave::encode(const CBWaveformOut & cbWaveformOut) {
         }
 
         ++frameId;
+        offset += samplesPerFrameOut;
     }
 
     switch (m_sampleFormatOut) {
         case GGWAVE_SAMPLE_FORMAT_UNDEFINED: break;
         case GGWAVE_SAMPLE_FORMAT_I16:
             {
-                cbWaveformOut(m_outputBlockI16.data(), frameId*samplesPerFrameOut*m_sampleSizeBytesOut);
+                cbWaveformOut(m_outputBlockI16.data(), offset*m_sampleSizeBytesOut);
             } break;
         case GGWAVE_SAMPLE_FORMAT_U8:
         case GGWAVE_SAMPLE_FORMAT_I8:
         case GGWAVE_SAMPLE_FORMAT_U16:
         case GGWAVE_SAMPLE_FORMAT_F32:
             {
-                cbWaveformOut(m_outputBlockTmp.data(), frameId*samplesPerFrameOut*m_sampleSizeBytesOut);
+                cbWaveformOut(m_outputBlockTmp.data(), offset*m_sampleSizeBytesOut);
             } break;
     }
 
-    m_txAmplitudeDataI16.resize(frameId*samplesPerFrameOut);
+    m_txAmplitudeDataI16.resize(offset);
-    for (int i = 0; i < frameId*samplesPerFrameOut; ++i) {
+    for (uint32_t i = 0; i < offset; ++i) {
         m_txAmplitudeDataI16[i] = m_outputBlockI16[i];
     }
 
@@ -697,7 +697,8 @@ void GGWave::decode(const CBWaveformInp & cbWaveformInp) {
         uint32_t nBytesNeeded = m_samplesNeeded*m_sampleSizeBytesInp;
 
         if (m_sampleRateInp != kBaseSampleRate) {
-            nBytesNeeded = m_impl->resampler.resample(1.0/factor, m_samplesNeeded, m_sampleAmplitudeResampled.data(), nullptr)*m_sampleSizeBytesInp;
+            // note : predict 4 extra samples just to make sure we have enough data
+            nBytesNeeded = (m_impl->resampler.resample(1.0f/factor, m_samplesNeeded, m_sampleAmplitudeResampled.data(), nullptr) + 4)*m_sampleSizeBytesInp;
         }
 
         uint32_t nBytesRecorded = 0;
@@ -770,13 +771,23 @@ void GGWave::decode(const CBWaveformInp & cbWaveformInp) {
             case GGWAVE_SAMPLE_FORMAT_F32: break;
         }
 
-        if (nBytesRecorded == 0) {
+        if (nSamplesRecorded == 0) {
             break;
         }
 
         uint32_t offset = m_samplesPerFrame - m_samplesNeeded;
 
         if (m_sampleRateInp != kBaseSampleRate) {
+            if (nSamplesRecorded <= 2*Resampler::kWidth) {
+                m_samplesNeeded = m_samplesPerFrame;
+                break;
+            }
+
+            // reset resampler state every minute
+            if (!m_receivingData && m_impl->resampler.nSamplesTotal() > 60.0f*factor*kBaseSampleRate) {
+                m_impl->resampler.reset();
+            }
+
             int nSamplesResampled = offset + m_impl->resampler.resample(factor, nSamplesRecorded, m_sampleAmplitudeResampled.data(), m_sampleAmplitude.data() + offset);
             nSamplesRecorded = nSamplesResampled;
         } else {

src/resampler.cpp

@@ -1,5 +1,6 @@
 #include "resampler.h"
 
+#include <cassert>
 #include <cmath>
 #include <cstdio>
 
@@ -9,8 +10,20 @@ double linear_interp(double first_number, double second_number, double fraction)
 }
 }
 
-Resampler::Resampler() {
+Resampler::Resampler() :
+    m_sincTable(kWidth*kSamplesPerZeroCrossing),
+    m_delayBuffer(3*kWidth),
+    m_edgeSamples(kWidth),
+    m_samplesInp(2048) {
     make_sinc();
+    reset();
+}
+
+void Resampler::reset() {
+    m_state = {};
+    std::fill(m_edgeSamples.begin(), m_edgeSamples.end(), 0.0f);
+    std::fill(m_delayBuffer.begin(), m_delayBuffer.end(), 0.0f);
+    std::fill(m_samplesInp.begin(), m_samplesInp.end(), 0.0f);
 }
 
 int Resampler::resample(
@@ -18,59 +31,96 @@ int Resampler::resample(
         int nSamples,
         const float * samplesInp,
         float * samplesOut) {
-    int idxInp = 0;
+    int idxInp = -1;
     int idxOut = 0;
     int notDone = 1;
-    double time_now = 0.0;
+    float data_in = 0.0f;
-    long num_samples = nSamples;
+    float data_out = 0.0f;
-    long int_time = 0;
-    long last_time = 0;
-    float data_in = samplesInp[idxInp];
-    float data_out;
     double one_over_factor = 1.0;
+
+    auto stateSave = m_state;
+
+    m_state.nSamplesTotal += nSamples;
+
+    if (samplesOut) {
+        assert(nSamples > kWidth);
+        if ((int) m_samplesInp.size() < nSamples + kWidth) {
+            m_samplesInp.resize(nSamples + kWidth);
+        }
+        for (int i = 0; i < kWidth; ++i) {
+            m_samplesInp[i] = m_edgeSamples[i];
+            m_edgeSamples[i] = samplesInp[nSamples - kWidth + i];
+        }
+        for (int i = 0; i < nSamples; ++i) {
+            m_samplesInp[i + kWidth] = samplesInp[i];
+        }
+        samplesInp = m_samplesInp.data();
+    }
+
     while (notDone) {
+        while (m_state.timeLast < m_state.timeInt) {
+            if (++idxInp >= nSamples) {
+                notDone = 0;
+                break;
+            } else {
+                data_in = samplesInp[idxInp];
+            }
+            //printf("xxxx idxInp = %d\n", idxInp);
+            if (samplesOut) new_data(data_in);
+            m_state.timeLast += 1;
+        }
+
+        if (notDone == false) break;
+
         double temp1 = 0.0;
-        long left_limit = time_now - kWidth + 1; /* leftmost neighboring sample used for interp.*/
+        int left_limit = m_state.timeNow - kWidth + 1; /* leftmost neighboring sample used for interp.*/
-        long right_limit = time_now + kWidth;    /* rightmost leftmost neighboring sample used for interp.*/
+        int right_limit = m_state.timeNow + kWidth;    /* rightmost leftmost neighboring sample used for interp.*/
-        if (left_limit<0) left_limit = 0;
+        if (left_limit < 0) left_limit = 0;
-        if (right_limit>num_samples) right_limit = num_samples;
+        if (right_limit > m_state.nSamplesTotal + kWidth) right_limit = m_state.nSamplesTotal + kWidth;
-        if (factor<1.0) {
+        if (factor < 1.0) {
-            for (int j=left_limit;j<right_limit;j++) {
+            for (int j = left_limit; j < right_limit; j++) {
-                temp1 += gimme_data(j-int_time)*sinc(time_now - (double) j);
+                temp1 += gimme_data(j - m_state.timeInt)*sinc(m_state.timeNow - (double) j);
             }
             data_out = temp1;
         }
         else {
             one_over_factor = 1.0 / factor;
-            for (int j=left_limit;j<right_limit;j++) {
+            for (int j = left_limit; j < right_limit; j++) {
-                temp1 += gimme_data(j-int_time)*one_over_factor*sinc(one_over_factor * (time_now - (double) j));
+                temp1 += gimme_data(j - m_state.timeInt)*one_over_factor*sinc(one_over_factor*(m_state.timeNow - (double) j));
             }
             data_out = temp1;
         }
 
         if (samplesOut) {
+            //printf("inp = %d, l = %d, r = %d, n = %d, a = %d, b = %d\n", idxInp, left_limit, right_limit, m_state.nSamplesTotal, left_limit - m_state.timeInt, right_limit - m_state.timeInt - 1);
             samplesOut[idxOut] = data_out;
         }
         ++idxOut;
 
-        time_now += factor;
+        m_state.timeNow += factor;
-        last_time = int_time;
+        m_state.timeLast = m_state.timeInt;
-        int_time = time_now;
+        m_state.timeInt = m_state.timeNow;
-        while (last_time<int_time) {
+        while (m_state.timeLast < m_state.timeInt) {
-            if (++idxInp == nSamples) {
+            if (++idxInp >= nSamples) {
                 notDone = 0;
+                break;
             } else {
                 data_in = samplesInp[idxInp];
             }
-            new_data(data_in);
+            if (samplesOut) new_data(data_in);
-            last_time += 1;
+            m_state.timeLast += 1;
         }
+        //printf("last idxInp = %d, nSamples = %d\n", idxInp, nSamples);
+    }
+
+    if (samplesOut == nullptr) {
+        m_state = stateSave;
     }
 
     return idxOut;
 }
 
-float Resampler::gimme_data(long j) const {
+float Resampler::gimme_data(int j) const {
     return m_delayBuffer[(int) j + kWidth];
 }
 
@@ -85,7 +135,7 @@ void Resampler::make_sinc() {
     double temp, win_freq, win;
     win_freq = M_PI/kWidth/kSamplesPerZeroCrossing;
     m_sincTable[0] = 1.0;
-    for (int i = 1; i < kWidth*kSamplesPerZeroCrossing;i++) {
+    for (int i = 1; i < kWidth*kSamplesPerZeroCrossing; i++) {
         temp = (double) i*M_PI/kSamplesPerZeroCrossing;
         m_sincTable[i] = sin(temp)/temp;
         win = 0.5 + 0.5*cos(win_freq*i);
@@ -99,7 +149,7 @@ double Resampler::sinc(double x) const {
     if (fabs(x) >= kWidth - 1) {
         return 0.0;
     } else {
-        temp = fabs(x) * (double) kSamplesPerZeroCrossing;
+        temp = fabs(x)*(double) kSamplesPerZeroCrossing;
         low = temp;          /* these are interpolation steps */
         delta = temp - low;  /* and can be ommited if desired */
         return linear_interp(m_sincTable[low], m_sincTable[low + 1], delta);

src/resampler.h

@@ -1,9 +1,21 @@
 #pragma once
 
+#include <vector>
+#include <cstdint>
+
 class Resampler {
 public:
+    // this controls the number of neighboring samples
+    // which are used to interpolate the new samples. The
+    // processing time is linearly related to this width
+    static const int kWidth = 64;
+
     Resampler();
 
+    void reset();
+
+    int nSamplesTotal() const { return m_state.nSamplesTotal; }
+
     int resample(
             float factor,
             int nSamples,
@@ -11,23 +23,27 @@ public:
             float * samplesOut);
 
 private:
-    float gimme_data(long j) const;
+    float gimme_data(int j) const;
     void new_data(float data);
     void make_sinc();
     double sinc(double x) const;
 
-    /* this controls the number of neighboring samples
-       which are used to interpolate the new samples.  The
-       processing time is linearly related to this width */
-    static const int kWidth = 64;
-
     static const int kDelaySize = 140;
 
-    /* this defines how finely the sinc function
+    // this defines how finely the sinc function is sampled for storage in the table
-       is sampled for storage in the table  */
     static const int kSamplesPerZeroCrossing = 32;
 
-    float m_sincTable[kWidth*kSamplesPerZeroCrossing] = { 0.0 };
+    std::vector<float> m_sincTable;
+    std::vector<float> m_delayBuffer;
+    std::vector<float> m_edgeSamples;
+    std::vector<float> m_samplesInp;
 
-    float m_delayBuffer[3*kWidth] = { 0 };
+    struct State {
+        int nSamplesTotal = 0;
+        int timeInt = 0;
+        int timeLast = 0;
+        double timeNow = 0.0;
+    };
+
+    State m_state;
 };

tests/test-ggwave.cpp

@@ -8,6 +8,8 @@
 #include <vector>
 #include <set>
 
+float frand() { return float(rand()%RAND_MAX)/RAND_MAX; }
+
 #define CHECK(cond) \
     if (!(cond)) { \
         fprintf(stderr, "[%s:%d] Check failed: %s\n", __FILE__, __LINE__, #cond); \
@@ -145,7 +147,42 @@ int main(int argc, char ** argv) {
                     };
                 } break;
         };
+    };
 
+    auto addNoiseHelper = [&](float level, GGWave::SampleFormat format) {
+        switch (format) {
+            case GGWAVE_SAMPLE_FORMAT_UNDEFINED: CHECK(false); break;
+            case GGWAVE_SAMPLE_FORMAT_U8:
+                {
+                    for (auto & s : bufferU8) {
+                        s = std::max(0.0f, std::min(255.0f, (float) s + (frand() - 0.5f)*(level*256)));
+                    }
+                } break;
+            case GGWAVE_SAMPLE_FORMAT_I8:
+                {
+                    for (auto & s : bufferI8) {
+                        s = std::max(-128.0f, std::min(127.0f, (float) s + (frand() - 0.5f)*(level*256)));
+                    }
+                } break;
+            case GGWAVE_SAMPLE_FORMAT_U16:
+                {
+                    for (auto & s : bufferU16) {
+                        s = std::max(0.0f, std::min(65535.0f, (float) s + (frand() - 0.5f)*(level*65536)));
+                    }
+                } break;
+            case GGWAVE_SAMPLE_FORMAT_I16:
+                {
+                    for (auto & s : bufferI16) {
+                        s = std::max(-32768.0f, std::min(32767.0f, (float) s + (frand() - 0.5f)*(level*65536)));
+                    }
+                } break;
+            case GGWAVE_SAMPLE_FORMAT_F32:
+                {
+                    for (auto & s : bufferF32) {
+                        s = std::max(-1.0f, std::min(1.0f, (float) s + (frand() - 0.5f)*(level)));
+                    }
+                } break;
+        };
     };
 
     uint32_t nSamples = 0;
@@ -194,7 +231,7 @@ int main(int argc, char ** argv) {
         printf("Testing: sample rate = %d\n", srInp);
 
         auto parameters = GGWave::getDefaultParameters();
-        parameters.soundMarkerThreshold = 1.1f;
+        parameters.soundMarkerThreshold = 3.0f;
 
         std::string payload = "hello123";
 
@@ -207,7 +244,8 @@ int main(int argc, char ** argv) {
             auto expectedSize = instanceOut.encodeSize_samples();
             instanceOut.encode(kCBWaveformOut.at(parameters.sampleFormatOut));
             printf("Expected = %d, actual = %d\n", expectedSize, nSamples);
-            CHECK(expectedSize == nSamples);
+            CHECK(expectedSize >= nSamples);
+            addNoiseHelper(0.01, parameters.sampleFormatOut); // add some artificial noise
             convertHelper(parameters.sampleFormatOut, parameters.sampleFormatInp);
         }