v1.0

.gitignore

+./CMakeFiles
+./CMakeCache.txt
+./Makefile
+./cmake_install.cmake
+

BDS-3-Clause

+Valid-License-Identifier: BSD-3-Clause
+SPDX-URL: https://spdx.org/licenses/BSD-3-Clause.html
+Usage-Guide:
+  To use the BSD 3-clause "New" or "Revised" License put the following SPDX
+  tag/value pair into a comment according to the placement guidelines in
+  the licensing rules documentation:
+    SPDX-License-Identifier: BSD-3-Clause
+License-Text:
+
+Copyright (c) <year> <owner> . All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice,
+this list of conditions and the following disclaimer.
+
+2. 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.
+
+3. 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 THE COPYRIGHT HOLDER OR CONTRIBUTORS 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

CMakeLists.txt

+cmake_minimum_required(VERSION 3.1)
+project(NS)
+
+set(CMAKE_C_FLAGS "-std=c99")
+set(LIBRARY_OUTPUT_PATH     ${PROJECT_SOURCE_DIR}/bin)
+set(EXECUTABLE_OUTPUT_PATH  ${PROJECT_SOURCE_DIR}/bin)
+
+#add_definitions(-g -ggdb -Wall)
+#add_definitions(-Ddebug)
+
+aux_source_directory(${PROJECT_SOURCE_DIR}/src/kiss_fft KISSFFT_SRC_FILES)
+
+add_executable(NS_demo2 src/context.c src/blms.c src/blms_fdaf_main.c ${KISSFFT_SRC_FILES})
+target_link_libraries(NS_demo2 sndfile  m)
+
+#add_executable(test_kissfft src/test_kissfft.c ${KISSFFT_SRC_FILES})

README.md

-# blms
+# BLMS 
 
-block lms filter实现
\ No newline at end of file
+block lms filter implementation
+
+## v1.0 version
+implementation for the BLMS Frequence-domian Adapt filtering
+
+# Build
+The inner build system supported by BLMS:
++ CMAKE >=3.1
+
+## Dependency
++ Kissfft (BSD-3)
+
+To generate the Makefile
+```
+cmake .
+```
+and then run the make
+```
+make .
+```
+
+#License
+BSD-3-Clause

src/blms.c

+/*
+ *   Copyright (c) 2021-2022, Chris Zhang Jin. All rights reserved.
+ *   This file is part of blms - https://codechina.csdn.net/mimiduck/blms
+ *  
+ *   SPDX-License-Identifier: BSD-3 License
+ *   See COPYING file for more information.
+ *       
+ */
+#include "blms.h"
+
+/***
+ *   fft_buffer construct   [old] [new]
+ *   adapt                  [w]   [0....0]
+ *   total size N = 2*w_len | <-  N   -> |
+ * @param fft_buffer
+ * @param input
+ * @param dx
+ * @param w
+ * @param w_len
+ * @param out
+ * @return
+ */
+int BLMS_FDAF(context_t *context, float *input, float *dx, float *w, int w_len, float *out) {
+    // move the old to the left part
+    memmove(context->fft_buffer, context->fft_buffer+w_len, w_len* sizeof(float));
+
+    // copy the new data
+    memcpy(context->fft_buffer+w_len, input, w_len* sizeof(float));
+
+    // copy the w data to left part
+    memcpy(context->w_buffer, w, w_len* sizeof(float));
+
+    const int kLen = w_len*2;
+    // FFT for x
+    kiss_fft_cpx x[kLen];
+    memset(x, 0, sizeof(kiss_fft_cpx)*kLen);
+    for (int i=0; i<kLen; i++) {
+        x[i].r = context->fft_buffer[i];
+    }
+    kiss_fft_cpx X[kLen];
+    kiss_fft(context->fft_cfg, x, X);
+
+    // FFT for w
+    kiss_fft_cpx w_h[kLen];
+    for (int i=0; i<kLen; i++) {
+        w_h[i].r = context->w_buffer[i];
+        w_h[i].i = 0;
+    }
+    kiss_fft_cpx W[kLen];
+    kiss_fft(context->fft_cfg, w_h, W);
+
+    // complex multiply
+    kiss_fft_cpx Y[kLen];
+    for (int i=0; i<kLen; i++) {
+        Y[i].r = (X[i].r * W[i].r - X[i].i * W[i].i) / (kLen * kLen);
+        Y[i].i = (X[i].r * W[i].i + X[i].i * W[i].r) / (kLen * kLen);
+    }
+
+    kiss_fft_cpx o[kLen];
+    kiss_fft(context->ifft_cfg, Y, o);
+
+    // get the last N of ifft result
+    for (int i=w_len; i<kLen; i++) {
+        out[i-w_len] = o[i].r;
+        out[i-w_len] *= kLen;
+#ifdef debug
+        fprintf(context->file_o, "%4.4f,  ", out[i-w_len]);
+#endif
+    }
+#ifdef debug
+    fprintf(context->file_o, "\n");
+#endif
+
+    // conjugate X for corelation
+    for (int i=0; i<kLen; i++) {
+        X[i].i = -X[i].i;
+    }
+
+    // padding with 0 for length 2N
+    kiss_fft_cpx e[kLen];
+    for (int i=0; i<w_len; i++) {
+        e[i].r = 0;
+        e[i].i = 0;
+    }
+    for (int i=w_len; i<kLen; i++) {
+        e[i].r = dx[i-w_len] - o[i].r;
+#ifdef debug
+        fprintf(context->file_e, "%4.4f,  ", e[i].r);
+#endif
+        e[i].i = 0;
+    }
+#ifdef debug
+    fprintf(context->file_e, "\n");
+#endif
+
+    kiss_fft_cpx E[kLen];
+    kiss_fft(context->fft_cfg, e, E);
+
+    // complex multiply2
+    kiss_fft_cpx M2[kLen];
+    for (int i=0; i<kLen; i++) {
+        M2[i].r = (X[i].r*E[i].r - X[i].i*E[i].i)/(kLen*kLen);
+        M2[i].i = (X[i].r*E[i].i + X[i].i*E[i].r)/(kLen*kLen);
+    }
+
+    // ifft for the w's gradient
+    kiss_fft_cpx o2[kLen];
+    kiss_fft(context->ifft_cfg, M2, o2);
+
+    // keep the first N of output
+    for (int i=0; i<w_len; i++) {
+        o2[i].i = 0;
+    }
+    for (int i=w_len; i<kLen; i++) {
+        o2[i].i = 0;
+        o2[i].r = 0;
+    }
+
+    // move the first N backwards
+    for (int i=0; i<w_len; i++) {
+        o2[i+w_len].r = o2[i].r;
+        o2[i].r = 0;
+    }
+
+    // fft for the output of gradient
+    kiss_fft_cpx G[kLen];
+    kiss_fft(context->fft_cfg, o2, G);
+
+    // update for the w in freq domain
+    for (int i=0; i<kLen; i++) {
+        W[i].r = (W[i].r + 2*context->mu*G[i].r)/kLen;
+        W[i].i = (W[i].i + 2*context->mu*G[i].i)/kLen;
+    }
+
+    // ifft the updated w
+    kiss_fft_cpx wo[kLen];
+    kiss_fft(context->ifft_cfg, W, wo);
+    // last N
+    for (int i=w_len; i<kLen; i++) {
+        w[i-w_len] = wo[i].r;
+#ifdef debug
+        fprintf(context->file_w, "%4.4f,  ", wo[i].r);
+#endif
+    }
+#ifdef debug
+    fprintf(context->file_w, "\n");
+#endif
+    // first N
+//    for (int i=0; i<w_len; i++) {
+//        w[i] = wo[i].r;
+//        fprintf(context->file_w, "%4.4f,  ", wo[i].r);
+//    }
+//    fprintf(context->file_w, "\n");
+
+    return 0;
+}

src/blms.h

+/*
+ *  Copyright (c) 2021-2022, Chris Zhang Jin. All rights reserved.
+ *  This file is part of blms - https://codechina.csdn.net/mimiduck/blms
+ *  
+ *  SPDX-License-Identifier: BSD-3 License
+ *  See COPYING file for more information.
+ *       
+ */
+#ifndef NS_BLMS_H
+#define NS_BLMS_H
+
+#include "context.h"
+
+int BLMS_FDAF(context_t *context, float *input, float *dx, float* w, int w_len, float *out);
+
+#endif //NS_BLMS_H

src/blms_fdaf_main.c

+/*
+ *  Copyright (c) 2021-2022, Chris Zhang Jin. All rights reserved.
+ *  This file is part of blms - https://codechina.csdn.net/mimiduck/blms
+ *    
+ *  SPDX-License-Identifier: BSD-3 License
+ *  See COPYING file for more information.
+ *         
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <sndfile.h>
+#include "blms.h"
+
+
+#define FRAME_SIZE     W_LEN
+
+
+void process_file(context_t *context) {
+    SF_INFO info;
+    SNDFILE *in = sf_open(context->input_file_name, SFM_READ, &info);
+    if (!in) {
+        fprintf(stderr, "cannot open file %s\n", context->input_file_name);
+        return;
+    }
+
+    SF_INFO info2;
+    SNDFILE *dst = sf_open(context->destination_file_name, SFM_READ, &info2);
+    if (!dst) {
+        fprintf(stderr, "cannot open file %s\n", context->destination_file_name);
+        return;
+    }
+    int len = 1000;
+    int len_step = 1000;
+    float *dest_samples = (float *) calloc(len, sizeof(float));
+    sf_count_t cnt = 0;
+    do {
+        cnt = sf_read_float(dst, dest_samples+(len-len_step), len_step);
+        len += cnt;
+        if (cnt == 0) break;
+
+        dest_samples = realloc(dest_samples, len* sizeof(float));
+    } while (1);
+    sf_close(dst);
+
+    SF_INFO onfo;
+    onfo.channels = 1;
+    onfo.samplerate = 8000;
+    onfo.format = SF_FORMAT_PCM_16 | SF_FORMAT_WAV;
+    onfo.sections = 1;
+    onfo.seekable = 1;
+    SNDFILE *out = sf_open(context->output_file_name, SFM_WRITE, &onfo);
+    if (!out) {
+        fprintf(stderr, "cannot open file %s\n", context->output_file_name);
+        return;
+    }
+
+    // initialize adaptive filter coefficients, by zeros
+    float w[W_LEN];
+    for (int i = 0; i < W_LEN; i++) {
+        w[i] = 0.0f;
+    }
+
+    float frame[FRAME_SIZE] = {0};
+    sf_count_t read = 0;
+    int dni = 0;
+    int fi = 0;
+    do {
+        read = sf_read_float(in, frame, FRAME_SIZE);
+        if (read == 0) break;
+        printf("processing frame %d\n", fi++);
+
+        BLMS_FDAF(context, frame,  dest_samples+dni, w, W_LEN, frame);
+        dni = (dni+read)%len;
+
+        sf_write_float(out, frame, read);
+    } while (read != 0);
+
+
+    sf_close(in);
+    sf_close(out);
+    free(dest_samples);
+}
+
+
+int main(int argc, char* const argv[]) {
+    if (argc != 4) {
+        fprintf(stderr, "please input the file,output file,destination file \n");
+        exit(1);
+    }
+
+    printf("the input file name: %s\nthe output file name: %s\nthe destination file name: %s", argv[1], argv[2], argv[3]);
+
+    context_t *context = init_context(argv[1], argv[2], argv[3]);
+
+    process_file(context);
+
+    uninit_context(context);
+    printf("===end===\n");
+
+    return 0;
+}

src/context.c

+/*
+ *  Copyright (c) 2021-2022, Chris Zhang Jin. All rights reserved.
+ *  This file is part of blms - https://codechina.csdn.net/mimiduck/blms
+ *    
+ *  SPDX-License-Identifier: BSD-3 License
+ *  See COPYING file for more information.
+ *       
+ */
+
+#include "context.h"
+
+context_t *init_context(const char *input_file_name, const char *output_file_name, const char *destination_file_name) {
+    context_t *context = calloc(1, sizeof(context_t));
+
+    context->input_file_name = input_file_name;
+    context->output_file_name = output_file_name;
+    context->destination_file_name = destination_file_name;
+    context->fft_cfg = kiss_fft_alloc(W_LEN*2, 0, 0, 0);
+    context->ifft_cfg = kiss_fft_alloc(W_LEN*2, 1, 0, 0);
+    context->fft_buffer = calloc(W_LEN*2, sizeof(float));
+    context->w_buffer = calloc(W_LEN*2, sizeof(float));
+    context->mu = 0.2;
+#ifdef debug
+    context->file_w = fopen("file_w.txt", "w");
+    context->file_o = fopen("file_o.txt", "w");
+    context->file_e = fopen("file_e.txt", "w");
+#endif
+
+    return context;
+}
+
+void uninit_context(context_t *context) {
+    kiss_fft_free(context->fft_cfg);
+    kiss_fft_free(context->ifft_cfg);
+
+    free(context->fft_buffer);
+    free(context->w_buffer);
+#ifdef debug
+    fclose(context->file_w);
+    fclose(context->file_o);
+    fclose(context->file_e);
+#endif
+
+    free(context);
+}

src/context.h

+/*
+ *  Copyright (c) 2021-2022, Chris Zhang Jin. All rights reserved.
+ *  This file is part of blms - https://codechina.csdn.net/mimiduck/blms
+ *   
+ *  SPDX-License-Identifier: BSD-3 License
+ *  See COPYING file for more information.
+ *        
+ */
+
+#ifndef NS_CONTEXT_H
+#define NS_CONTEXT_H
+
+#include "kiss_fft/kiss_fft.h"
+
+#define W_LEN    120
+
+typedef struct {
+    const char              *input_file_name;
+    const char              *output_file_name;
+    const char              *destination_file_name;
+    float                   *fft_buffer;
+    float                   *w_buffer;
+    kiss_fft_cfg            fft_cfg;
+    kiss_fft_cfg            ifft_cfg;
+    float mu;
+#ifdef debug
+    FILE *file_w;
+    FILE *file_o;
+    FILE *file_e;
+#endif
+} context_t;
+
+context_t *init_context(const char *input_file_name, const char *output_file_name, const char *destination_file_name);
+
+void uninit_context(context_t *context);
+
+#endif //NS_CONTEXT_H

src/kiss_fft/_kiss_fft_guts.h

+/*
+ *  Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
+ *  This file is part of KISS FFT - https://github.com/mborgerding/kissfft
+ *
+ *  SPDX-License-Identifier: BSD-3-Clause
+ *  See COPYING file for more information.
+ */
+
+/* kiss_fft.h
+   defines kiss_fft_scalar as either short or a float type
+   and defines
+   typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */
+
+#ifndef _kiss_fft_guts_h
+#define _kiss_fft_guts_h
+
+#include "kiss_fft.h"
+#include <limits.h>
+
+#define MAXFACTORS 32
+/* e.g. an fft of length 128 has 4 factors 
+ as far as kissfft is concerned
+ 4*4*4*2
+ */
+
+struct kiss_fft_state{
+    int nfft;
+    int inverse;
+    int factors[2*MAXFACTORS];
+    kiss_fft_cpx twiddles[1];
+};
+
+/*
+  Explanation of macros dealing with complex math:
+
+   C_MUL(m,a,b)         : m = a*b
+   C_FIXDIV( c , div )  : if a fixed point impl., c /= div. noop otherwise
+   C_SUB( res, a,b)     : res = a - b
+   C_SUBFROM( res , a)  : res -= a
+   C_ADDTO( res , a)    : res += a
+ * */
+#ifdef FIXED_POINT
+#include <stdint.h>
+#if (FIXED_POINT==32)
+# define FRACBITS 31
+# define SAMPPROD int64_t
+#define SAMP_MAX INT32_MAX
+#define SAMP_MIN INT32_MIN
+#else
+# define FRACBITS 15
+# define SAMPPROD int32_t
+#define SAMP_MAX INT16_MAX
+#define SAMP_MIN INT16_MIN
+#endif
+
+#if defined(CHECK_OVERFLOW)
+#  define CHECK_OVERFLOW_OP(a,op,b)  \
+    if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
+        fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) );  }
+#endif
+
+
+#   define smul(a,b) ( (SAMPPROD)(a)*(b) )
+#   define sround( x )  (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS )
+
+#   define S_MUL(a,b) sround( smul(a,b) )
+
+#   define C_MUL(m,a,b) \
+      do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
+          (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
+
+#   define DIVSCALAR(x,k) \
+    (x) = sround( smul(  x, SAMP_MAX/k ) )
+
+#   define C_FIXDIV(c,div) \
+    do {    DIVSCALAR( (c).r , div);  \
+        DIVSCALAR( (c).i  , div); }while (0)
+
+#   define C_MULBYSCALAR( c, s ) \
+    do{ (c).r =  sround( smul( (c).r , s ) ) ;\
+        (c).i =  sround( smul( (c).i , s ) ) ; }while(0)
+
+#else  /* not FIXED_POINT*/
+
+#   define S_MUL(a,b) ( (a)*(b) )
+#define C_MUL(m,a,b) \
+    do{ (m).r = (a).r*(b).r - (a).i*(b).i;\
+        (m).i = (a).r*(b).i + (a).i*(b).r; }while(0)
+#   define C_FIXDIV(c,div) /* NOOP */
+#   define C_MULBYSCALAR( c, s ) \
+    do{ (c).r *= (s);\
+        (c).i *= (s); }while(0)
+#endif
+
+#ifndef CHECK_OVERFLOW_OP
+#  define CHECK_OVERFLOW_OP(a,op,b) /* noop */
+#endif
+
+#define  C_ADD( res, a,b)\
+    do { \
+        CHECK_OVERFLOW_OP((a).r,+,(b).r)\
+        CHECK_OVERFLOW_OP((a).i,+,(b).i)\
+        (res).r=(a).r+(b).r;  (res).i=(a).i+(b).i; \
+    }while(0)
+#define  C_SUB( res, a,b)\
+    do { \
+        CHECK_OVERFLOW_OP((a).r,-,(b).r)\
+        CHECK_OVERFLOW_OP((a).i,-,(b).i)\
+        (res).r=(a).r-(b).r;  (res).i=(a).i-(b).i; \
+    }while(0)
+#define C_ADDTO( res , a)\
+    do { \
+        CHECK_OVERFLOW_OP((res).r,+,(a).r)\
+        CHECK_OVERFLOW_OP((res).i,+,(a).i)\
+        (res).r += (a).r;  (res).i += (a).i;\
+    }while(0)
+
+#define C_SUBFROM( res , a)\
+    do {\
+        CHECK_OVERFLOW_OP((res).r,-,(a).r)\
+        CHECK_OVERFLOW_OP((res).i,-,(a).i)\
+        (res).r -= (a).r;  (res).i -= (a).i; \
+    }while(0)
+
+
+#ifdef FIXED_POINT
+#  define KISS_FFT_COS(phase)  floor(.5+SAMP_MAX * cos (phase))
+#  define KISS_FFT_SIN(phase)  floor(.5+SAMP_MAX * sin (phase))
+#  define HALF_OF(x) ((x)>>1)
+#elif defined(USE_SIMD)
+#  define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) )
+#  define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) )
+#  define HALF_OF(x) ((x)*_mm_set1_ps(.5))
+#else
+#  define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase)
+#  define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase)
+#  define HALF_OF(x) ((x)*.5)
+#endif
+
+#define  kf_cexp(x,phase) \
+    do{ \
+        (x)->r = KISS_FFT_COS(phase);\
+        (x)->i = KISS_FFT_SIN(phase);\
+    }while(0)
+
+
+/* a debugging function */
+#define pcpx(c)\
+    fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) )
+
+
+#ifdef KISS_FFT_USE_ALLOCA
+// define this to allow use of alloca instead of malloc for temporary buffers
+// Temporary buffers are used in two case: 
+// 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5
+// 2. "in-place" FFTs.  Notice the quotes, since kissfft does not really do an in-place transform.
+#include <alloca.h>
+#define  KISS_FFT_TMP_ALLOC(nbytes) alloca(nbytes)
+#define  KISS_FFT_TMP_FREE(ptr) 
+#else
+#define  KISS_FFT_TMP_ALLOC(nbytes) KISS_FFT_MALLOC(nbytes)
+#define  KISS_FFT_TMP_FREE(ptr) KISS_FFT_FREE(ptr)
+#endif
+
+#endif /* _kiss_fft_guts_h */
+

src/kiss_fft/kiss_fft.c

+/*
+ *  Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
+ *  This file is part of KISS FFT - https://github.com/mborgerding/kissfft
+ *
+ *  SPDX-License-Identifier: BSD-3-Clause
+ *  See COPYING file for more information.
+ */
+
+
+#include "_kiss_fft_guts.h"
+/* The guts header contains all the multiplication and addition macros that are defined for
+ fixed or floating point complex numbers.  It also delares the kf_ internal functions.
+ */
+
+static void kf_bfly2(
+        kiss_fft_cpx * Fout,
+        const size_t fstride,
+        const kiss_fft_cfg st,
+        int m
+        )
+{
+    kiss_fft_cpx * Fout2;
+    kiss_fft_cpx * tw1 = st->twiddles;
+    kiss_fft_cpx t;
+    Fout2 = Fout + m;
+    do{
+        C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2);
+
+        C_MUL (t,  *Fout2 , *tw1);
+        tw1 += fstride;
+        C_SUB( *Fout2 ,  *Fout , t );
+        C_ADDTO( *Fout ,  t );
+        ++Fout2;
+        ++Fout;
+    }while (--m);
+}
+
+static void kf_bfly4(
+        kiss_fft_cpx * Fout,
+        const size_t fstride,
+        const kiss_fft_cfg st,
+        const size_t m
+        )
+{
+    kiss_fft_cpx *tw1,*tw2,*tw3;
+    kiss_fft_cpx scratch[6];
+    size_t k=m;
+    const size_t m2=2*m;
+    const size_t m3=3*m;
+
+
+    tw3 = tw2 = tw1 = st->twiddles;
+
+    do {
+        C_FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4);
+
+        C_MUL(scratch[0],Fout[m] , *tw1 );
+        C_MUL(scratch[1],Fout[m2] , *tw2 );
+        C_MUL(scratch[2],Fout[m3] , *tw3 );
+
+        C_SUB( scratch[5] , *Fout, scratch[1] );
+        C_ADDTO(*Fout, scratch[1]);
+        C_ADD( scratch[3] , scratch[0] , scratch[2] );
+        C_SUB( scratch[4] , scratch[0] , scratch[2] );
+        C_SUB( Fout[m2], *Fout, scratch[3] );
+        tw1 += fstride;
+        tw2 += fstride*2;
+        tw3 += fstride*3;
+        C_ADDTO( *Fout , scratch[3] );
+
+        if(st->inverse) {
+            Fout[m].r = scratch[5].r - scratch[4].i;
+            Fout[m].i = scratch[5].i + scratch[4].r;
+            Fout[m3].r = scratch[5].r + scratch[4].i;
+            Fout[m3].i = scratch[5].i - scratch[4].r;
+        }else{
+            Fout[m].r = scratch[5].r + scratch[4].i;
+            Fout[m].i = scratch[5].i - scratch[4].r;
+            Fout[m3].r = scratch[5].r - scratch[4].i;
+            Fout[m3].i = scratch[5].i + scratch[4].r;
+        }
+        ++Fout;
+    }while(--k);
+}
+
+static void kf_bfly3(
+         kiss_fft_cpx * Fout,
+         const size_t fstride,
+         const kiss_fft_cfg st,
+         size_t m
+         )
+{
+     size_t k=m;
+     const size_t m2 = 2*m;
+     kiss_fft_cpx *tw1,*tw2;
+     kiss_fft_cpx scratch[5];
+     kiss_fft_cpx epi3;
+     epi3 = st->twiddles[fstride*m];
+
+     tw1=tw2=st->twiddles;
+
+     do{
+         C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
+
+         C_MUL(scratch[1],Fout[m] , *tw1);
+         C_MUL(scratch[2],Fout[m2] , *tw2);
+
+         C_ADD(scratch[3],scratch[1],scratch[2]);
+         C_SUB(scratch[0],scratch[1],scratch[2]);
+         tw1 += fstride;
+         tw2 += fstride*2;
+
+         Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
+         Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
+
+         C_MULBYSCALAR( scratch[0] , epi3.i );
+
+         C_ADDTO(*Fout,scratch[3]);
+
+         Fout[m2].r = Fout[m].r + scratch[0].i;
+         Fout[m2].i = Fout[m].i - scratch[0].r;
+
+         Fout[m].r -= scratch[0].i;
+         Fout[m].i += scratch[0].r;
+
+         ++Fout;
+     }while(--k);
+}
+
+static void kf_bfly5(
+        kiss_fft_cpx * Fout,
+        const size_t fstride,
+        const kiss_fft_cfg st,
+        int m
+        )
+{
+    kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
+    int u;
+    kiss_fft_cpx scratch[13];
+    kiss_fft_cpx * twiddles = st->twiddles;
+    kiss_fft_cpx *tw;
+    kiss_fft_cpx ya,yb;
+    ya = twiddles[fstride*m];
+    yb = twiddles[fstride*2*m];
+
+    Fout0=Fout;
+    Fout1=Fout0+m;
+    Fout2=Fout0+2*m;
+    Fout3=Fout0+3*m;
+    Fout4=Fout0+4*m;
+
+    tw=st->twiddles;
+    for ( u=0; u<m; ++u ) {
+        C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
+        scratch[0] = *Fout0;
+
+        C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
+        C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
+        C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
+        C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
+
+        C_ADD( scratch[7],scratch[1],scratch[4]);
+        C_SUB( scratch[10],scratch[1],scratch[4]);
+        C_ADD( scratch[8],scratch[2],scratch[3]);
+        C_SUB( scratch[9],scratch[2],scratch[3]);
+
+        Fout0->r += scratch[7].r + scratch[8].r;
+        Fout0->i += scratch[7].i + scratch[8].i;
+
+        scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
+        scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
+
+        scratch[6].r =  S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
+        scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);
+
+        C_SUB(*Fout1,scratch[5],scratch[6]);
+        C_ADD(*Fout4,scratch[5],scratch[6]);
+
+        scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
+        scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
+        scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
+        scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);
+
+        C_ADD(*Fout2,scratch[11],scratch[12]);
+        C_SUB(*Fout3,scratch[11],scratch[12]);
+
+        ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
+    }
+}
+
+/* perform the butterfly for one stage of a mixed radix FFT */
+static void kf_bfly_generic(
+        kiss_fft_cpx * Fout,
+        const size_t fstride,
+        const kiss_fft_cfg st,
+        int m,
+        int p
+        )
+{
+    int u,k,q1,q;
+    kiss_fft_cpx * twiddles = st->twiddles;
+    kiss_fft_cpx t;
+    int Norig = st->nfft;
+
+    kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p);
+
+    for ( u=0; u<m; ++u ) {
+        k=u;
+        for ( q1=0 ; q1<p ; ++q1 ) {
+            scratch[q1] = Fout[ k  ];
+            C_FIXDIV(scratch[q1],p);
+            k += m;
+        }
+
+        k=u;
+        for ( q1=0 ; q1<p ; ++q1 ) {
+            int twidx=0;
+            Fout[ k ] = scratch[0];
+            for (q=1;q<p;++q ) {
+                twidx += fstride * k;
+                if (twidx>=Norig) twidx-=Norig;
+                C_MUL(t,scratch[q] , twiddles[twidx] );
+                C_ADDTO( Fout[ k ] ,t);
+            }
+            k += m;
+        }
+    }
+    KISS_FFT_TMP_FREE(scratch);
+}
+
+static
+void kf_work(
+        kiss_fft_cpx * Fout,
+        const kiss_fft_cpx * f,
+        const size_t fstride,
+        int in_stride,
+        int * factors,
+        const kiss_fft_cfg st
+        )
+{
+    kiss_fft_cpx * Fout_beg=Fout;
+    const int p=*factors++; /* the radix  */
+    const int m=*factors++; /* stage's fft length/p */
+    const kiss_fft_cpx * Fout_end = Fout + p*m;
+
+#ifdef _OPENMP
+    // use openmp extensions at the 
+    // top-level (not recursive)
+    if (fstride==1 && p<=5 && m!=1)
+    {
+        int k;
+
+        // execute the p different work units in different threads
+#       pragma omp parallel for
+        for (k=0;k<p;++k) 
+            kf_work( Fout +k*m, f+ fstride*in_stride*k,fstride*p,in_stride,factors,st);
+        // all threads have joined by this point
+
+        switch (p) {
+            case 2: kf_bfly2(Fout,fstride,st,m); break;
+            case 3: kf_bfly3(Fout,fstride,st,m); break; 
+            case 4: kf_bfly4(Fout,fstride,st,m); break;
+            case 5: kf_bfly5(Fout,fstride,st,m); break; 
+            default: kf_bfly_generic(Fout,fstride,st,m,p); break;
+        }
+        return;
+    }
+#endif
+
+    if (m==1) {
+        do{
+            *Fout = *f;
+            f += fstride*in_stride;
+        }while(++Fout != Fout_end );
+    }else{
+        do{
+            // recursive call:
+            // DFT of size m*p performed by doing
+            // p instances of smaller DFTs of size m, 
+            // each one takes a decimated version of the input
+            kf_work( Fout , f, fstride*p, in_stride, factors,st);
+            f += fstride*in_stride;
+        }while( (Fout += m) != Fout_end );
+    }
+
+    Fout=Fout_beg;
+
+    // recombine the p smaller DFTs 
+    switch (p) {
+        case 2: kf_bfly2(Fout,fstride,st,m); break;
+        case 3: kf_bfly3(Fout,fstride,st,m); break; 
+        case 4: kf_bfly4(Fout,fstride,st,m); break;
+        case 5: kf_bfly5(Fout,fstride,st,m); break; 
+        default: kf_bfly_generic(Fout,fstride,st,m,p); break;
+    }
+}
+
+/*  facbuf is populated by p1,m1,p2,m2, ...
+    where 
+    p[i] * m[i] = m[i-1]
+    m0 = n                  */
+static 
+void kf_factor(int n,int * facbuf)
+{
+    int p=4;
+    double floor_sqrt;
+    floor_sqrt = floor( sqrt((double)n) );
+
+    /*factor out powers of 4, powers of 2, then any remaining primes */
+    do {
+        while (n % p) {
+            switch (p) {
+                case 4: p = 2; break;
+                case 2: p = 3; break;
+                default: p += 2; break;
+            }
+            if (p > floor_sqrt)
+                p = n;          /* no more factors, skip to end */
+        }
+        n /= p;
+        *facbuf++ = p;
+        *facbuf++ = n;
+    } while (n > 1);
+}
+
+/*
+ *
+ * User-callable function to allocate all necessary storage space for the fft.
+ *
+ * The return value is a contiguous block of memory, allocated with malloc.  As such,
+ * It can be freed with free(), rather than a kiss_fft-specific function.
+ * */
+kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem )
+{
+    KISS_FFT_ALIGN_CHECK(mem)
+
+    kiss_fft_cfg st=NULL;
+    size_t memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fft_state)
+        + sizeof(kiss_fft_cpx)*(nfft-1)); /* twiddle factors*/
+
+    if ( lenmem==NULL ) {
+        st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
+    }else{
+        if (mem != NULL && *lenmem >= memneeded)
+            st = (kiss_fft_cfg)mem;
+        *lenmem = memneeded;
+    }
+    if (st) {
+        int i;
+        st->nfft=nfft;
+        st->inverse = inverse_fft;
+
+        for (i=0;i<nfft;++i) {
+            const double pi=3.141592653589793238462643383279502884197169399375105820974944;
+            double phase = -2*pi*i / nfft;
+            if (st->inverse)
+                phase *= -1;
+            kf_cexp(st->twiddles+i, phase );
+        }
+
+        kf_factor(nfft,st->factors);
+    }
+    return st;
+}
+
+
+void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
+{
+    if (fin == fout) {
+        //NOTE: this is not really an in-place FFT algorithm.
+        //It just performs an out-of-place FFT into a temp buffer
+        kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft);
+        kf_work(tmpbuf,fin,1,in_stride, st->factors,st);
+        memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft);
+        KISS_FFT_TMP_FREE(tmpbuf);
+    }else{
+        kf_work( fout, fin, 1,in_stride, st->factors,st );
+    }
+}
+
+void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
+{
+    kiss_fft_stride(cfg,fin,fout,1);
+}
+
+
+void kiss_fft_cleanup(void)
+{
+    // nothing needed any more
+}
+
+int kiss_fft_next_fast_size(int n)
+{
+    while(1) {
+        int m=n;
+        while ( (m%2) == 0 ) m/=2;
+        while ( (m%3) == 0 ) m/=3;
+        while ( (m%5) == 0 ) m/=5;
+        if (m<=1)
+            break; /* n is completely factorable by twos, threes, and fives */
+        n++;
+    }
+    return n;
+}

src/kiss_fft/kiss_fft.h

+/*
+ *  Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
+ *  This file is part of KISS FFT - https://github.com/mborgerding/kissfft
+ *
+ *  SPDX-License-Identifier: BSD-3-Clause
+ *  See COPYING file for more information.
+ */
+
+#ifndef KISS_FFT_H
+#define KISS_FFT_H
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*
+ ATTENTION!
+ If you would like a :
+ -- a utility that will handle the caching of fft objects
+ -- real-only (no imaginary time component ) FFT
+ -- a multi-dimensional FFT
+ -- a command-line utility to perform ffts
+ -- a command-line utility to perform fast-convolution filtering
+
+ Then see kfc.h kiss_fftr.h kiss_fftnd.h fftutil.c kiss_fastfir.c
+  in the tools/ directory.
+*/
+
+/* User may override KISS_FFT_MALLOC and/or KISS_FFT_FREE. */
+#ifdef USE_SIMD
+# include <xmmintrin.h>
+# define kiss_fft_scalar __m128
+# ifndef KISS_FFT_MALLOC
+#  define KISS_FFT_MALLOC(nbytes) _mm_malloc(nbytes,16)
+#  define KISS_FFT_ALIGN_CHECK(ptr) 
+#  define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 15UL) & ~0xFUL)
+# endif
+# ifndef KISS_FFT_FREE
+#  define KISS_FFT_FREE _mm_free
+# endif
+#else
+# define KISS_FFT_ALIGN_CHECK(ptr)
+# define KISS_FFT_ALIGN_SIZE_UP(size) (size)
+# ifndef KISS_FFT_MALLOC
+#  define KISS_FFT_MALLOC malloc
+# endif
+# ifndef KISS_FFT_FREE
+#  define KISS_FFT_FREE free
+# endif
+#endif
+
+
+#ifdef FIXED_POINT
+#include <stdint.h>
+# if (FIXED_POINT == 32)
+#  define kiss_fft_scalar int32_t
+# else	
+#  define kiss_fft_scalar int16_t
+# endif
+#else
+# ifndef kiss_fft_scalar
+/*  default is float */
+#   define kiss_fft_scalar float
+# endif
+#endif
+
+typedef struct {
+    kiss_fft_scalar r;
+    kiss_fft_scalar i;
+}kiss_fft_cpx;
+
+typedef struct kiss_fft_state* kiss_fft_cfg;
+
+/* 
+ *  kiss_fft_alloc
+ *  
+ *  Initialize a FFT (or IFFT) algorithm's cfg/state buffer.
+ *
+ *  typical usage:      kiss_fft_cfg mycfg=kiss_fft_alloc(1024,0,NULL,NULL);
+ *
+ *  The return value from fft_alloc is a cfg buffer used internally
+ *  by the fft routine or NULL.
+ *
+ *  If lenmem is NULL, then kiss_fft_alloc will allocate a cfg buffer using malloc.
+ *  The returned value should be free()d when done to avoid memory leaks.
+ *  
+ *  The state can be placed in a user supplied buffer 'mem':
+ *  If lenmem is not NULL and mem is not NULL and *lenmem is large enough,
+ *      then the function places the cfg in mem and the size used in *lenmem
+ *      and returns mem.
+ *  
+ *  If lenmem is not NULL and ( mem is NULL or *lenmem is not large enough),
+ *      then the function returns NULL and places the minimum cfg 
+ *      buffer size in *lenmem.
+ * */
+
+kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem); 
+
+/*
+ * kiss_fft(cfg,in_out_buf)
+ *
+ * Perform an FFT on a complex input buffer.
+ * for a forward FFT,
+ * fin should be  f[0] , f[1] , ... ,f[nfft-1]
+ * fout will be   F[0] , F[1] , ... ,F[nfft-1]
+ * Note that each element is complex and can be accessed like
+    f[k].r and f[k].i
+ * */
+void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout);
+
+/*
+ A more generic version of the above function. It reads its input from every Nth sample.
+ * */
+void kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride);
+
+/* If kiss_fft_alloc allocated a buffer, it is one contiguous 
+   buffer and can be simply free()d when no longer needed*/
+#define kiss_fft_free KISS_FFT_FREE
+
+/*
+ Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up 
+ your compiler output to call this before you exit.
+*/
+void kiss_fft_cleanup(void);
+	
+
+/*
+ * Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5)
+ */
+int kiss_fft_next_fast_size(int n);
+
+/* for real ffts, we need an even size */
+#define kiss_fftr_next_fast_size_real(n) \
+        (kiss_fft_next_fast_size( ((n)+1)>>1)<<1)
+
+#ifdef __cplusplus
+} 
+#endif
+
+#endif

src/kiss_fft/kiss_fftr.c

+/*
+ *  Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
+ *  This file is part of KISS FFT - https://github.com/mborgerding/kissfft
+ *
+ *  SPDX-License-Identifier: BSD-3-Clause
+ *  See COPYING file for more information.
+ */
+
+#include "kiss_fftr.h"
+#include "_kiss_fft_guts.h"
+
+struct kiss_fftr_state{
+    kiss_fft_cfg substate;
+    kiss_fft_cpx * tmpbuf;
+    kiss_fft_cpx * super_twiddles;
+#ifdef USE_SIMD
+    void * pad;
+#endif
+};
+
+kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)
+{
+	KISS_FFT_ALIGN_CHECK(mem)
+
+    int i;
+    kiss_fftr_cfg st = NULL;
+    size_t subsize = 0, memneeded;
+
+    if (nfft & 1) {
+        fprintf(stderr,"Real FFT optimization must be even.\n");
+        return NULL;
+    }
+    nfft >>= 1;
+
+    kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize);
+    memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2);
+
+    if (lenmem == NULL) {
+        st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded);
+    } else {
+        if (*lenmem >= memneeded)
+            st = (kiss_fftr_cfg) mem;
+        *lenmem = memneeded;
+    }
+    if (!st)
+        return NULL;
+
+    st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */
+    st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize);
+    st->super_twiddles = st->tmpbuf + nfft;
+    kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize);
+
+    for (i = 0; i < nfft/2; ++i) {
+        double phase =
+            -3.14159265358979323846264338327 * ((double) (i+1) / nfft + .5);
+        if (inverse_fft)
+            phase *= -1;
+        kf_cexp (st->super_twiddles+i,phase);
+    }
+    return st;
+}
+
+void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata)
+{
+    /* input buffer timedata is stored row-wise */
+    int k,ncfft;
+    kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc;
+
+    if ( st->substate->inverse) {
+        fprintf(stderr,"kiss fft usage error: improper alloc\n");
+        exit(1);
+    }
+
+    ncfft = st->substate->nfft;
+
+    /*perform the parallel fft of two real signals packed in real,imag*/
+    kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf );
+    /* The real part of the DC element of the frequency spectrum in st->tmpbuf
+     * contains the sum of the even-numbered elements of the input time sequence
+     * The imag part is the sum of the odd-numbered elements
+     *
+     * The sum of tdc.r and tdc.i is the sum of the input time sequence. 
+     *      yielding DC of input time sequence
+     * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... 
+     *      yielding Nyquist bin of input time sequence
+     */
+ 
+    tdc.r = st->tmpbuf[0].r;
+    tdc.i = st->tmpbuf[0].i;
+    C_FIXDIV(tdc,2);
+    CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i);
+    CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i);
+    freqdata[0].r = tdc.r + tdc.i;
+    freqdata[ncfft].r = tdc.r - tdc.i;
+#ifdef USE_SIMD    
+    freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0);
+#else
+    freqdata[ncfft].i = freqdata[0].i = 0;
+#endif
+
+    for ( k=1;k <= ncfft/2 ; ++k ) {
+        fpk    = st->tmpbuf[k]; 
+        fpnk.r =   st->tmpbuf[ncfft-k].r;
+        fpnk.i = - st->tmpbuf[ncfft-k].i;
+        C_FIXDIV(fpk,2);
+        C_FIXDIV(fpnk,2);
+
+        C_ADD( f1k, fpk , fpnk );
+        C_SUB( f2k, fpk , fpnk );
+        C_MUL( tw , f2k , st->super_twiddles[k-1]);
+
+        freqdata[k].r = HALF_OF(f1k.r + tw.r);
+        freqdata[k].i = HALF_OF(f1k.i + tw.i);
+        freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r);
+        freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i);
+    }
+}
+
+void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
+{
+    /* input buffer timedata is stored row-wise */
+    int k, ncfft;
+
+    if (st->substate->inverse == 0) {
+        fprintf (stderr, "kiss fft usage error: improper alloc\n");
+        exit (1);
+    }
+
+    ncfft = st->substate->nfft;
+
+    st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r;
+    st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r;
+    C_FIXDIV(st->tmpbuf[0],2);
+
+    for (k = 1; k <= ncfft / 2; ++k) {
+        kiss_fft_cpx fk, fnkc, fek, fok, tmp;
+        fk = freqdata[k];
+        fnkc.r = freqdata[ncfft - k].r;
+        fnkc.i = -freqdata[ncfft - k].i;
+        C_FIXDIV( fk , 2 );
+        C_FIXDIV( fnkc , 2 );
+
+        C_ADD (fek, fk, fnkc);
+        C_SUB (tmp, fk, fnkc);
+        C_MUL (fok, tmp, st->super_twiddles[k-1]);
+        C_ADD (st->tmpbuf[k],     fek, fok);
+        C_SUB (st->tmpbuf[ncfft - k], fek, fok);
+#ifdef USE_SIMD        
+        st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0);
+#else
+        st->tmpbuf[ncfft - k].i *= -1;
+#endif
+    }
+    kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata);
+}

src/kiss_fft/kiss_fftr.h

+/*
+ *  Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
+ *  This file is part of KISS FFT - https://github.com/mborgerding/kissfft
+ *
+ *  SPDX-License-Identifier: BSD-3-Clause
+ *  See COPYING file for more information.
+ */
+
+#ifndef KISS_FTR_H
+#define KISS_FTR_H
+
+#include "kiss_fft.h"
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+    
+/* 
+ 
+ Real optimized version can save about 45% cpu time vs. complex fft of a real seq.
+
+ 
+ 
+ */
+
+typedef struct kiss_fftr_state *kiss_fftr_cfg;
+
+
+kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenmem);
+/*
+ nfft must be even
+
+ If you don't care to allocate space, use mem = lenmem = NULL 
+*/
+
+
+void kiss_fftr(kiss_fftr_cfg cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata);
+/*
+ input timedata has nfft scalar points
+ output freqdata has nfft/2+1 complex points
+*/
+
+void kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata);
+/*
+ input freqdata has  nfft/2+1 complex points
+ output timedata has nfft scalar points
+*/
+
+#define kiss_fftr_free KISS_FFT_FREE
+
+#ifdef __cplusplus
+}
+#endif
+#endif

src/test_kissfft.c

+/*
+ *  Copyright (c) 2021-2022, Chris Zhang Jin. All rights reserved.
+ *  This file is part of blms - https://codechina.csdn.net/mimiduck/blms
+ *   
+ *  SPDX-License-Identifier: BSD-3 License
+ *  See COPYING file for more information.
+ *        
+ */
+
+#include "kiss_fft/kiss_fft.h"
+#include "kiss_fft/kiss_fftr.h"
+
+static kiss_fft_scalar rand_scalar(void) {
+    kiss_fft_scalar s = (kiss_fft_scalar) ((rand() + 10) % 256);
+    return s / 256.;
+}
+
+static void print_fft_result(kiss_fft_cpx *x, int n) {
+    int l = 0;
+    for (int i = 0; i < n; i++) {
+        printf("(%4.4f+%4.4fi), ", x[i].r, x[i].i);
+        if (l++ == 7) {
+            l = 0;
+            printf("\n");
+        }
+    }
+    printf("\n");
+}
+
+void test_fft() {
+    int i = 0;
+    int nfft = 20;
+    kiss_fft_cpx cin[nfft];
+    kiss_fft_cpx cout[nfft];
+    kiss_fft_cpx sout[nfft];
+    kiss_fft_cfg kiss_fft_state;
+
+    kiss_fft_scalar zero;
+    memset(&zero, 0, sizeof(zero));
+
+    for (i = 0; i < nfft; ++i) {
+        cin[i].r = rand_scalar();
+        cin[i].i = zero;
+    }
+    printf("\n");
+
+    printf(" init data for kiss_fft (cin): \n");
+    print_fft_result(cin, nfft);
+
+    memset(cout, 0, sizeof(short) * nfft);
+    memset(sout, 0, sizeof(short) * nfft);
+
+    kiss_fft_state = kiss_fft_alloc(nfft, 0, 0, 0);
+
+    kiss_fft(kiss_fft_state, cin, cout);
+    kiss_fft_free(kiss_fft_state);
+
+    printf(" results from kiss_fft (cout): \n");
+    print_fft_result(cout, nfft);
+
+    kiss_fft_state = kiss_fft_alloc(nfft, 1, 0, 0);
+
+    for (i = 0; i < nfft; i++) {
+        cout[i].r /= nfft;
+        cout[i].i /= nfft;
+    }
+    // case A.
+    kiss_fft(kiss_fft_state, cout, sout);
+
+    // end case
+    kiss_fft_free(kiss_fft_state);
+
+    printf(" results from kiss_ifft (sout): \n");
+    print_fft_result(sout, nfft);
+    printf("\n");
+}
+
+void test_fftr() {
+    int i = 0;
+    int nfft = 20;
+    kiss_fftr_cfg kiss_fftr_state;
+
+    kiss_fft_scalar rin[nfft];
+    kiss_fft_scalar rout[nfft];
+    kiss_fft_cpx out[nfft];
+    memset(rin, 0, sizeof(short) * nfft);
+    memset(rout, 0, sizeof(short) * nfft);
+    memset(out, 0, sizeof(kiss_fft_cpx) * nfft);
+    kiss_fft_scalar zero;
+    memset(&zero, 0, sizeof(zero));
+
+    for (i = 0; i < nfft; ++i) {
+        rin[i] = rand_scalar();
+    }
+//    for (; i < nfft * 2; ++i) {
+//        rin[i] = 0;
+//    }
+
+    printf(" init data for kiss_fft (rin): \n");
+    for (i = 0; i < nfft; i++) {
+        printf("%4.4f,", rin[i]);
+    }
+    printf("\n");
+
+
+    kiss_fftr_state = kiss_fftr_alloc(nfft, 0, 0, 0);
+
+    kiss_fftr(kiss_fftr_state, rin, out);
+    kiss_fftr_free(kiss_fftr_state);
+
+    printf(" results from kiss_fft (out): \n");
+    print_fft_result(out, nfft);
+
+    kiss_fftr_state = kiss_fftr_alloc(nfft, 1, 0, 0);
+
+    for (i = 0; i < nfft; i++) {
+        out[i].r /= nfft;
+        out[i].i /= nfft;
+    }
+    // case A.
+    kiss_fftri(kiss_fftr_state, out, rout);
+
+    // end case
+    kiss_fftr_free(kiss_fftr_state);
+
+
+    printf(" results from kiss_ifft (rout): \n");
+    for (i = 0; i < nfft; i++) {
+        printf("%4.4f,", rout[i]);
+    }
+    printf("\n");
+}
+
+int main(int argc, char *const argv[]) {
+
+    test_fftr();
+
+    return 0;
+}