diff --git a/CMakeLists.txt b/CMakeLists.txt index fb985b4..6ac6335 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -116,7 +116,7 @@ endif() ######################################################################## # select the release build type by default to get optimization flags ######################################################################## -if(NOT CMAKE_BUILD_TYPE) +if (NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE "Release") message(STATUS "Build type not specified: defaulting to release.") endif(NOT CMAKE_BUILD_TYPE) @@ -196,7 +196,13 @@ endif() ###################################################### if (USE_TYPE_FLOAT) - add_library(PFDSP STATIC pf_mixer.cpp pf_mixer.h pf_cplx.h pf_carrier.cpp pf_carrier.h pf_cic.cpp pf_cic.h fmv.h ) + add_library(PFDSP STATIC + pf_windowing.cpp pf_windowing.h + pf_mixer.cpp pf_mixer.h + pf_carrier.cpp pf_carrier.h + pf_cic.cpp pf_cic.h + pf_cplx.h fmv.h + ) set_property(TARGET PFDSP PROPERTY CXX_STANDARD 11) set_property(TARGET PFDSP PROPERTY CXX_STANDARD_REQUIRED ON) set_target_properties(PFDSP PROPERTIES OUTPUT_NAME "pfdsp") @@ -216,7 +222,7 @@ if (USE_TYPE_FLOAT) ) if (INSTALL_PFDSP) set(INSTALL_TARGETS ${INSTALL_TARGETS} PFDSP) - set(INSTALL_HEADERS ${INSTALL_HEADERS} pf_mixer.h pf_cplx.h pf_carrier.h pf_cic.h) + set(INSTALL_HEADERS ${INSTALL_HEADERS} pf_windowing.h pf_mixer.h pf_carrier.h pf_cic.h pf_cplx.h) endif() endif() @@ -593,6 +599,18 @@ if (USE_TYPE_FLOAT) endif() target_link_libraries( bench_pf_conv_float pf_conv_dispatcher PFDSP $<$:stdc++> ) + ############################################################################ + + add_executable(bench_pf_window_float bench_windowing.cpp pf_helper_uclock.c pf_helper_uclock.h) + target_compile_definitions(bench_pf_window_float PRIVATE _USE_MATH_DEFINES) + target_compile_definitions(bench_pf_window_float PRIVATE PFFFT_ENABLE_FLOAT) + target_link_libraries( bench_pf_window_float PFDSP ${STDCXXLIB} ${ASANLIB} ) + + add_executable(test_windowing_float test_windowing.cpp pf_helper_uclock.c pf_helper_uclock.h) + target_compile_definitions(test_windowing_float PRIVATE _USE_MATH_DEFINES) + target_compile_definitions(test_windowing_float PRIVATE PFFFT_ENABLE_FLOAT) + target_link_libraries( test_windowing_float PFDSP ${STDCXXLIB} ${ASANLIB} ) + endif() ###################################################### @@ -659,5 +677,10 @@ if (USE_TYPE_FLOAT) WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} ) + add_test(NAME test_windowing_float + COMMAND "${CMAKE_CURRENT_BINARY_DIR}/test_windowing_float" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + ) + endif() diff --git a/bench_windowing.cpp b/bench_windowing.cpp new file mode 100644 index 0000000..2b34b97 --- /dev/null +++ b/bench_windowing.cpp @@ -0,0 +1,302 @@ +/* + Copyright (c) 2021 Hayati Ayguen ( h_ayguen@web.de ) + + bench for windowing algorithm/implementations + + */ + +#include +#include + +#include +#include +#include +#include + + +template +double bench_windowing_trig_cached(int N, int normalize, pf_windowing_T w) { + double t0, t1, tstop, T, nI, windSumI; + int iter, off, k; + int M = 512 * 1024 * 1024 / ( N * sizeof(TYPE) ); + float windScale; + TYPE *data = (TYPE *)malloc(M * N * sizeof(TYPE)); + float *wind = (float *)malloc(N * sizeof(TYPE)); + + for (k =0; k < M * N; ++k) + data[k] = 1.0F; + + iter = 0; + t0 = pf_uclock_sec(1); + tstop = t0 + 0.5; /* benchmark duration: 500 ms */ + + for (k =0; k < N; ++k) + wind[k] = 1.0F; + + /* calculate window */ + pf_win_trig(wind, w, N, 1.0F, 1, 0.16F); /* pf_win_trig(.., int N, float scale, int w0, float alpha) */ + + windScale = 1.0F; + if (normalize) + { + /* window sum - for normalization */ + windSumI = 0.0; + for (k =0; k < N; ++k) + windSumI += wind[k]; + windScale = (float)N / windSumI; + /* normalize window */ + for (k =0; k < N; ++k) + wind[k] *= windScale; + } + + int datapos = 0; + do { + TYPE *dataptr = data + datapos; + /* work: multiply with cached and scaled window wind[] */ + for (k =0; k < N; ++k) + dataptr[k] *= wind[k]; + datapos += N; + if (datapos >= N*M) + datapos = 0; + ++iter; + t1 = pf_uclock_sec(0); + } while ( t1 < tstop ); + + printf("windowing scale is %f\n", windScale); + T = ( t1 - t0 ); /* duration per fft() */ + printf("processed %d windows in %f ms\n", iter, T*1E3); + nI = iter / (T*1E3); + printf(" %f windows per ms\n", nI); + + free(wind); + free(data); + return nI; +} + + +template +double bench_windowing_trig_uncached(int N, int normalize, pf_windowing_T w) { + double t0, t1, tstop, T, nI; + SUMTYPE windSumI; + int iter, off, k; + int M = 512 * 1024 * 1024 / ( N * sizeof(TYPE) ); + float windScale; + TYPE *data = (TYPE *)malloc(M * N * sizeof(TYPE)); + for (k =0; k < M * N; ++k) + data[k] = 1.0F; + + iter = 0; + t0 = pf_uclock_sec(1); + tstop = t0 + 0.5; /* benchmark duration: 500 ms */ + + /* calculate window */ + pf_win_trig(data, w, N, 1.0F, 1, 0.16F); /* pf_win_trig(.., int N, float scale, int w0, float alpha) */ + /* window sum - for normalization */ + windSumI = 0.0; + for (k =0; k < N; ++k) + windSumI += data[k]; + windScale = normalize ? ((float)N / real(windSumI)) : 1.0F; + + int datapos = 0; + do { + TYPE *dataptr = data + datapos; + /* work */ + pf_win_trig(dataptr, w, N, windScale, 1, 0.16F); + datapos += N; + if (datapos >= N*M) + datapos = 0; + ++iter; + t1 = pf_uclock_sec(0); + } while ( t1 < tstop ); + + printf("windowing scale is %f\n", windScale); + T = ( t1 - t0 ); /* duration per fft() */ + printf("processed %d windows in %f ms\n", iter, T*1E3); + nI = iter / (T*1E3); + printf(" %f windows per ms\n", nI); + + free(data); + return nI; +} + + +template +double bench_windowing_uncached(int N, int normalize, pf_windowing_T w) { + double t0, t1, tstop, T, nI; + SUMTYPE windSumI; + int iter, off, k; + int M = 512 * 1024 * 1024 / ( N * sizeof(TYPE) ); + float windScale; + TYPE *data = (TYPE *)malloc(M * N * sizeof(TYPE)); + for (k =0; k < M * N; ++k) + data[k] = 1.0F; + + iter = 0; + t0 = pf_uclock_sec(1); + tstop = t0 + 0.5; /* benchmark duration: 500 ms */ + + /* calculate window */ + const struct pf_windowing_param_tag *wp = pf_window_alloc(w, N, 1, 0.16F); + /* window sum - for normalization */ + windSumI = pf_window_sum(wp); + windScale = normalize ? ((float)N / real(windSumI)) : 1.0F; + + int datapos = 0; + do { + TYPE *dataptr = data + datapos; + /* work */ + pf_win(dataptr, wp, windScale); + datapos += N; + if (datapos >= N*M) + datapos = 0; + ++iter; + t1 = pf_uclock_sec(0); + } while ( t1 < tstop ); + + printf("windowing scale is %f\n", windScale); + T = ( t1 - t0 ); /* duration per fft() */ + printf("processed %d windows in %f ms\n", iter, T*1E3); + nI = iter / (T*1E3); + printf(" %f windows per ms\n", nI); + + pf_window_free(wp); + free(data); + return nI; +} + + +template +double bench_windowing_cached(int N, int normalize, pf_windowing_T w) { + double t0, t1, tstop, T, nI; + SUMTYPE windSumI; + int iter, off, k; + int M = 512 * 1024 * 1024 / ( N * sizeof(TYPE) ); + float windScale; + TYPE *data = (TYPE *)malloc(M * N * sizeof(TYPE)); + float *wind = (float *)malloc(N * sizeof(TYPE)); + for (k =0; k < M * N; ++k) + data[k] = 1.0F; + + iter = 0; + t0 = pf_uclock_sec(1); + tstop = t0 + 0.5; /* benchmark duration: 500 ms */ + + for (k =0; k < N; ++k) + wind[k] = 1.0F; + + /* calculate window */ + const struct pf_windowing_param_tag *wp = pf_window_alloc(w, N, 1, 0.16F); +#if 0 + /* window sum - for normalization */ + windSumI = pf_window_sum(wp); + windScale = normalize ? ((float)N / real(windSumI)) : 1.0F; + pf_window(wind, wp, windScale); +#else + pf_window(wind, wp); + + windScale = 1.0F; + if (normalize) + { + /* window sum - for normalization */ + windSumI = 0.0; + for (k =0; k < N; ++k) + windSumI += wind[k]; + windScale = (float)N / windSumI; + /* normalize window */ + for (k =0; k < N; ++k) + wind[k] *= windScale; + } +#endif + + int datapos = 0; + do { + TYPE *dataptr = data + datapos; + /* work: multiply with cached and scaled window wind[] */ + for (k =0; k < N; ++k) + dataptr[k] *= wind[k]; + datapos += N; + if (datapos >= N*M) + datapos = 0; + ++iter; + t1 = pf_uclock_sec(0); + } while ( t1 < tstop ); + + printf("windowing scale is %f\n", windScale); + T = ( t1 - t0 ); /* duration per fft() */ + printf("processed %d windows in %f ms\n", iter, T*1E3); + nI = iter / (T*1E3); + printf(" %f windows per ms\n", nI); + + pf_window_free(wp); + free(wind); + free(data); + return nI; +} + + +template +void bench(const char * type_str, const char * sum_str, const int N, const int normalize) { + const pf_windowing_T aBenchWins[] = { pf_winHamming, pf_winBlackmanHarris, pf_winFlatTop_HFT90D }; + const char * acWinStrs[] = { "pf_winHamming", "pf_winBlackmanHarris", "pf_winFlatTop_HFT90D" }; + double rt, rtA, rtB, rtC, rtD; + + for ( int wno = 0; wno < 3; ++wno ) + { + printf("\n==============================================================\n"); + printf("\nstarting bench of bench_windowing_trig_uncached<%s, %s>() with %s\n", type_str, "complexd", acWinStrs[wno]); + rtA = bench_windowing_trig_uncached(N, normalize, aBenchWins[wno]); + /* printf(" %f MSamples/sec\n\n", rt * 1E-6); */ + + printf("\nstarting bench of bench_windowing_trig_cached<%s>() with %s\n", type_str, acWinStrs[wno]); + rtB = bench_windowing_trig_cached(N, normalize, aBenchWins[wno]); + + printf("\nstarting bench of bench_windowing_uncached<%s, %s>() with %s\n", type_str, sum_str, acWinStrs[wno]); + rtC = bench_windowing_uncached(N, normalize, aBenchWins[wno]); + + printf("\nstarting bench of bench_windowing_cached<%s, %s>() with %s\n", type_str, sum_str, acWinStrs[wno]); + rtD = bench_windowing_cached(N, normalize, aBenchWins[wno]); + + printf("\nwindowing_trig_cached() / windowing() = %f\n", rtB / rtC); + printf("windowing() / windowing_trig_uncached() = %f\n", rtC / rtA); + printf("bench_windowing_cached() / windowing_trig_cached() = %f\n", rtD / rtB); + } + printf("\n==============================================================\n"); +} + + +int main(int argc, char **argv) +{ + int N = 1024; + int showUsage = 0; + int normalize = 1; + int cx = 1; + + if (argc == 1) + showUsage = 1; + if (1 < argc) + N = atoi(argv[1]); + if (2 < argc) + normalize = atoi(argv[2]); + if (3 < argc) + cx = atoi(argv[3]); + if ( !N || showUsage ) + { + fprintf(stderr, "%s [ [ [] ] ]\n", argv[0]); + return 0; + } + + fprintf(stderr, "will benchmark window size N = %d with%s for %s data\n", N, normalize?"":"out", cx ? "complex":"real"); + pf_uclock_sec(0); + + if (!cx) + { + bench("float", "double", N, normalize); + } + else + { + bench("complexf", "double", N, normalize); + } + + return 0; +} + diff --git a/pf_cplx.h b/pf_cplx.h index 61d8486..7f16fe1 100644 --- a/pf_cplx.h +++ b/pf_cplx.h @@ -29,6 +29,11 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #pragma once +#ifdef __cplusplus +extern "C" { +#endif + + /* _____ _ / ____| | | @@ -40,5 +45,99 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |_| */ -typedef struct complexf_s { float i; float q; } complexf; +typedef struct complexf_s { + float i; + float q; + +#ifdef __cplusplus + inline complexf_s & operator=(float r) { + this->i = r; + this->q = 0.0F; + return *this; + } +#endif + +} complexf; + + +typedef struct complexd_s { + double i; + double q; + +#ifdef __cplusplus + inline complexd_s & operator=(double r) { + this->i = r; + this->q = 0.0; + return *this; + } +#endif + +} complexd; + + + +#ifdef __cplusplus +} + +/* add/define some helper functions */ + +static inline float real(const float self) { + return self; +} + +static inline float real(const complexf &self) { + return self.i; +} + +static inline float imag(const complexf &self) { + return self.q; +} + +static inline void operator+=(complexf &self, float r) { + self.i += r; +} + +static inline void operator+=(complexf &self, const complexf sum) { + self.i += sum.i; + self.q += sum.q; +} + +static inline void operator*=(complexf &self, float r) { + self.i *= r; + self.q *= r; +} + + + +static inline double real(const double self) { + return self; +} + +static inline double real(const complexd &self) { + return self.i; +} + +static inline double imag(const complexd &self) { + return self.q; +} + +static inline void operator+=(complexd &self, double r) { + self.i += r; +} + +static inline void operator+=(complexd &self, const complexd sum) { + self.i += sum.i; + self.q += sum.q; +} + +static inline void operator+=(complexd &self, const complexf sum) { + self.i += sum.i; + self.q += sum.q; +} + +static inline void operator*=(complexd &self, double r) { + self.i *= r; + self.q *= r; +} +#endif diff --git a/pf_helper_uclock.c b/pf_helper_uclock.c new file mode 100644 index 0000000..33f769a --- /dev/null +++ b/pf_helper_uclock.c @@ -0,0 +1,76 @@ +/* + Copyright (c) 2021 Hayati Ayguen ( h_ayguen@web.de ) + + helper function(s) + + */ + +#include "pf_helper_uclock.h" + +/* #include */ +#include +#include + +#if defined(__linux__) +#define HAVE_SYS_TIMES +#endif + +#ifdef HAVE_SYS_TIMES +# include +# include +#endif + +#ifdef WIN32 +#define WIN32_LEAN_AND_MEAN +#define VC_EXTRALEAN +#include +#endif + + +#if defined(HAVE_SYS_TIMES) + static double ttclk = 0.; + + double pf_uclock_sec(int find_start) + { + struct tms t0, t; + if (ttclk == 0.) + { + ttclk = sysconf(_SC_CLK_TCK); + fprintf(stderr, "sysconf(_SC_CLK_TCK) => %f\n", ttclk); + } + times(&t); + if (find_start) + { + t0 = t; + while (t0.tms_utime == t.tms_utime) + times(&t); + } + /* use only the user time of this process - not realtime, which depends on OS-scheduler .. */ + return ((double)t.tms_utime) / ttclk; + } + +#elif defined(WIN32) + // https://docs.microsoft.com/en-us/windows/win32/api/processthreadsapi/nf-processthreadsapi-getprocesstimes + double pf_uclock_sec(int find_start) + { + FILETIME a, b, c, d; + if (GetProcessTimes(GetCurrentProcess(), &a, &b, &c, &d) != 0) + { + // Returns total user time. + // Can be tweaked to include kernel times as well. + return + (double)(d.dwLowDateTime | + ((unsigned long long)d.dwHighDateTime << 32)) * 0.0000001; + } + else { + // Handle error + return 0; + } + } + +#else + double pf_uclock_sec(int find_start) + { return (double)clock()/(double)CLOCKS_PER_SEC; } +#endif + + diff --git a/pf_helper_uclock.h b/pf_helper_uclock.h new file mode 100644 index 0000000..a5d81c7 --- /dev/null +++ b/pf_helper_uclock.h @@ -0,0 +1,41 @@ +/* +This software is part of pffft/pfdsp, a set of simple DSP routines. + +Copyright (c) 2021 Hayati Ayguen +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the copyright holder nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL ANDRAS RETZLER BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#pragma once + +#ifdef __cplusplus +extern "C" { +#endif + +double pf_uclock_sec(int find_start); + +#ifdef __cplusplus +} +#endif + diff --git a/pf_mixer.h b/pf_mixer.h index e153ad0..f8b3f96 100644 --- a/pf_mixer.h +++ b/pf_mixer.h @@ -30,11 +30,11 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #pragma once +#include "pf_cplx.h" + #include #include -#include "pf_cplx.h" - #ifdef __cplusplus extern "C" { #endif diff --git a/pf_windowing.cpp b/pf_windowing.cpp new file mode 100644 index 0000000..bb6c6a4 --- /dev/null +++ b/pf_windowing.cpp @@ -0,0 +1,1212 @@ +/* +This software is part of pffft/pfdsp, a set of simple DSP routines. + +Copyright (c) 2021 Hayati Ayguen +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the copyright holder nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL ANDRAS RETZLER BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#include "pf_windowing.h" +#include "fmv.h" + +#include +#include +#include + + +typedef float freal; +typedef double frealsum; + +typedef freal freal2[2]; + +#define rcos(X) cosf(X) +#define rsin(X) sinf(X) +#define rsqrt(X) sqrtf(X) + +#define TRIG_TAB_SIZE 4 + +//they dropped M_PI in C99, so we define it: +#define PI ((freal)3.14159265358979323846) +#define TWO_PI ((freal)(2.0 * 3.14159265358979323846)) + +#define c0 0.0F +#define c1 1.0F +#define c2 2.0F +#define c3 3.0F +#define c4 4.0F + +#define c_0p5 0.5F + + +typedef struct pf_windowing_param_tag +{ + freal trig_tab[TRIG_TAB_SIZE +1][2]; + freal a[5]; // 0 .. 4 + freal N_div; + int N_div_i; // N_div as integer + int num_a_coeff; + int N; + int w0; +} pf_windowing_param_t; + + +const char * pf_window_text(pf_windowing_T win) +{ + switch (win) + { + case pf_winRect: return "Rectangle / Boxcar"; + case pf_winHann: return "von Hann"; + case pf_winHamming: return "Hamming"; + case pf_winBlackman: return "Blackman"; + case pf_winBlackmanExact: return "Blackman (Exact)"; + case pf_winNutall: return "Nutall"; + case pf_winBlackmanHarris: return "Blackman-Harris"; + case pf_winBlackmanNutall: return "Blackman-Nutall"; + case pf_winFlatTop_MATLAB: return "FlatTop (Matlab)"; + case pf_winFlatTop_SR785: return "FlatTop (SR785)"; + case pf_winFlatTop_HFT70: return "FlatTop (HFT70)"; + case pf_winFlatTop_HFT95: return "FlatTop (HFT95)"; + case pf_winFlatTop_HFT90D: return "FlatTop (HFT90D)"; + case pf_num_windows: + default: ; + } + return 0; +} + + +static +pf_windowing_param_t pf_window_param(int calc_cordic, pf_windowing_T win, int N, int w0, freal alpha) +{ + pf_windowing_param_t wp; + //pf_windowing_param_t *p = ℘ + // classical alpha for 3-term Blackman: 0.16 + const freal bl_alpha = (alpha >= 0) ? alpha : (freal)(0.16); + wp.a[0] = (freal)(1.0); + wp.a[1] = (freal)(0.0); + wp.a[2] = (freal)(0.0); + wp.a[3] = (freal)(0.0); + wp.a[4] = (freal)(0.0); + wp.N_div = (freal)(N -1); + wp.N_div_i = N - 1; + wp.num_a_coeff = 1; + wp.N = N; + wp.w0 = w0; + + switch( win ) + { + case pf_winRect: + break; + default: assert(0); + case pf_winHann: + wp.a[0] = (freal)(0.5); + wp.a[1] = (freal)(-0.5); + wp.num_a_coeff = 2; + break; + case pf_winHamming: + wp.a[0] = (freal)( 25.0 / 46.0); // ~= 0.54 + wp.a[1] = (freal)(-21.0 / 46.0); // ~= 0.46 + wp.num_a_coeff = 2; + break; + case pf_winBlackman: + wp.a[0] = c_0p5 * (c1 - bl_alpha); // default: 0.42 + wp.a[1] = (freal)(-0.5); + wp.a[2] = c_0p5 * bl_alpha; // default: 0.08 + wp.num_a_coeff = 3; + break; + case pf_winBlackmanExact: + wp.a[0] = (freal)( 7938.0/18608.0); // = 0.42659 + wp.a[1] = (freal)(-9240.0/18608.0); // = 0.49656 + wp.a[2] = (freal)( 1430.0/18608.0); // = 0.076849 + wp.num_a_coeff = 3; + break; + case pf_winNutall: + wp.a[0] = (freal)( 0.355768); + wp.a[1] = (freal)(-0.487396); + wp.a[2] = (freal)( 0.144232); + wp.a[3] = (freal)(-0.012604); + wp.num_a_coeff = 4; + break; + case pf_winBlackmanHarris: + wp.a[0] = (freal)( 0.35875); + wp.a[1] = (freal)(-0.48829); + wp.a[2] = (freal)( 0.14128); + wp.a[3] = (freal)(-0.01168); + wp.num_a_coeff = 4; + break; + case pf_winBlackmanNutall: + wp.a[0] = (freal)( 0.3635819); + wp.a[1] = (freal)(-0.4891775); + wp.a[2] = (freal)( 0.1365995); + wp.a[3] = (freal)(-0.0106411); + wp.num_a_coeff = 4; + break; + case pf_winFlatTop_MATLAB: + wp.a[0] = (freal)( 0.21557895); + wp.a[1] = (freal)(-0.41663158); + wp.a[2] = (freal)( 0.277263158); + wp.a[3] = (freal)(-0.083578947); + wp.a[4] = (freal)( 0.006947368); + //wp.N_div = (freal)(N -1); + //wp.N_div_i = N -1; + wp.num_a_coeff = 5; + break; + case pf_winFlatTop_SR785: + wp.a[0] = (freal)( 1.0); + wp.a[1] = (freal)(-1.93); + wp.a[2] = (freal)( 1.29); + wp.a[3] = (freal)(-0.388); + wp.a[4] = (freal)( 0.028); + wp.num_a_coeff = 5; + break; + case pf_winFlatTop_HFT70: + wp.a[0] = (freal)( 1.0); + wp.a[1] = (freal)(-1.90796); + wp.a[2] = (freal)( 1.07349); + wp.a[3] = (freal)(-0.18199); + wp.num_a_coeff = 4; + break; + case pf_winFlatTop_HFT95: + wp.a[0] = (freal)( 1.0); + wp.a[1] = (freal)(-1.9383379); + wp.a[2] = (freal)( 1.3045202); + wp.a[3] = (freal)(-0.4028270); + wp.a[4] = (freal)( 0.0350665); + wp.num_a_coeff = 5; + break; + case pf_winFlatTop_HFT90D: + wp.a[0] = (freal)( 1.0); + wp.a[1] = (freal)(-1.942604); + wp.a[2] = (freal)( 1.340318); + wp.a[3] = (freal)(-0.440811); + wp.a[4] = (freal)( 0.043097); + wp.num_a_coeff = 5; + break; + } + + + if (!calc_cordic) + return wp; + // pre-calculate small trig table - hope this allow SIMD parallelization + freal2 *trig_tab = wp.trig_tab; + const freal phi_inc = TWO_PI / wp.N_div; + const freal c_inc_cos = rcos(phi_inc); + const freal c_inc_sin = rsin(phi_inc); + freal cos_phi = trig_tab[0][0] = c1; + freal sin_phi = trig_tab[0][1] = c0; +#if 0 + printf("pf_window_param(N = %d, %s)\n", N, pf_window_text(win)); +#endif + for ( int k = 1; k <= TRIG_TAB_SIZE; ++k ) + { + // rotate + const freal cos_sav = cos_phi; + cos_phi = c_inc_cos * cos_sav - c_inc_sin * sin_phi; + sin_phi = c_inc_sin * cos_sav + c_inc_cos * sin_phi; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi * cos_phi + sin_phi * sin_phi); + cos_phi /= mag; + sin_phi /= mag; + trig_tab[k][0] = cos_phi; + trig_tab[k][1] = sin_phi; +#if 0 + printf("pf_window_param() trig_tab[%d] = %.3g + i * %.3g. %.3g + i %.3g. dphi = %g\n", k + , trig_tab[k][0], trig_tab[k][1] + , rcos(k*phi_inc), rsin(k*phi_inc) + , phi_inc + ); +#endif + } + + return wp; +} + + +const struct pf_windowing_param_tag * pf_window_alloc(pf_windowing_T win, int N, int w0, freal alpha) +{ + struct pf_windowing_param_tag * param = new pf_windowing_param_t(); + *param = pf_window_param(1, win, N, w0, alpha); + return param; +} + + +void pf_get_window_param(const struct pf_windowing_param_tag * param, int *N_div, int *num_alpha_coeff) +{ + if (N_div) + *N_div = param->N_div_i; + if (num_alpha_coeff) + *num_alpha_coeff = param->num_a_coeff; +} + + +void pf_window_free(const struct pf_windowing_param_tag * param) +{ + delete param; +} + + + +PF_TARGET_CLONES +void pf_window_trig(freal * pafWin, pf_windowing_T win, int N, freal scale, int w0, freal alpha) +{ + const freal phi0 = w0 ? c0 : PI; + const pf_windowing_param_t wp = pf_window_param(0, win, N, w0, alpha); + + switch( wp.num_a_coeff ) + { + case 1: + { + if (scale == (freal)(1.0)) + break; + const freal a0 = scale; + for ( int n = 0; n < N; ++n ) + pafWin[n] *= a0; + } + break; + default: assert(0); + case 2: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + for ( int n = 0; n < N; ++n ) + { + const freal phi = phi0 + ( TWO_PI * n ) / wp.N_div; + const freal wn = a0 + a1 * rcos(phi); + pafWin[n] *= scale * wn; + } + } + break; + case 3: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + for ( int n = 0; n < N; ++n ) + { + const freal phi = phi0 + ( TWO_PI * n ) / wp.N_div; + const freal wn = a0 + a1 * rcos(phi) + + a2 * rcos(c2 * phi); + pafWin[n] *= scale * wn; + } + } + break; + case 4: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + for ( int n = 0; n < N; ++n ) + { + const freal phi = phi0 + ( TWO_PI * n ) / wp.N_div; + const freal wn = a0 + a1 * rcos(phi) + + a2 * rcos(c2 * phi) + + a3 * rcos(c3 * phi); + pafWin[n] *= scale * wn; + } + } + break; + case 5: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + const freal a4 = wp.a[4]; + for ( int n = 0; n < N; ++n ) + { + const freal phi = phi0 + ( TWO_PI * n ) / wp.N_div; + const freal wn = a0 + a1 * rcos(phi) + + a2 * rcos(c2 * phi) + + a3 * rcos(c3 * phi) + + a4 * rcos(c4 * phi); + pafWin[n] *= scale * wn; + } + } + break; + } +} + + +PF_TARGET_CLONES +void pf_window_trig_cx(complexf * pacWin, pf_windowing_T win, int N, freal scale, int w0, freal alpha) +{ + const freal phi0 = w0 ? c0 : PI; + const pf_windowing_param_t wp = pf_window_param(0, win, N, w0, alpha); + + switch( wp.num_a_coeff ) + { + case 1: + { + if (scale == (freal)(1.0)) + break; + const freal a0 = scale; + for ( int n = 0; n < N; ++n ) + { + pacWin[n].i *= a0; + pacWin[n].q *= a0; + } + } + break; + default: assert(0); + case 2: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + for ( int n = 0; n < N; ++n ) + { + const freal phi = phi0 + ( TWO_PI * n ) / wp.N_div; + const freal wn = a0 + a1 * rcos(phi); + pacWin[n].i *= scale * wn; + pacWin[n].q *= scale * wn; + } + } + break; + case 3: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + for ( int n = 0; n < N; ++n ) + { + const freal phi = phi0 + ( TWO_PI * n ) / wp.N_div; + const freal wn = a0 + a1 * rcos(phi) + + a2 * rcos(c2 * phi); + pacWin[n].i *= scale * wn; + pacWin[n].q *= scale * wn; + } + } + break; + case 4: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + for ( int n = 0; n < N; ++n ) + { + const freal phi = phi0 + ( TWO_PI * n ) / wp.N_div; + const freal wn = a0 + a1 * rcos(phi) + + a2 * rcos(c2 * phi) + + a3 * rcos(c3 * phi); + pacWin[n].i *= scale * wn; + pacWin[n].q *= scale * wn; + } + } + break; + case 5: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + const freal a4 = wp.a[4]; + for ( int n = 0; n < N; ++n ) + { + const freal phi = phi0 + ( TWO_PI * n ) / wp.N_div; + const freal wn = a0 + a1 * rcos(phi) + + a2 * rcos(c2 * phi) + + a3 * rcos(c3 * phi) + + a4 * rcos(c4 * phi); + pacWin[n].i *= scale * wn; + pacWin[n].q *= scale * wn; + } + } + break; + } +} + + +frealsum pf_window_sum(const struct pf_windowing_param_tag * param) +{ + const pf_windowing_param_t &wp = *param; + const int N = wp.N; + const freal2 *trig_tab = wp.trig_tab; + // w0 != 0: start with phi = 0 + // w0 == 0: start with phi = ( TWO_PI * (-N/2) ) / N == -PI + freal cos_phi0 = wp.w0 ? c1 : -c1; + freal sin_phi0 = c0; + frealsum wsum = (frealsum)(0.0); + + switch( wp.num_a_coeff ) + { + case 1: + wsum += N; + break; + default: assert(0); + case 2: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + //const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal wn = a0 + a1 * cos_phi; + wsum += wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 3: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi; + wsum += wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 4: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi + a3 * cos_3phi; + wsum += wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 5: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + const freal a4 = wp.a[4]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal cos_4phi = cos_2phi * cos_2phi - sin_2phi * sin_2phi; + const freal wn = a0 +a1 * cos_phi +a2 * cos_2phi + +a3 * cos_3phi +a4 * cos_4phi; + wsum += wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + } + return wsum; +} + + +PF_TARGET_CLONES +void pf_window_no_scale(freal * pafWin, const struct pf_windowing_param_tag * param) +{ + const pf_windowing_param_t &wp = *param; + const int N = wp.N; + const freal2 *trig_tab = wp.trig_tab; + // w0 != 0: start with phi = 0 + // w0 == 0: start with phi = ( TWO_PI * (-N/2) ) / N == -PI + freal cos_phi0 = wp.w0 ? c1 : -c1; + freal sin_phi0 = c0; + + switch( wp.num_a_coeff ) + { + case 1: + break; + default: assert(0); + case 2: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + //const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal wn = a0 + a1 * cos_phi; + pafWin[n+k] *= wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 3: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi; + pafWin[n+k] *= wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 4: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi + a3 * cos_3phi; + pafWin[n+k] *= wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 5: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + const freal a4 = wp.a[4]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal cos_4phi = cos_2phi * cos_2phi - sin_2phi * sin_2phi; + const freal wn = a0 +a1 * cos_phi +a2 * cos_2phi + +a3 * cos_3phi +a4 * cos_4phi; + pafWin[n+k] *= wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + } +} + + +PF_TARGET_CLONES +void pf_window(freal * pafWin, const struct pf_windowing_param_tag * param, freal scale) +{ + if (scale == c1) + { + pf_window_no_scale(pafWin, param); + return; + } + + const pf_windowing_param_t &wp = *param; + const int N = wp.N; + const freal2 *trig_tab = wp.trig_tab; + // w0 != 0: start with phi = 0 + // w0 == 0: start with phi = ( TWO_PI * (-N/2) ) / N == -PI + freal cos_phi0 = wp.w0 ? c1 : -c1; + freal sin_phi0 = c0; + + switch( wp.num_a_coeff ) + { + case 1: + { + if (scale == (freal)(1.0)) + break; + const freal a0 = scale; + for ( int n =0; n < N; ++n ) + pafWin[n] *= a0; + } + break; + default: assert(0); + case 2: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + //const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal wn = a0 + a1 * cos_phi; + pafWin[n+k] *= scale * wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 3: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi; + pafWin[n+k] *= scale * wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 4: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi + a3 * cos_3phi; + pafWin[n+k] *= scale * wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 5: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + const freal a4 = wp.a[4]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal cos_4phi = cos_2phi * cos_2phi - sin_2phi * sin_2phi; + const freal wn = a0 +a1 * cos_phi +a2 * cos_2phi + +a3 * cos_3phi +a4 * cos_4phi; + pafWin[n+k] *= scale * wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + } +} + + +PF_TARGET_CLONES +void pf_window_cx_no_scale(complexf * pacWin, const struct pf_windowing_param_tag * param) +{ + const pf_windowing_param_t &wp = *param; + const int N = wp.N; + const freal2 *trig_tab = wp.trig_tab; + // w0 != 0: start with phi = 0 + // w0 == 0: start with phi = ( TWO_PI * (-N/2) ) / N == -PI + freal cos_phi0 = wp.w0 ? c1 : -c1; + freal sin_phi0 = c0; + + switch( wp.num_a_coeff ) + { + case 1: + break; + default: assert(0); + case 2: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + //const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal wn = a0 + a1 * cos_phi; + pacWin[n+k].i *= wn; + pacWin[n+k].q *= wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 3: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi; + pacWin[n+k].i *= wn; + pacWin[n+k].q *= wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 4: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi + a3 * cos_3phi; + pacWin[n+k].i *= wn; + pacWin[n+k].q *= wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 5: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + const freal a4 = wp.a[4]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal cos_4phi = cos_2phi * cos_2phi - sin_2phi * sin_2phi; + const freal wn = a0 +a1 * cos_phi +a2 * cos_2phi + +a3 * cos_3phi +a4 * cos_4phi; + pacWin[n+k].i *= wn; + pacWin[n+k].q *= wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + } +} + + +PF_TARGET_CLONES +void pf_window_cx(complexf * pacWin, const struct pf_windowing_param_tag * param, freal scale) +{ + if (scale == c1) + { + pf_window_cx_no_scale(pacWin, param); + return; + } + + const pf_windowing_param_t &wp = *param; + const int N = wp.N; + const freal2 *trig_tab = wp.trig_tab; + // w0 != 0: start with phi = 0 + // w0 == 0: start with phi = ( TWO_PI * (-N/2) ) / N == -PI + freal cos_phi0 = wp.w0 ? c1 : -c1; + freal sin_phi0 = c0; + + switch( wp.num_a_coeff ) + { + case 1: + { + if (scale == (freal)(1.0)) + break; + const freal a0 = scale; + for ( int n =0; n < N; ++n ) + { + pacWin[n].i *= a0; + pacWin[n].q *= a0; + } + } + break; + default: assert(0); + case 2: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + //const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal wn = a0 + a1 * cos_phi; + pacWin[n+k].i *= scale * wn; + pacWin[n+k].q *= scale * wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 3: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi; + pacWin[n+k].i *= scale * wn; + pacWin[n+k].q *= scale * wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 4: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal wn = a0 + a1 * cos_phi + a2 * cos_2phi + a3 * cos_3phi; + pacWin[n+k].i *= scale * wn; + pacWin[n+k].q *= scale * wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + case 5: + { + const freal a0 = wp.a[0]; + const freal a1 = wp.a[1]; + const freal a2 = wp.a[2]; + const freal a3 = wp.a[3]; + const freal a4 = wp.a[4]; + for ( int n =0; n < N; ) + { + const int M = (N - n >= TRIG_TAB_SIZE) ? TRIG_TAB_SIZE : (N - n); + for ( int k = 0; k < M; ++k ) + { + const freal cos_phi = trig_tab[k][0] * cos_phi0 - trig_tab[k][1] * sin_phi0; + const freal sin_phi = trig_tab[k][1] * cos_phi0 + trig_tab[k][0] * sin_phi0; + const freal cos_2phi = cos_phi * cos_phi - sin_phi * sin_phi; + const freal sin_2phi = c2 * sin_phi * cos_phi; + const freal cos_3phi = cos_2phi * cos_phi - sin_2phi * sin_phi; + const freal cos_4phi = cos_2phi * cos_2phi - sin_2phi * sin_2phi; + const freal wn = a0 +a1 * cos_phi +a2 * cos_2phi + +a3 * cos_3phi +a4 * cos_4phi; + pacWin[n+k].i *= scale * wn; + pacWin[n+k].q *= scale * wn; + } + n += TRIG_TAB_SIZE; + if ( n >= N ) + break; + // rotate + const freal cos_sav = cos_phi0; + cos_phi0 = trig_tab[TRIG_TAB_SIZE][0] * cos_sav - trig_tab[TRIG_TAB_SIZE][1] * sin_phi0; + sin_phi0 = trig_tab[TRIG_TAB_SIZE][1] * cos_sav + trig_tab[TRIG_TAB_SIZE][0] * sin_phi0; + // re-normalize cos_phi/sin_phi to unit-circle + const freal mag = rsqrt(cos_phi0 * cos_phi0 + sin_phi0 * sin_phi0); + cos_phi0 /= mag; + sin_phi0 /= mag; + } + } + break; + } +} + diff --git a/pf_windowing.h b/pf_windowing.h new file mode 100644 index 0000000..9a8850e --- /dev/null +++ b/pf_windowing.h @@ -0,0 +1,161 @@ +/* +This software is part of pffft/pfdsp, a set of simple DSP routines. + +Copyright (c) 2021 Hayati Ayguen +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the copyright holder nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL ANDRAS RETZLER BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#pragma once + +#include "pf_cplx.h" + +#ifdef __cplusplus +extern "C" { +#endif + +struct pf_windowing_param_tag; + +/* https://en.wikipedia.org/wiki/Window_function */ + +typedef enum +{ + pf_winRect =0 + , pf_winHann + , pf_winHamming + , pf_winBlackman + , pf_winBlackmanExact + , pf_winNutall + , pf_winBlackmanHarris + , pf_winBlackmanNutall + , pf_winFlatTop_MATLAB /* https://en.wikipedia.org/wiki/Window_function + * https://www.mathworks.com/help/signal/ref/flattopwin.html + */ + , pf_winFlatTop_SR785 /* Flat Top of spectrum-analyzer SR785 of Stanford Research Systems + * see https://de.wikipedia.org/wiki/Fensterfunktion + */ + , pf_winFlatTop_HFT70 /* "HFT70" from D.3.1 page 45 of + * https://pure.mpg.de/rest/items/item_152164_1/component/file_152163/content + * first zero at f = +/- 4.00 bins + * highest sidelobe: -70.4 dB at f = +/- 4.65 bins + * sidelobe drop at rate f^-1 + * optimal overlap: 72.2% + * amplitude flatness: 0.964 + * power flatness: 0.637 + * overlap correlation: 0.041 + */ + , pf_winFlatTop_HFT95 /* "HFT95" from D.3.2 page 46 of + * https://pure.mpg.de/rest/items/item_152164_1/component/file_152163/content + * first zero at f = +/- 5.00 bins + * highest sidelobe: -95.0 dB at f = +/- 7.49 bins + * sidelobe drop at rate f^-1 + * optimal overlap: 75.6% + * amplitude flatness: 0.952 + * power flatness: 0.647 + * overlap correlation: 0.056 + */ + , pf_winFlatTop_HFT90D /* "HFT90D" from D.3.3 page 46 of + * https://pure.mpg.de/rest/items/item_152164_1/component/file_152163/content + * first zero at f = +/- 5.00 bins + * highest sidelobe: -90.2 dB at f = +/- 5.58 bins + * sidelobe drop at rate f^-3 + * optimal overlap: 76.0% + * amplitude flatness: 0.953 + * power flatness: 0.646 + * overlap correlation: 0.054 + */ + , pf_num_windows +} pf_windowing_T; + + +const char * pf_window_text(pf_windowing_T win); + +void pf_window_free(const struct pf_windowing_param_tag * param); + +/* + * scale: for additional scaling, e.g. normalization of window or FFT scale compensation + * + * w0: != 0 => top at center of coefficients; towards 0 at corners + * == 0 => top at corners; towards 0 at center + * + * alpha: for blackman window; negative value uses default of 0.16 + * + */ + +#ifdef __cplusplus +// add default values for some parameters + +// _trig functions utilize the (slow) trigonometric functions +void pf_window_trig( float * pafWin, pf_windowing_T win, int N, float scale =1.0F, int w0 =1, float alpha =-1.0F); +void pf_window_trig_cx(complexf * pacWin, pf_windowing_T win, int N, float scale =1.0F, int w0 =1, float alpha =-1.0F); + +const struct pf_windowing_param_tag * pf_window_alloc(pf_windowing_T win, int N, int w0 = 1, float alpha =-1.0F); +void pf_get_window_param(const struct pf_windowing_param_tag * param, int *N_div =0, int *num_alpha_coeff =0); + +double pf_window_sum(const struct pf_windowing_param_tag * param); +void pf_window( float * pafWin, const struct pf_windowing_param_tag * param, float scale =1.0F); +void pf_window_cx( complexf * pacWin, const struct pf_windowing_param_tag * param, float scale =1.0F); + +#else +/* no overloads and no default parameters with pure C */ + +/* _trig functions utilize the (slow) trigonometric functions */ +void pf_window_trig( float * pafWin, pf_windowing_T win, int N, float scale, int w0, float alpha); +void pf_window_trig_cx(complexf * pacWin, pf_windowing_T win, int N, float scale, int w0, float alpha); + +const struct pf_windowing_param_tag * pf_window_alloc(pf_windowing_T win, int N, int w0, float alpha); +void pf_get_window_param(const struct pf_windowing_param_tag * param, int *N_div, int *num_alpha_coeff); + +double pf_window_sum(const struct pf_windowing_param_tag * param); +void pf_window( float * pafWin, const struct pf_windowing_param_tag * param, float scale); +void pf_window_cx( complexf * pacWin, const struct pf_windowing_param_tag * param, float scale); +#endif + + + +#ifdef __cplusplus +} + + +// additional C++ inline function pf_win_trig() with overloads: + +inline void pf_win_trig(float * pafWin, pf_windowing_T win, int N, float scale =1.0F, int w0 =1, float alpha =-1.0F) { + pf_window_trig(pafWin, win, N, scale, w0, alpha); +} + +inline void pf_win_trig(complexf * pacWin, pf_windowing_T win, int N, float scale =1.0F, int w0 =1, float alpha =-1.0F) { + pf_window_trig_cx(pacWin, win, N, scale, w0, alpha); +} + + +inline void pf_win(float * pafWin, const struct pf_windowing_param_tag * param, float scale = 1.0F ) { + pf_window(pafWin, param, scale); +} + +inline void pf_win(complexf * pacWin, const struct pf_windowing_param_tag * param, float scale = 1.0F ) { + pf_window_cx(pacWin, param, scale); +} + +#endif + diff --git a/test_windowing.cpp b/test_windowing.cpp new file mode 100644 index 0000000..04c00e5 --- /dev/null +++ b/test_windowing.cpp @@ -0,0 +1,100 @@ + +#include + +#include +#include + +#define WIN_SIZE 64 +#define ERR_LIMIT 1E-5 + +int main(int argc, char** argv) +{ + float win_coeff_trig[WIN_SIZE], win_coeff_cordic[WIN_SIZE], delta[WIN_SIZE]; + int showUsage = 0; + + if (argc == 1) + showUsage = 1; + + for (int winNo = 0; winNo < pf_num_windows; ++winNo) + { + for (int w0 = 0; w0 <= 1; ++w0 ) + { + pf_windowing_T win = (pf_windowing_T)winNo; + const int N = WIN_SIZE; + + for (int k = 0; k < N; ++k) + win_coeff_trig[k] = 1.0F; + //printf("calculate trig win ..\n"); + pf_window_trig(win_coeff_trig, win, N, 1.0F, w0); + + for (int k = 0; k < N; ++k) + win_coeff_cordic[k] = 1.0F; + //printf("calculate cordic win ..\n"); + const struct pf_windowing_param_tag *wp = pf_window_alloc(win, N, w0, 0.16F); + pf_window(win_coeff_cordic, wp); + + double wind_sum = pf_window_sum(wp); + int trig_max_idx = 0, cordic_max_idx = 0, delta_max_idx = 0; + float trig_max = win_coeff_trig[0]; + float cordic_max = win_coeff_cordic[0]; + float delta_max = win_coeff_cordic[0] - win_coeff_trig[0]; + for (int k = 0; k < N; ++k) + { + delta[k] = std::abs(win_coeff_cordic[k] - win_coeff_trig[k]); + if (delta_max < delta[k]) { + delta_max = delta[k]; + delta_max_idx = k; + } + if (trig_max < win_coeff_trig[k]) { + trig_max = win_coeff_trig[k]; + trig_max_idx = k; + } + if (cordic_max < win_coeff_cordic[k]) { + cordic_max = win_coeff_cordic[k]; + cordic_max_idx = k; + } + } + + printf("\ntrig_max: %f @ %d\tcordic_max: %f @ %d\tdelta_max: %g @ %d\tfor %s w0 = %d\n", + trig_max, trig_max_idx, cordic_max, cordic_max_idx, delta_max, delta_max_idx, + pf_window_text(win), w0 ); + printf("sum: %g\t0: %g / %g\tN-1: %g / %g\n", wind_sum + , win_coeff_trig[0], win_coeff_cordic[0] + , win_coeff_trig[N-1], win_coeff_cordic[N-1] ); + + printf("trig: "); + for (int k=0; k < 8; ++k) + printf("%d: %.2f\t", k * (N-1) / 7, win_coeff_trig[k * (N-1) / 7]); + printf("\n"); + printf("cordic: "); + for (int k=0; k < 8; ++k) + printf("%d: %.2f\t", k * (N-1) / 7, win_coeff_cordic[k * (N-1) / 7]); + printf("\n"); + + if (delta_max > ERR_LIMIT) + { + printf("trig and cordic windows coefficient differ too much: %g @ idx %d -> test failed!\n", delta_max, delta_max_idx ); + return 1; + } + + int N_div_i = 0; + pf_get_window_param(wp, &N_div_i); + + int expected_peak_idx = (w0 == 0) ? 0 : (N / 2 -1); + int expected_peak_idx2 = (w0 == 0) ? (N -1) : expected_peak_idx; + // && N_div_i == (N -1) + if (trig_max_idx != expected_peak_idx && trig_max_idx != expected_peak_idx2 && win != pf_winRect) + { + printf("trig: expected maximum at index %d for w0 = %d: failed!\n", expected_peak_idx, w0); + return 1; + } + if (cordic_max_idx != expected_peak_idx && cordic_max_idx != expected_peak_idx2 && win != pf_winRect) + { + printf("cordic: expected maximum at index %d for w0 = %d: failed!\n", expected_peak_idx, w0); + return 1; + } + } + } + + return 0; +}