From 0ab93dfad3bdc5fe63e7ea6f1de41e6a1ec96e72 Mon Sep 17 00:00:00 2001 From: Georgi Gerganov Date: Tue, 12 Jul 2022 22:31:15 +0300 Subject: [PATCH] ggwave : add filter function Currently support the following filter: - Hann window - Hamming window - First-order high-pass --- include/ggwave/ggwave.h | 43 ++++++++++++++++++--- src/ggwave.cpp | 86 +++++++++++++++++++++++++++++++++++++---- 2 files changed, 117 insertions(+), 12 deletions(-) diff --git a/include/ggwave/ggwave.h b/include/ggwave/ggwave.h index 7690cc6..95a4acb 100644 --- a/include/ggwave/ggwave.h +++ b/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 { diff --git a/src/ggwave.cpp b/src/ggwave.cpp index 3abed39..3ca8eb1 100644 --- a/src/ggwave.cpp +++ b/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; } //