Merge pull request #35 from Digitelektro/develop

Develop

Author: Digitelektro

CMakeLists.txt

@@ -4,7 +4,7 @@ include(ExternalProject)
 
 project(Meteordemod LANGUAGES CXX)
 
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pthread")
 
@@ -93,6 +93,12 @@ add_executable(meteordemod
     DSP/meteordemodulator.cpp
     DSP/agc.cpp
     DSP/pll.cpp
+	DSP/meteorcostas.cpp
+	DSP/phasecontrolloop.cpp
+	DSP/polyphasebank.cpp
+	DSP/mm.cpp
+	DSP/window.cpp
+	DSP/windowedsinc.cpp
     DSP/filter.cpp
     DSP/iqsource.cpp
     DSP/wavreader.cpp

DSP/agc.cpp

@@ -2,32 +2,16 @@
 
 namespace DSP {
 
-Agc::Agc()
-    : mWindowSize(AGC_WINSIZE)
-    , mAvg(AGC_TARGET)
+Agc::Agc(float targetAmplitude, float maxGain, float windowSize, float biasWindowSize)
+    : mWindowSize(windowSize)
+    , mAvg(targetAmplitude)
     , mGain(1)
-    , mTargetAmplitude(AGC_TARGET)
+    , mMaxGain(maxGain)
+    , mTargetAmplitude(targetAmplitude)
+    , mBiasWindowSize(biasWindowSize)
     , mBias(0)
 {
 
 }
 
-Agc::complex Agc::process(complex sample)
-{
-    float rho;
-
-    mBias = (mBias * static_cast<float>(AGC_BIAS_WINSIZE - 1) + sample) / static_cast<float>(AGC_BIAS_WINSIZE);
-    sample -= mBias;
-
-    // Update the sample magnitude average
-    rho = sqrtf(std::real(sample) * std::real(sample) + std::imag(sample) * std::imag(sample));
-    mAvg = (mAvg * (mWindowSize - 1) + rho) / mWindowSize;
-
-    mGain = mTargetAmplitude / mAvg;
-    if (mGain > AGC_MAX_GAIN) {
-        mGain = AGC_MAX_GAIN;
-    }
-    return sample * mGain;
-}
-
 } //namespace DSP

DSP/agc.h

@@ -11,9 +11,31 @@ class Agc
     typedef std::complex<float> complex;
 
 public:
-    Agc();
+    Agc(float targetAmplitude, float maxGain = 20, float windowSize = 1024*64, float biasWindowSize = 256*1024);
 
-    complex process(complex sample);
+    inline complex process(complex sample){
+        float rho;
+
+        mBias = (mBias * static_cast<float>(mBiasWindowSize - 1) + sample) / static_cast<float>(mBiasWindowSize);
+        sample -= mBias;
+
+        // Update the sample magnitude average
+        rho = std::abs(sample);
+        mAvg = (mAvg * (mWindowSize - 1) + rho) / mWindowSize;
+
+        mGain = mTargetAmplitude / mAvg;
+        if (mGain > mMaxGain) {
+            mGain = mMaxGain;
+        }
+        return sample * mGain;
+    }
+
+    inline void process(const complex *inSamples, complex *outsamples, unsigned int count) {
+        for(unsigned int i = 0; i < count; i++) {
+            *outsamples = process(*inSamples++);
+            outsamples++;
+        }
+    }
 
 public:
     float getGain() const {
@@ -24,14 +46,10 @@ class Agc
     uint32_t mWindowSize;
     float mAvg;
     float mGain;
+    float mMaxGain;
     float mTargetAmplitude;
+    float mBiasWindowSize;
     complex mBias;
-
-private:
-    static constexpr uint32_t AGC_WINSIZE = 1024*64;
-    static constexpr uint32_t AGC_TARGET = 180;
-    static constexpr uint32_t AGC_MAX_GAIN = 20;
-    static constexpr uint32_t AGC_BIAS_WINSIZE = 256*1024;
 };
 
 } //namespace DSP

DSP/filter.cpp

@@ -9,49 +9,31 @@ FilterBase::FilterBase(int taps)
 
 }
 
-FilterBase::complex FilterBase::process(const FilterBase::complex &in)
+RRCFilter::RRCFilter(int taps, float beta, float symbolrate, float samplerate)
+    : FilterBase(taps)
 {
-    complex out = 0.0f;
-
-    mMemory.push_front(in);
-    mMemory.pop_back();
-
-    std::list<complex>::const_reverse_iterator it = mMemory.crbegin();
-
-    for(int i = mTaps - 1; i >=0; --i, ++it) {
-        out += (*it) * mCoeffs[i];
-    }
-    return out;
-}
-
-RRCFilter::RRCFilter(unsigned order, unsigned factor, float osf, float alpha)
-    :FilterBase(order * 2 + 1)
-{
-    for (uint16_t i=0; i < mTaps; i++) {
-        mCoeffs.emplace_back(computeCoeffs(i, mTaps, osf * factor, alpha));
-    }
-
+    computeCoeffs(taps, beta, samplerate / symbolrate);
     mMemory.resize(mTaps, 0);
 }
 
-float RRCFilter::computeCoeffs(int stageNo, uint16_t taps, float osf, float alpha)
+void RRCFilter::computeCoeffs(int taps, float beta, float Ts)
 {
     float coeff;
-    float t;
-    float interm;
-    int order;
-
-    order = (taps - 1)/2;
-
-    if (order == stageNo) {
-        return 1 - alpha + 4 * alpha / M_PI;
+    float half = taps / 2.0f;
+    float limit = Ts / (4.0 * beta);
+    for (int i = 0; i < taps; i++) {
+        float t = (float)i - half + 0.5;
+        if (t == 0.0) {
+            coeff = (1.0f + beta * (4.0f / M_PI - 1.0f)) / Ts;
+        }
+        else if (t == limit || t == -limit) {
+            coeff = ((1.0f + 2.0f / M_PI) * sinf(M_PI / (4.0f * beta)) + (1.0f - 2.0f / M_PI)*cos(M_PI / (4.0f * beta))) * beta / (Ts * sqrtf(2.0f));
+        }
+        else {
+            coeff = ((sinf((1.0f - beta) * M_PI * t / Ts) + cosf((1.0f + beta) * M_PI * t / Ts) * 4.0f * beta * t / Ts) / ((1.0f - (4.0f * beta * t / Ts) * (4.0f * beta * t / Ts)) * M_PI * t / Ts)) / Ts;
+        }
+        mCoeffs.emplace_back(coeff);
     }
-
-    t = std::abs(order - stageNo) / osf;
-    coeff = std::sin(M_PI * t * (1 - alpha)) + 4 * alpha * t * std::cos(M_PI * t * (1 + alpha));
-    interm = M_PI * t * (1 - (4 * alpha * t) * (4 * alpha * t));
-
-    return coeff / interm;
 }
 
 } //namespace DSP

DSP/filter.h

@@ -14,7 +14,27 @@ class FilterBase
 public:
     FilterBase(int taps);
     virtual ~FilterBase() {}
-    complex process(const complex &in);
+
+    inline complex process(const complex &in) {
+        complex out = 0.0f;
+
+        mMemory.push_front(in);
+        mMemory.pop_back();
+
+        std::list<complex>::const_reverse_iterator it = mMemory.crbegin();
+
+        for(int i = mTaps - 1; i >=0; --i, ++it) {
+            out += (*it) * mCoeffs[i];
+        }
+        return out;
+    }
+
+    inline void process(const complex *inSamples, complex *outSamples, unsigned int count) {
+        for(unsigned int i = 0; i < count; i++) {
+            *outSamples = process(*inSamples++);
+            outSamples++;
+        }
+    }
 
 protected:
     std::vector<float> mCoeffs;
@@ -24,10 +44,10 @@ class FilterBase
 
 class RRCFilter : public FilterBase {
 public:
-    RRCFilter(unsigned order, unsigned factor, float osf, float alpha);
+    RRCFilter(int taps, float beta, float symbolrate, float samplerate);
 
 private:
-    float computeCoeffs(int stageNo, uint16_t taps, float osf, float alpha);
+    void computeCoeffs(int taps, float beta, float Ts);
 };
 
 } //namespace DSP

DSP/meteorcostas.cpp

@@ -0,0 +1,17 @@
+#include "meteorcostas.h"
+
+namespace DSP {
+
+MeteorCostas::MeteorCostas(float bandWidth, float initPhase, float initFreq, float minFreq, float maxFreq, bool brokenModulation)
+    : PLL(bandWidth, initPhase, initFreq, minFreq, maxFreq)
+    , mPllOriginalBandwidth(bandWidth)
+    , mBrokenModulation(brokenModulation)
+    , mError(0.0f)
+    , mLockDetector(0.0f)
+    , mIsLocked(false)
+    , mIsLockedOnce(false)
+{
+
+}
+
+} // namespace DSP

DSP/meteorcostas.h

@@ -0,0 +1,99 @@
+#ifndef DSP_METEORCOSTAS_H
+#define DSP_METEORCOSTAS_H
+
+#include "pll.h"
+#include <algorithm>
+
+namespace DSP {
+
+class MeteorCostas : public PLL
+{    
+private:
+    static constexpr float cLockDetectionTreshold = 0.22;
+    static constexpr float cUnLockDetectionTreshold = 0.25;
+    static constexpr float cLockFilterCoeff = 0.00001f;
+
+public:
+    MeteorCostas(float bandWidth, float initPhase = 0.0f, float initFreq = 0.0f, float minFreq = -M_PI, float maxFreq = M_PI, bool brokenModulation = false);
+
+    inline virtual complex process(const complex &sample) override {
+        complex retval;
+        retval = sample * complex(cosf(-mPhase), sinf(-mPhase));
+        mError = errorFunction(retval);
+        advance(mError);
+        return retval;
+    }
+
+    inline virtual void process(const complex *insamples, complex *outsampes, unsigned int count) {
+        for(unsigned int i = 0; i < count; i++) {
+            *outsampes++ = process(*insamples++);
+        }
+    }
+
+protected:
+    float errorFunction(complex value) {
+        float error;
+        if(mBrokenModulation) {
+            const float PHASE1 = 0.47439988279190737;
+            const float PHASE2 = 2.1777839908413044;
+            const float PHASE3 = 3.8682349942715186;
+            const float PHASE4 = -0.29067248091319986;
+
+            float phase = atan2f(value.imag(), value.real());
+            float dp1 = normalizePhase(phase - PHASE1);
+            float dp2 = normalizePhase(phase - PHASE2);
+            float dp3 = normalizePhase(phase - PHASE3);
+            float dp4 = normalizePhase(phase - PHASE4);
+            float lowest = dp1;
+            if (fabsf(dp2) < fabsf(lowest)) { lowest = dp2; }
+            if (fabsf(dp3) < fabsf(lowest)) { lowest = dp3; }
+            if (fabsf(dp4) < fabsf(lowest)) { lowest = dp4; }
+            error = lowest * std::abs(value);
+        } else {
+            error = (step(value.real()) * value.imag()) - (step(value.imag()) * value.real());
+        }
+
+        mLockDetector = std::abs(error) * cLockFilterCoeff + mLockDetector * (1.0f - cLockFilterCoeff);
+
+        if(mLockDetector < cLockDetectionTreshold && !mIsLocked) {
+            mIsLocked = true;
+            setBandWidth(mPllOriginalBandwidth/5.0f);
+        } else if(mLockDetector > cUnLockDetectionTreshold && mIsLocked) {
+            mIsLocked = false;
+            setBandWidth(mPllOriginalBandwidth);
+        }
+
+        if(mLockDetector < cLockDetectionTreshold) {
+            mIsLockedOnce = true;
+        }
+
+        return std::clamp(error, -1.0f, 1.0f);
+    }
+
+    inline float step(float val) {
+        return val > 0 ? 1.0f : -1.0f;
+    }
+
+public:
+    inline float getError() const {
+        return mLockDetector;
+    }
+    inline bool isLocked() const {
+        return mIsLocked;
+    }
+    inline bool isLockedOnce() const {
+        return mIsLockedOnce;
+    }
+
+protected:
+    float mPllOriginalBandwidth;
+    bool mBrokenModulation;
+    float mError;
+    float mLockDetector;
+    bool mIsLocked;
+    bool mIsLockedOnce;
+};
+
+} // namespace DSP
+
+#endif // DSP_METEORCOSTAS_H