ggwave : add filter function

Currently support the following filter:

- Hann window
- Hamming window
- First-order high-pass

include/ggwave/ggwave.h

@@ -72,6 +72,12 @@ extern "C" {
         GGWAVE_PROTOCOL_COUNT,
     } ggwave_ProtocolId;
 
+    typedef enum {
+        GGWAVE_FILTER_HANN,
+        GGWAVE_FILTER_HAMMING,
+        GGWAVE_FILTER_FIRST_ORDER_HIGH_PASS,
+    } ggwave_Filter;
+
     // Operating modes of ggwave
     //
     //   GGWAVE_OPERATING_MODE_RX:
@@ -770,14 +776,41 @@ public:
     //
     //   src - input real-valued data, size is N
     //   dst - output complex-valued data, size is 2*N
-    //   ip  - work buffer, with size 2*N
-    //   w   - work buffer, with size 3 + sqrt(N/2)
+    //   wi  - work buffer, with size 2*N
+    //   wf  - work buffer, with size 3 + sqrt(N/2)
     //
-    //   First time calling thid function, make sure that ip[0] == 0
-    //   This will initialize some internal coefficients and store them in ip and w for
+    //   First time calling this function, make sure that wi[0] == 0
+    //   This will initialize some internal coefficients and store them in wi and wf for
     //   future usage.
     //
-    static bool computeFFTR(const float * src, float * dst, int N, int * ip, float * w);
+    //   If wi == nullptr                   - returns the needed size for wi
+    //   If wi != nullptr and wf == nullptr - returns the needed size for wf
+    //   If wi != nullptr and wf != nullptr - returns 1 on success, 0 on failure
+    //
+    static int computeFFTR(const float * src, float * dst, int N, int * wi, float * wf);
+
+    // Filter the waveform
+    //
+    //   filter   - filter to use
+    //   waveform - input waveform, size is N
+    //   N        - number of samples in the waveform
+    //   p0       - parameter
+    //   p1       - parameter
+    //   w        - work buffer
+    //
+    //   Filter is applied in-place.
+    //   First time calling this function, make sure that w[0] == 0 and w[1] == 0
+    //   This will initialize some internal coefficients and store them in w for
+    //   future usage.
+    //
+    //   For GGWAVE_FILTER_FIRST_ORDER_HIGH_PASS:
+    //     - p0 = cutoff frequency in Hz
+    //     - p1 = sample rate in Hz
+    //
+    //   If w == nullptr - returns the needed size for w for the specified filter
+    //   If w != nullptr - returns 1 on success, 0 on failure
+    //
+    static int filter(ggwave_Filter filter, float * waveform, int N, float p0, float p1, float * w);
 
     // Resample audio waveforms from one sample rate to another using sinc interpolation
     class Resampler {

src/ggwave.cpp

@@ -224,14 +224,14 @@ uint8_t getDSSMagic(int i) {
 #endif
 }
 
-void FFT(float * f, int N, int * ip, float * w) {
-    rdft(N, 1, f, ip, w);
+void FFT(float * f, int N, int * wi, float * wf) {
+    rdft(N, 1, f, wi, wf);
 }
 
-void FFT(const float * src, float * dst, int N, int * ip, float * w) {
+void FFT(const float * src, float * dst, int N, int * wi, float * wf) {
     memcpy(dst, src, N * sizeof(float));
 
-    FFT(dst, N, ip, w);
+    FFT(dst, N, wi, wf);
 }
 
 inline void addAmplitudeSmooth(
@@ -1300,10 +1300,82 @@ bool GGWave::computeFFTR(const float * src, float * dst, int N) {
     return true;
 }
 
-bool GGWave::computeFFTR(const float * src, float * dst, int N, int * ip, float * w) {
-    FFT(src, dst, N, ip, w);
+int GGWave::computeFFTR(const float * src, float * dst, int N, int * wi, float * wf) {
+    if (wi == nullptr) return 2*N;
+    if (wf == nullptr) return 3 + sqrt(N/2);
 
-    return true;
+    FFT(src, dst, N, wi, wf);
+
+    return 1;
+}
+
+int GGWave::filter(ggwave_Filter filter, float * waveform, int N, float p0, float p1, float * w) {
+    if (w == nullptr) {
+        switch (filter) {
+            case GGWAVE_FILTER_HANN:                  return N;
+            case GGWAVE_FILTER_HAMMING:               return N;
+            case GGWAVE_FILTER_FIRST_ORDER_HIGH_PASS: return 11;
+        };
+    }
+
+    if (w[0] == 0.0f && w[1] == 0.0f) {
+        switch (filter) {
+            case GGWAVE_FILTER_HANN:
+                {
+                    const float f = 2.0f*M_PI/(float)N;
+                    for (int i = 0; i < N; i++) {
+                        w[i] = 0.5f - 0.5f*cosf(f*(float)i);
+                    }
+                } break;
+            case GGWAVE_FILTER_HAMMING:
+                {
+                    const float f = 2.0f*M_PI/(float)N;
+                    for (int i = 0; i < N; i++) {
+                        w[i] = 0.54f - 0.46f*cosf(f*(float)i);
+                    }
+                } break;
+            case GGWAVE_FILTER_FIRST_ORDER_HIGH_PASS:
+                {
+                    const float th = 2.0f*M_PI*p0/p1;
+                    const float g = cos(th)/(1.0f + sin(th));
+                    w[0] = (1.0f + g)/2.0f;
+                    w[1] = -((1.0f + g)/2.0f);
+                    w[2] = 0.0f;
+                    w[3] = -g;
+                    w[4] = 0.0f;
+
+                    w[5] = 0.0f;
+                    w[6] = 0.0f;
+                    w[7] = 0.0f;
+                    w[8] = 0.0f;
+                } break;
+        };
+    }
+
+    switch (filter) {
+        case GGWAVE_FILTER_HANN:
+        case GGWAVE_FILTER_HAMMING:
+            {
+                for (int i = 0; i < N; i++) {
+                    waveform[i] *= w[i];
+                }
+            } break;
+        case GGWAVE_FILTER_FIRST_ORDER_HIGH_PASS:
+            {
+                for (int i = 0; i < N; i++) {
+                    float xn = waveform[i];
+                    float yn = w[0]*xn + w[1]*w[5] + w[2]*w[6] + w[3]*w[7] + w[4]*w[8];
+                    w[6] = w[5];
+                    w[5] = xn;
+                    w[8] = w[7];
+                    w[7] = yn;
+
+                    waveform[i] = yn;
+                }
+            } break;
+    };
+
+    return 1;
 }
 
 //