First commit

- Method to get parameters of one harmonic

16 files changed, 1703 insertions, 0 deletions

diff --git a/src/3rd/kissfft/test/Makefile b/src/3rd/kissfft/test/Makefile
new file mode 100644
index 0000000..c204511
--- /dev/null
+++ b/src/3rd/kissfft/test/Makefile
@@ -0,0 +1,108 @@
+
+WARNINGS=-W -Wall -Wstrict-prototypes -Wmissing-prototypes -Waggregate-return \
+    -Wcast-align -Wcast-qual -Wnested-externs -Wshadow -Wbad-function-cast \
+    -Wwrite-strings
+
+CFLAGS=-O3 -I.. -I../tools $(WARNINGS)
+CFLAGS+=-ffast-math -fomit-frame-pointer 
+#CFLAGS+=-funroll-loops
+#CFLAGS+=-march=prescott 
+#CFLAGS+= -mtune=native 
+# TIP: try adding -openmp or -fopenmp  to enable OPENMP directives and use of multiple cores
+#CFLAGS+=-fopenmp
+CFLAGS+= $(CFLAGADD)
+
+
+ifeq "$(NFFT)" ""
+ NFFT=1800
+endif
+ifeq "$(NUMFFTS)" ""
+ NUMFFTS=10000
+endif
+
+ifeq "$(DATATYPE)" ""
+ DATATYPE=float
+endif
+
+BENCHKISS=bm_kiss_$(DATATYPE)
+BENCHFFTW=bm_fftw_$(DATATYPE)
+SELFTEST=st_$(DATATYPE)
+TESTREAL=tr_$(DATATYPE)
+TESTKFC=tkfc_$(DATATYPE)
+FASTFILTREAL=ffr_$(DATATYPE)
+SELFTESTSRC=twotonetest.c
+
+
+TYPEFLAGS=-Dkiss_fft_scalar=$(DATATYPE)
+
+ifeq  "$(DATATYPE)" "int16_t" 
+ TYPEFLAGS=-DFIXED_POINT=16
+endif
+
+ifeq  "$(DATATYPE)" "int32_t"
+ TYPEFLAGS=-DFIXED_POINT=32 
+endif
+
+ifeq  "$(DATATYPE)" "simd"
+ TYPEFLAGS=-DUSE_SIMD=1 -msse
+endif
+
+
+ifeq  "$(DATATYPE)" "float"
+ # fftw needs to be built with --enable-float to build this lib
+ FFTWLIB=-lfftw3f
+else
+ FFTWLIB=-lfftw3
+endif
+
+FFTWLIBDIR=-L/usr/local/lib/
+
+SRCFILES=../kiss_fft.c ../tools/kiss_fftnd.c ../tools/kiss_fftr.c pstats.c ../tools/kfc.c ../tools/kiss_fftndr.c
+
+all: tools $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL) $(TESTKFC)
+
+tools:
+	cd ../tools && make all
+
+
+$(SELFTEST): $(SELFTESTSRC) $(SRCFILES)
+	$(CC) -o $@ $(CFLAGS) $(TYPEFLAGS) $+ -lm 
+
+$(TESTKFC): $(SRCFILES) 
+	$(CC) -o $@ $(CFLAGS)  -DKFC_TEST $(TYPEFLAGS) $+ -lm
+	
+$(TESTREAL): test_real.c $(SRCFILES)
+	$(CC) -o $@ $(CFLAGS) $(TYPEFLAGS) $+ -lm
+
+$(BENCHKISS): benchkiss.c $(SRCFILES)
+	$(CC) -o $@ $(CFLAGS) $(TYPEFLAGS)  $+ -lm
+
+$(BENCHFFTW): benchfftw.c pstats.c
+	@echo "======attempting to build FFTW benchmark"
+	@$(CC) -o $@ $(CFLAGS) -DDATATYPE$(DATATYPE) $+ $(FFTWLIB) $(FFTWLIBDIR) -lm || echo "FFTW not available for comparison"
+
+test: all
+	@./$(TESTKFC)
+	@echo "======1d & 2-d complex fft self test (type= $(DATATYPE) )"
+	@./$(SELFTEST)
+	@echo "======real FFT (type= $(DATATYPE) )"
+	@./$(TESTREAL)
+	@echo "======timing test (type=$(DATATYPE))"
+	@./$(BENCHKISS) -x $(NUMFFTS) -n $(NFFT) 
+	@[ -x ./$(BENCHFFTW) ] && ./$(BENCHFFTW) -x $(NUMFFTS) -n $(NFFT) ||true
+	@echo "======higher dimensions type=$(DATATYPE))"
+	@./testkiss.py
+
+selftest.c:
+	./mk_test.py 10 12 14 > selftest.c
+selftest_short.c:
+	./mk_test.py -s 10 12 14 > selftest_short.c
+
+
+CXXFLAGS=-O3 -ffast-math -fomit-frame-pointer  -I.. -I../tools -W -Wall
+testcpp: testcpp.cc ../kissfft.hh
+	$(CXX) -o $@ $(CXXFLAGS) testcpp.cc -lm
+
+
+clean:
+	rm -f *~ bm_* st_* tr_* kf_* tkfc_* ff_* ffr_* *.pyc *.pyo *.dat testcpp
diff --git a/src/3rd/kissfft/test/benchfftw.c b/src/3rd/kissfft/test/benchfftw.c
new file mode 100644
index 0000000..0bb2ef1
--- /dev/null
+++ b/src/3rd/kissfft/test/benchfftw.c
@@ -0,0 +1,101 @@
+/*
+ *  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 <stdio.h>
+#include <stdlib.h>
+#include <fftw3.h>
+#include <unistd.h>
+#include "pstats.h"
+
+#ifdef DATATYPEdouble
+
+#define CPXTYPE fftw_complex
+#define PLAN fftw_plan
+#define FFTMALLOC fftw_malloc
+#define MAKEPLAN fftw_plan_dft_1d
+#define DOFFT fftw_execute
+#define DESTROYPLAN fftw_destroy_plan
+#define FFTFREE fftw_free
+
+#elif defined(DATATYPEfloat)
+
+#define CPXTYPE fftwf_complex
+#define PLAN fftwf_plan
+#define FFTMALLOC fftwf_malloc
+#define MAKEPLAN fftwf_plan_dft_1d
+#define DOFFT fftwf_execute
+#define DESTROYPLAN fftwf_destroy_plan
+#define FFTFREE fftwf_free
+
+#endif
+
+#ifndef CPXTYPE
+int main(void)
+{
+    fprintf(stderr,"Datatype not available in FFTW\n" );
+    return 0;
+}
+#else
+int main(int argc,char ** argv)
+{
+    int nfft=1024;
+    int isinverse=0;
+    int numffts=1000,i;
+
+    CPXTYPE * in=NULL;
+    CPXTYPE * out=NULL;
+    PLAN p;
+
+    pstats_init();
+
+    while (1) {
+      int c = getopt (argc, argv, "n:ix:h");
+      if (c == -1)
+        break;
+      switch (c) {
+      case 'n':
+        nfft = atoi (optarg);
+        break;
+      case 'x':
+        numffts = atoi (optarg);
+        break;
+      case 'i':
+        isinverse = 1;
+        break;
+      case 'h':
+      case '?':
+      default:
+        fprintf(stderr,"options:\n-n N: complex fft length\n-i: inverse\n-x N: number of ffts to compute\n"
+                "");
+      }
+    }
+
+    in=FFTMALLOC(sizeof(CPXTYPE) * nfft);
+    out=FFTMALLOC(sizeof(CPXTYPE) * nfft);
+    for (i=0;i<nfft;++i ) {
+        in[i][0] = rand() - RAND_MAX/2;
+        in[i][1] = rand() - RAND_MAX/2;
+    }
+
+    if ( isinverse )
+        p = MAKEPLAN(nfft, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
+    else    
+        p = MAKEPLAN(nfft, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
+
+    for (i=0;i<numffts;++i)
+        DOFFT(p);
+
+    DESTROYPLAN(p);
+
+    FFTFREE(in); FFTFREE(out);
+
+    fprintf(stderr,"fftw\tnfft=%d\tnumffts=%d\n", nfft,numffts);
+    pstats_report();
+
+    return 0;
+}
+#endif
diff --git a/src/3rd/kissfft/test/benchkiss.c b/src/3rd/kissfft/test/benchkiss.c
new file mode 100644
index 0000000..7dae65d
--- /dev/null
+++ b/src/3rd/kissfft/test/benchkiss.c
@@ -0,0 +1,129 @@
+/*
+ *  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 <stdio.h>
+#include <stdlib.h>
+#include <sys/times.h>
+#include <unistd.h>
+#include "kiss_fft.h"
+#include "kiss_fftr.h"
+#include "kiss_fftnd.h"
+#include "kiss_fftndr.h"
+
+#include "pstats.h"
+
+static
+int getdims(int * dims, char * arg)
+{
+    char *s;
+    int ndims=0;
+    while ( (s=strtok( arg , ",") ) ) {
+        dims[ndims++] = atoi(s);
+        //printf("%s=%d\n",s,dims[ndims-1]);
+        arg=NULL;
+    }
+    return ndims;
+}
+
+int main(int argc,char ** argv)
+{
+    int k;
+    int nfft[32];
+    int ndims = 1;
+    int isinverse=0;
+    int numffts=1000,i;
+    kiss_fft_cpx * buf;
+    kiss_fft_cpx * bufout;
+    int real = 0;
+
+    nfft[0] = 1024;// default
+
+    while (1) {
+        int c = getopt (argc, argv, "n:ix:r");
+        if (c == -1)
+            break;
+        switch (c) {
+            case 'r':
+                real = 1;
+                break;
+            case 'n':
+                ndims = getdims(nfft, optarg );
+                if (nfft[0] != kiss_fft_next_fast_size(nfft[0]) ) {
+                    int ng = kiss_fft_next_fast_size(nfft[0]);
+                    fprintf(stderr,"warning: %d might be a better choice for speed than %d\n",ng,nfft[0]);
+                }
+                break;
+            case 'x':
+                numffts = atoi (optarg);
+                break;
+            case 'i':
+                isinverse = 1;
+                break;
+        }
+    }
+    int nbytes = sizeof(kiss_fft_cpx);
+    for (k=0;k<ndims;++k)
+        nbytes *= nfft[k];
+
+#ifdef USE_SIMD        
+    numffts /= 4;
+    fprintf(stderr,"since SIMD implementation does 4 ffts at a time, numffts is being reduced to %d\n",numffts);
+#endif
+
+    buf=(kiss_fft_cpx*)KISS_FFT_MALLOC(nbytes);
+    bufout=(kiss_fft_cpx*)KISS_FFT_MALLOC(nbytes);
+    memset(buf,0,nbytes);
+
+    pstats_init();
+
+    if (ndims==1) {
+        if (real) {
+            kiss_fftr_cfg st = kiss_fftr_alloc( nfft[0] ,isinverse ,0,0);
+            if (isinverse)
+                for (i=0;i<numffts;++i)
+                    kiss_fftri( st ,(kiss_fft_cpx*)buf,(kiss_fft_scalar*)bufout );
+            else
+                for (i=0;i<numffts;++i)
+                    kiss_fftr( st ,(kiss_fft_scalar*)buf,(kiss_fft_cpx*)bufout );
+            free(st);
+        }else{
+            kiss_fft_cfg st = kiss_fft_alloc( nfft[0] ,isinverse ,0,0);
+            for (i=0;i<numffts;++i)
+                kiss_fft( st ,buf,bufout );
+            free(st);
+        }
+    }else{
+        if (real) {
+            kiss_fftndr_cfg st = kiss_fftndr_alloc( nfft,ndims ,isinverse ,0,0);
+            if (isinverse)
+                for (i=0;i<numffts;++i)
+                    kiss_fftndri( st ,(kiss_fft_cpx*)buf,(kiss_fft_scalar*)bufout );
+            else
+                for (i=0;i<numffts;++i)
+                    kiss_fftndr( st ,(kiss_fft_scalar*)buf,(kiss_fft_cpx*)bufout );
+            free(st);
+        }else{
+            kiss_fftnd_cfg st= kiss_fftnd_alloc(nfft,ndims,isinverse ,0,0);
+            for (i=0;i<numffts;++i)
+                kiss_fftnd( st ,buf,bufout );
+            free(st);
+        }
+    }
+
+    free(buf); free(bufout);
+
+    fprintf(stderr,"KISS\tnfft=");
+    for (k=0;k<ndims;++k)
+        fprintf(stderr, "%d,",nfft[k]);
+    fprintf(stderr,"\tnumffts=%d\n" ,numffts);
+    pstats_report();
+
+    kiss_fft_cleanup();
+
+    return 0;
+}
+
diff --git a/src/3rd/kissfft/test/compfft.py b/src/3rd/kissfft/test/compfft.py
new file mode 100755
index 0000000..d2671c1
--- /dev/null
+++ b/src/3rd/kissfft/test/compfft.py
@@ -0,0 +1,97 @@
+#!/usr/bin/env python
+#  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.
+
+# use FFTPACK as a baseline
+import FFT
+from Numeric import *
+import math
+import random
+import sys
+import struct
+import fft
+
+pi=math.pi
+e=math.e
+j=complex(0,1)
+lims=(-32768,32767)
+
+def randbuf(n,cpx=1):
+    res = array( [ random.uniform( lims[0],lims[1] ) for i in range(n) ] )
+    if cpx:
+        res = res + j*randbuf(n,0)
+    return res
+
+def main():
+    from getopt import getopt
+    import popen2
+    opts,args = getopt( sys.argv[1:],'u:n:Rt:' )
+    opts=dict(opts)
+    exitcode=0
+
+    util = opts.get('-u','./kf_float')
+
+    try:
+        dims = [ int(d) for d in opts['-n'].split(',')]
+        cpx = opts.get('-R') is None
+        fmt=opts.get('-t','f')
+    except KeyError:
+        sys.stderr.write("""
+        usage: compfft.py 
+        -n d1[,d2,d3...]  : FFT dimension(s)
+        -u utilname : see sample_code/fftutil.c, default = ./kf_float
+        -R : real-optimized version\n""")
+        sys.exit(1)
+
+    x = fft.make_random( dims )
+
+    cmd = '%s -n %s ' % ( util, ','.join([ str(d) for d in dims]) )
+    if cpx:
+        xout = FFT.fftnd(x)
+        xout = reshape(xout,(size(xout),))
+    else:
+        cmd += '-R '
+        xout = FFT.real_fft(x)
+
+    proc = popen2.Popen3( cmd , bufsize=len(x) )
+
+    proc.tochild.write( dopack( x , fmt ,cpx ) )
+    proc.tochild.close()
+    xoutcomp = dounpack( proc.fromchild.read( ) , fmt ,1 )
+    #xoutcomp = reshape( xoutcomp , dims )
+
+    sig = xout * conjugate(xout)
+    sigpow = sum( sig )
+
+    diff = xout-xoutcomp
+    noisepow = sum( diff * conjugate(diff) )
+
+    snr = 10 * math.log10(abs( sigpow / noisepow ) )
+    if snr<100:
+        print xout
+        print xoutcomp
+        exitcode=1
+    print 'NFFT=%s,SNR = %f dB' % (str(dims),snr)
+    sys.exit(exitcode)
+
+def dopack(x,fmt,cpx):
+    x = reshape( x, ( size(x),) )
+    if cpx:
+        s = ''.join( [ struct.pack('ff',c.real,c.imag) for c in x ] )
+    else:
+        s = ''.join( [ struct.pack('f',c) for c in x ] )
+    return s 
+
+def dounpack(x,fmt,cpx):
+    uf = fmt * ( len(x) / 4 )
+    s = struct.unpack(uf,x)
+    if cpx:
+        return array(s[::2]) + array( s[1::2] )*j
+    else:    
+        return array(s )
+
+if __name__ == "__main__":
+    main()
diff --git a/src/3rd/kissfft/test/doit.c b/src/3rd/kissfft/test/doit.c
new file mode 100644
index 0000000..36e42e8
--- /dev/null
+++ b/src/3rd/kissfft/test/doit.c
@@ -0,0 +1,134 @@
+/*
+ *  This program is in the public domain 
+ *  A program that helps the authors of the fine fftw library benchmark kiss 
+ *  This file is part of KISS FFT - https://github.com/mborgerding/kissfft
+ *
+ *  SPDX-License-Identifier: Unlicense
+ *  See LICENSES/Unlicense for more information.
+ */
+
+#include "bench-user.h"
+#include <math.h>
+
+#include "kiss_fft.h"
+#include "kiss_fftnd.h"
+#include "kiss_fftr.h"
+
+BEGIN_BENCH_DOC
+BENCH_DOC("name", "kissfft")
+BENCH_DOC("version", "1.0.1")
+BENCH_DOC("year", "2004")
+BENCH_DOC("author", "Mark Borgerding")
+BENCH_DOC("language", "C")
+BENCH_DOC("url", "http://sourceforge.net/projects/kissfft/")
+BENCH_DOC("copyright",
+"Copyright (c) 2003,4 Mark Borgerding\n"
+"\n"
+"All rights reserved.\n"
+"\n"
+"Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:\n"
+"\n"
+"    * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.\n"
+"    * 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.\n"
+"    * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.\n"
+"\n"
+	  "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 OWNER 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.\n")
+END_BENCH_DOC
+
+int can_do(struct problem *p)
+{
+    if (p->rank == 1) {
+        if (p->kind == PROBLEM_REAL) {
+            return (p->n[0] & 1) == 0;  /* only even real is okay */
+        } else {
+            return 1;
+        }
+    } else {
+        return p->kind == PROBLEM_COMPLEX;
+    }
+}
+
+static kiss_fft_cfg cfg=NULL;
+static kiss_fftr_cfg cfgr=NULL;
+static kiss_fftnd_cfg cfgnd=NULL;
+
+#define FAILIF( c ) \
+    if (  c ) do {\
+    fprintf(stderr,\
+    "kissfft: " #c " (file=%s:%d errno=%d %s)\n",\
+        __FILE__,__LINE__ , errno,strerror( errno ) ) ;\
+    exit(1);\
+    }while(0)
+
+
+
+void setup(struct problem *p)
+{
+    size_t i;
+
+    /*
+    fprintf(stderr,"%s %s %d-d ",
+            (p->sign == 1)?"Inverse":"Forward",
+            p->kind == PROBLEM_COMPLEX?"Complex":"Real",
+            p->rank);
+    */
+    if (p->rank == 1) {
+        if (p->kind == PROBLEM_COMPLEX) {
+            cfg = kiss_fft_alloc (p->n[0], (p->sign == 1), 0, 0);
+            FAILIF(cfg==NULL);
+        }else{
+            cfgr = kiss_fftr_alloc (p->n[0], (p->sign == 1), 0, 0);
+            FAILIF(cfgr==NULL);
+        }
+    }else{
+        int dims[5];
+        for (i=0;i<p->rank;++i){
+            dims[i] = p->n[i];
+        }
+        /* multi-dimensional */
+        if (p->kind == PROBLEM_COMPLEX) {
+            cfgnd = kiss_fftnd_alloc( dims , p->rank, (p->sign == 1), 0, 0 );
+            FAILIF(cfgnd==NULL);
+        }
+    }
+}
+
+void doit(int iter, struct problem *p)
+{
+    int i;
+    void *in = p->in;
+    void *out = p->out;
+
+    if (p->in_place)
+        out = p->in;
+
+    if (p->rank == 1) {
+        if (p->kind == PROBLEM_COMPLEX){
+            for (i = 0; i < iter; ++i)
+                kiss_fft (cfg, in, out);
+        } else {
+            /* PROBLEM_REAL */
+            if (p->sign == -1)   /* FORWARD */
+                for (i = 0; i < iter; ++i)
+                    kiss_fftr (cfgr, in, out);
+            else
+                for (i = 0; i < iter; ++i)
+                    kiss_fftri (cfgr, in, out);
+        }
+    }else{
+        /* multi-dimensional */
+        for (i = 0; i < iter; ++i)
+            kiss_fftnd(cfgnd,in,out);
+    }
+}
+
+void done(struct problem *p)
+{
+     free(cfg);
+     cfg=NULL;
+     free(cfgr);
+     cfgr=NULL;
+     free(cfgnd);
+     cfgnd=NULL;
+     UNUSED(p);
+}
diff --git a/src/3rd/kissfft/test/fastfir.py b/src/3rd/kissfft/test/fastfir.py
new file mode 100755
index 0000000..18662d4
--- /dev/null
+++ b/src/3rd/kissfft/test/fastfir.py
@@ -0,0 +1,107 @@
+#!/usr/bin/env python
+#  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.
+
+from Numeric import *
+from FFT import *
+
+def make_random(len):
+    import random
+    res=[]
+    for i in range(int(len)):
+        r=random.uniform(-1,1)
+        i=random.uniform(-1,1)
+        res.append( complex(r,i) )
+    return res
+
+def slowfilter(sig,h):
+    translen = len(h)-1
+    return convolve(sig,h)[translen:-translen]
+
+def nextpow2(x):
+    return 2 ** math.ceil(math.log(x)/math.log(2))
+
+def fastfilter(sig,h,nfft=None):
+    if nfft is None:
+        nfft = int( nextpow2( 2*len(h) ) )
+    H = fft( h , nfft )
+    scraplen = len(h)-1
+    keeplen = nfft-scraplen
+    res=[]
+    isdone = 0
+    lastidx = nfft
+    idx0 = 0
+    while not isdone:
+        idx1 = idx0 + nfft
+        if idx1 >= len(sig):
+            idx1 = len(sig)
+            lastidx = idx1-idx0
+            if lastidx <= scraplen:
+                break
+            isdone = 1
+        Fss = fft(sig[idx0:idx1],nfft)
+        fm = Fss * H
+        m = inverse_fft(fm)
+        res.append( m[scraplen:lastidx] )
+        idx0 += keeplen
+    return concatenate( res )
+
+def main():
+    import sys
+    from getopt import getopt
+    opts,args = getopt(sys.argv[1:],'rn:l:')
+    opts=dict(opts)
+
+    siglen = int(opts.get('-l',1e4 ) )
+    hlen =50 
+ 
+    nfft = int(opts.get('-n',128) )
+    usereal = opts.has_key('-r')
+
+    print 'nfft=%d'%nfft
+    # make a signal
+    sig = make_random( siglen )
+    # make an impulse response
+    h = make_random( hlen )
+    #h=[1]*2+[0]*3
+    if usereal:
+        sig=[c.real for c in sig]
+        h=[c.real for c in h]
+
+    # perform MAC filtering
+    yslow = slowfilter(sig,h)
+    #print '<YSLOW>',yslow,'</YSLOW>'
+    #yfast = fastfilter(sig,h,nfft)
+    yfast = utilfastfilter(sig,h,nfft,usereal)
+    #print yfast
+    print 'len(yslow)=%d'%len(yslow)
+    print 'len(yfast)=%d'%len(yfast)
+    diff = yslow-yfast
+    snr = 10*log10( abs( vdot(yslow,yslow) / vdot(diff,diff) ) )
+    print 'snr=%s' % snr
+    if snr < 10.0:
+        print 'h=',h
+        print 'sig=',sig[:5],'...'
+        print 'yslow=',yslow[:5],'...'
+        print 'yfast=',yfast[:5],'...'
+
+def utilfastfilter(sig,h,nfft,usereal):
+    import compfft
+    import os
+    open( 'sig.dat','w').write( compfft.dopack(sig,'f',not usereal) )
+    open( 'h.dat','w').write( compfft.dopack(h,'f',not usereal) )
+    if usereal: 
+        util = './fastconvr' 
+    else:
+        util = './fastconv'
+    cmd = 'time %s -n %d -i sig.dat -h h.dat -o out.dat' % (util, nfft)
+    print cmd
+    ec  = os.system(cmd)
+    print 'exited->',ec
+    return compfft.dounpack(open('out.dat').read(),'f',not usereal)
+
+if __name__ == "__main__":
+    main()
diff --git a/src/3rd/kissfft/test/fft.py b/src/3rd/kissfft/test/fft.py
new file mode 100755
index 0000000..4208a20
--- /dev/null
+++ b/src/3rd/kissfft/test/fft.py
@@ -0,0 +1,201 @@
+#!/usr/bin/env python
+#  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.
+
+import math
+import sys
+import random
+
+pi=math.pi
+e=math.e
+j=complex(0,1)
+
+def fft(f,inv):
+    n=len(f)
+    if n==1:
+        return f
+
+    for p in 2,3,5:
+        if n%p==0:
+            break
+    else:
+        raise Exception('%s not factorable ' % n)
+
+    m = n/p
+    Fout=[]
+    for q in range(p): # 0,1
+        fp = f[q::p]  # every p'th time sample
+        Fp = fft( fp ,inv)
+        Fout.extend( Fp )
+
+    for u in range(m):
+        scratch = Fout[u::m] # u to end in strides of m
+        for q1 in range(p):
+            k = q1*m + u  # indices to Fout above that became scratch
+            Fout[ k ] = scratch[0] # cuz e**0==1 in loop below
+            for q in range(1,p):
+                if inv:
+                    t = e ** ( j*2*pi*k*q/n )
+                else:                    
+                    t = e ** ( -j*2*pi*k*q/n )
+                Fout[ k ] += scratch[q] * t
+
+    return Fout
+
+def rifft(F):
+    N = len(F) - 1
+    Z = [0] * (N)
+    for k in range(N):
+        Fek = ( F[k] + F[-k-1].conjugate() )
+        Fok = ( F[k] - F[-k-1].conjugate() ) * e ** (j*pi*k/N)
+        Z[k] = Fek + j*Fok
+
+    fp = fft(Z , 1)
+
+    f = []
+    for c in fp:
+        f.append(c.real)
+        f.append(c.imag)
+    return f
+
+def real_fft( f,inv ):
+    if inv:
+        return rifft(f)
+
+    N = len(f) / 2
+
+    res = f[::2]
+    ims = f[1::2]
+
+    fp = [ complex(r,i) for r,i in zip(res,ims) ]
+    print 'fft input ', fp
+    Fp = fft( fp ,0 )
+    print 'fft output ', Fp
+
+    F = [ complex(0,0) ] * ( N+1 )
+    
+    F[0] = complex( Fp[0].real + Fp[0].imag , 0 ) 
+
+    for k in range(1,N/2+1):
+        tw = e ** ( -j*pi*(.5+float(k)/N ) )
+        
+        F1k = Fp[k] + Fp[N-k].conjugate()
+        F2k = Fp[k] - Fp[N-k].conjugate()
+        F2k *= tw
+        F[k] = ( F1k + F2k ) * .5
+        F[N-k] = ( F1k - F2k ).conjugate() * .5
+        #F[N-k] = ( F1kp + e ** ( -j*pi*(.5+float(N-k)/N ) ) * F2kp ) * .5
+        #F[N-k] = ( F1k.conjugate() - tw.conjugate() * F2k.conjugate() ) * .5
+
+    F[N] = complex( Fp[0].real - Fp[0].imag , 0 ) 
+    return F
+
+def main():
+    #fft_func = fft
+    fft_func = real_fft
+
+    tvec = [0.309655,0.815653,0.768570,0.591841,0.404767,0.637617,0.007803,0.012665]
+    Ftvec = [ complex(r,i) for r,i in zip(
+                [3.548571,-0.378761,-0.061950,0.188537,-0.566981,0.188537,-0.061950,-0.378761],
+                [0.000000,-1.296198,-0.848764,0.225337,0.000000,-0.225337,0.848764,1.296198] ) ]
+
+    F = fft_func( tvec,0 )
+
+    nerrs= 0
+    for i in range(len(Ftvec)/2 + 1):
+        if abs( F[i] - Ftvec[i] )> 1e-5:
+            print 'F[%d]: %s != %s' % (i,F[i],Ftvec[i])
+            nerrs += 1
+
+    print '%d errors in forward fft' % nerrs
+    if nerrs:
+        return
+
+    trec = fft_func( F , 1 )
+
+    for i in range(len(trec) ):
+        trec[i] /= len(trec)
+
+    for i in range(len(tvec) ):
+        if abs( trec[i] - tvec[i] )> 1e-5:
+            print 't[%d]: %s != %s' % (i,tvec[i],trec[i])
+            nerrs += 1
+
+    print '%d errors in reverse fft' % nerrs
+
+
+def make_random(dims=[1]):
+    import Numeric 
+    res = []
+    for i in range(dims[0]):
+        if len(dims)==1:
+            r=random.uniform(-1,1)
+            i=random.uniform(-1,1)
+            res.append( complex(r,i) )
+        else:
+            res.append( make_random( dims[1:] ) )
+    return Numeric.array(res)
+
+def flatten(x):
+    import Numeric
+    ntotal = Numeric.product(Numeric.shape(x))
+    return Numeric.reshape(x,(ntotal,))
+
+def randmat( ndims ):
+    dims=[]
+    for i in range( ndims ):
+        curdim = int( random.uniform(2,4) )
+        dims.append( curdim )
+    return make_random(dims )
+
+def test_fftnd(ndims=3):
+    import FFT
+    import Numeric
+
+    x=randmat( ndims )
+    print 'dimensions=%s' % str( Numeric.shape(x) )
+    #print 'x=%s' %str(x)
+    xver = FFT.fftnd(x)
+    x2=myfftnd(x)
+    err = xver - x2
+    errf = flatten(err)
+    xverf = flatten(xver)
+    errpow = Numeric.vdot(errf,errf)+1e-10
+    sigpow = Numeric.vdot(xverf,xverf)+1e-10
+    snr = 10*math.log10(abs(sigpow/errpow) )
+    if snr<80:
+        print xver
+        print x2
+    print 'SNR=%sdB' % str( snr )
+ 
+def myfftnd(x):
+    import Numeric
+    xf = flatten(x)
+    Xf = fftndwork( xf , Numeric.shape(x) )
+    return Numeric.reshape(Xf,Numeric.shape(x) )
+
+def fftndwork(x,dims):
+    import Numeric
+    dimprod=Numeric.product( dims )
+
+    for k in range( len(dims) ):
+        cur_dim=dims[ k ]
+        stride=dimprod/cur_dim
+        next_x = [complex(0,0)]*len(x)
+        for i in range(stride):
+            next_x[i*cur_dim:(i+1)*cur_dim] = fft(x[i:(i+cur_dim)*stride:stride],0)
+        x = next_x
+    return x
+
+if __name__ == "__main__":
+    try:
+        nd = int(sys.argv[1])
+    except:
+        nd=None
+    if nd:    
+        test_fftnd( nd )
+    else:    
+        sys.exit(0)
diff --git a/src/3rd/kissfft/test/mk_test.py b/src/3rd/kissfft/test/mk_test.py
new file mode 100755
index 0000000..ab6acdf
--- /dev/null
+++ b/src/3rd/kissfft/test/mk_test.py
@@ -0,0 +1,122 @@
+#!/usr/bin/env python
+#  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.
+
+import FFT
+import sys
+import random
+import re
+j=complex(0,1)
+
+def randvec(n,iscomplex):
+    if iscomplex:
+        return [
+                int(random.uniform(-32768,32767) ) + j*int(random.uniform(-32768,32767) )
+                for i in range(n) ]
+    else:                
+        return [ int(random.uniform(-32768,32767) ) for i in range(n) ]
+    
+def c_format(v,round=0):
+    if round:
+        return ','.join( [ '{%d,%d}' %(int(c.real),int(c.imag) ) for c in v ] ) 
+    else:
+        s= ','.join( [ '{%.60f ,%.60f }' %(c.real,c.imag) for c in v ] ) 
+        return re.sub(r'\.?0+ ',' ',s)
+
+def test_cpx( n,inverse ,short):
+    v = randvec(n,1)
+    scale = 1
+    if short:
+        minsnr=30
+    else:
+        minsnr=100
+
+    if inverse:
+        tvecout = FFT.inverse_fft(v)
+        if short:
+            scale = 1
+        else:            
+            scale = len(v)
+    else:
+        tvecout = FFT.fft(v)
+        if short:
+            scale = 1.0/len(v)
+
+    tvecout = [ c * scale for c in tvecout ]
+
+
+    s="""#define NFFT %d""" % len(v) + """
+    {
+        double snr;
+        kiss_fft_cpx test_vec_in[NFFT] = { """  + c_format(v) + """};
+        kiss_fft_cpx test_vec_out[NFFT] = {"""  + c_format( tvecout ) + """};
+        kiss_fft_cpx testbuf[NFFT];
+        void * cfg = kiss_fft_alloc(NFFT,%d,0,0);""" % inverse + """
+
+        kiss_fft(cfg,test_vec_in,testbuf);
+        snr = snr_compare(test_vec_out,testbuf,NFFT);
+        printf("DATATYPE=" xstr(kiss_fft_scalar) ", FFT n=%d, inverse=%d, snr = %g dB\\n",NFFT,""" + str(inverse) + """,snr);
+        if (snr<""" + str(minsnr) + """)
+            exit_code++;
+        free(cfg);
+    }
+#undef NFFT    
+"""
+    return s
+
+def compare_func():
+    s="""
+#define xstr(s) str(s)
+#define str(s) #s
+double snr_compare( kiss_fft_cpx * test_vec_out,kiss_fft_cpx * testbuf, int n)
+{
+    int k;
+    double sigpow,noisepow,err,snr,scale=0;
+    kiss_fft_cpx err;
+    sigpow = noisepow = .000000000000000000000000000001; 
+
+    for (k=0;k<n;++k) {
+        sigpow += test_vec_out[k].r * test_vec_out[k].r + 
+                  test_vec_out[k].i * test_vec_out[k].i;
+        C_SUB(err,test_vec_out[k],testbuf[k].r);
+        noisepow += err.r * err.r + err.i + err.i;
+
+        if (test_vec_out[k].r)
+            scale += testbuf[k].r / test_vec_out[k].r;
+    }
+    snr = 10*log10( sigpow / noisepow );
+    scale /= n;
+    if (snr<10)
+        printf( "\\npoor snr, try a scaling factor %f\\n" , scale );
+    return snr;
+}
+"""
+    return s
+
+def main():
+
+    from getopt import getopt
+    opts,args = getopt(sys.argv[1:],'s')
+    opts = dict(opts)
+    short = int( opts.has_key('-s') )
+
+    fftsizes = args
+    if not fftsizes:
+        fftsizes = [ 1800 ]
+    print '#include "kiss_fft.h"'
+    print compare_func()
+    print "int main() { int exit_code=0;\n"
+    for n in fftsizes:
+        n = int(n)
+        print test_cpx(n,0,short)
+        print test_cpx(n,1,short)
+    print """
+    return exit_code;
+}
+"""
+
+if __name__ == "__main__":
+    main()
diff --git a/src/3rd/kissfft/test/pstats.c b/src/3rd/kissfft/test/pstats.c
new file mode 100644
index 0000000..5082bda
--- /dev/null
+++ b/src/3rd/kissfft/test/pstats.c
@@ -0,0 +1,56 @@
+/*
+ *  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 <stdio.h>
+#include <stdlib.h>
+#include <sys/times.h>
+#include <sys/types.h>
+#include <unistd.h>
+
+#include "pstats.h"
+
+static struct tms tms_beg;
+static struct tms tms_end;
+static int has_times = 0;
+
+
+void pstats_init(void)
+{
+    has_times = times(&tms_beg) != -1;
+}
+
+static void tms_report(void)
+{
+    double cputime;
+    if (! has_times )
+        return;
+    times(&tms_end);
+    cputime = ( ((float)tms_end.tms_utime + tms_end.tms_stime + tms_end.tms_cutime + tms_end.tms_cstime ) -
+                ((float)tms_beg.tms_utime + tms_beg.tms_stime + tms_beg.tms_cutime + tms_beg.tms_cstime ) )
+               / sysconf(_SC_CLK_TCK);
+    fprintf(stderr,"\tcputime=%.3f\n" , cputime);
+}
+
+static void ps_report(void)
+{
+    char buf[1024];
+#ifdef __APPLE__ /*  MAC OS X */
+    sprintf(buf,"ps -o command,majflt,minflt,rss,pagein,vsz -p %d 1>&2",getpid() );
+#else /* GNU/Linux */
+    sprintf(buf,"ps -o comm,majflt,minflt,rss,drs,pagein,sz,trs,vsz %d 1>&2",getpid() );
+#endif    
+    if (system( buf )==-1) {
+        perror("system call to ps failed");
+    }
+}
+
+void pstats_report()
+{
+    ps_report();
+    tms_report();
+}
+
diff --git a/src/3rd/kissfft/test/pstats.h b/src/3rd/kissfft/test/pstats.h
new file mode 100644
index 0000000..04a304b
--- /dev/null
+++ b/src/3rd/kissfft/test/pstats.h
@@ -0,0 +1,14 @@
+/*
+ *  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 PSTATS_H
+#define PSTATS_H
+
+void pstats_init(void);
+void pstats_report(void);
+ 
+#endif
diff --git a/src/3rd/kissfft/test/tailscrap.m b/src/3rd/kissfft/test/tailscrap.m
new file mode 100644
index 0000000..abf9046
--- /dev/null
+++ b/src/3rd/kissfft/test/tailscrap.m
@@ -0,0 +1,26 @@
+function maxabsdiff=tailscrap()
+% test code for circular convolution with the scrapped portion 
+% at the tail of the buffer, rather than the front
+%
+% The idea is to rotate the zero-padded h (impulse response) buffer
+% to the left nh-1 samples, rotating the junk samples as well.
+% This could be very handy in avoiding buffer copies during fast filtering.
+nh=10;
+nfft=256;
+
+h=rand(1,nh);
+x=rand(1,nfft);
+
+hpad=[ h(nh) zeros(1,nfft-nh) h(1:nh-1) ]; 
+
+% baseline comparison
+y1 = filter(h,1,x);
+y1_notrans = y1(nh:nfft);
+
+% fast convolution
+y2 = ifft( fft(hpad) .* fft(x) );
+y2_notrans=y2(1:nfft-nh+1);
+
+maxabsdiff = max(abs(y2_notrans - y1_notrans))
+
+end
diff --git a/src/3rd/kissfft/test/test_real.c b/src/3rd/kissfft/test/test_real.c
new file mode 100644
index 0000000..9e4bd58
--- /dev/null
+++ b/src/3rd/kissfft/test/test_real.c
@@ -0,0 +1,179 @@
+/*
+ *  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_fftr.h"
+#include "_kiss_fft_guts.h"
+#include <sys/times.h>
+#include <time.h>
+#include <unistd.h>
+
+static double cputime(void)
+{
+    struct tms t;
+    times(&t);
+    return (double)(t.tms_utime + t.tms_stime)/  sysconf(_SC_CLK_TCK) ;
+}
+
+static
+kiss_fft_scalar rand_scalar(void) 
+{
+#ifdef USE_SIMD
+    return _mm_set1_ps(rand()-RAND_MAX/2);
+#else
+    kiss_fft_scalar s = (kiss_fft_scalar)(rand() -RAND_MAX/2);
+    return s/2;
+#endif
+}
+
+static
+double snr_compare( kiss_fft_cpx * vec1,kiss_fft_cpx * vec2, int n)
+{
+    int k;
+    double sigpow=1e-10,noisepow=1e-10,err,snr,scale=0;
+
+#ifdef USE_SIMD
+    float *fv1 = (float*)vec1;
+    float *fv2 = (float*)vec2;
+    for (k=0;k<8*n;++k) {
+        sigpow += *fv1 * *fv1;
+        err = *fv1 - *fv2;
+        noisepow += err*err;
+        ++fv1;
+        ++fv2;
+    }
+#else
+    for (k=0;k<n;++k) {
+        sigpow += (double)vec1[k].r * (double)vec1[k].r + 
+                  (double)vec1[k].i * (double)vec1[k].i;
+        err = (double)vec1[k].r - (double)vec2[k].r;
+        noisepow += err * err;
+        err = (double)vec1[k].i - (double)vec2[k].i;
+        noisepow += err * err;
+
+        if (vec1[k].r)
+            scale +=(double) vec2[k].r / (double)vec1[k].r;
+    }
+#endif    
+    snr = 10*log10( sigpow / noisepow );
+    scale /= n;
+    if (snr<10) {
+        printf( "\npoor snr, try a scaling factor %f\n" , scale );
+        exit(1);
+    }
+    return snr;
+}
+
+#ifndef NUMFFTS
+#define NUMFFTS 10000
+#endif
+
+
+int main(int argc,char ** argv)
+{
+    int nfft = 8*3*5;
+    double ts,tfft,trfft;
+    int i;
+    if (argc>1)
+        nfft = atoi(argv[1]);
+    kiss_fft_cpx cin[nfft];
+    kiss_fft_cpx cout[nfft];
+    kiss_fft_cpx sout[nfft];
+    kiss_fft_cfg  kiss_fft_state;
+    kiss_fftr_cfg  kiss_fftr_state;
+
+    kiss_fft_scalar rin[nfft+2];
+    kiss_fft_scalar rout[nfft+2];
+    kiss_fft_scalar zero;
+    memset(&zero,0,sizeof(zero) ); // ugly way of setting short,int,float,double, or __m128 to zero
+
+    srand(time(0));
+
+    for (i=0;i<nfft;++i) {
+        rin[i] = rand_scalar();
+        cin[i].r = rin[i];
+        cin[i].i = zero;
+    }
+
+    kiss_fft_state = kiss_fft_alloc(nfft,0,0,0);
+    kiss_fftr_state = kiss_fftr_alloc(nfft,0,0,0);
+    kiss_fft(kiss_fft_state,cin,cout);
+    kiss_fftr(kiss_fftr_state,rin,sout);
+    /*
+    printf(" results from kiss_fft : (%f,%f), (%f,%f), (%f,%f) ...\n "
+            , (float)cout[0].r , (float)cout[0].i
+            , (float)cout[1].r , (float)cout[1].i
+            , (float)cout[2].r , (float)cout[2].i); 
+    printf(" results from kiss_fftr: (%f,%f), (%f,%f), (%f,%f) ...\n "
+            , (float)sout[0].r , (float)sout[0].i
+            , (float)sout[1].r , (float)sout[1].i
+            , (float)sout[2].r , (float)sout[2].i); 
+    */
+        
+    printf( "nfft=%d, inverse=%d, snr=%g\n",
+            nfft,0, snr_compare(cout,sout,(nfft/2)+1) );
+    ts = cputime();
+    for (i=0;i<NUMFFTS;++i) {
+        kiss_fft(kiss_fft_state,cin,cout);
+    }
+    tfft = cputime() - ts;
+    
+    ts = cputime();
+    for (i=0;i<NUMFFTS;++i) {
+        kiss_fftr( kiss_fftr_state, rin, cout );
+        /* kiss_fftri(kiss_fftr_state,cout,rin); */
+    }
+    trfft = cputime() - ts;
+
+    printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft);
+
+    free(kiss_fft_state);
+    free(kiss_fftr_state);
+
+    kiss_fft_state = kiss_fft_alloc(nfft,1,0,0);
+    kiss_fftr_state = kiss_fftr_alloc(nfft,1,0,0);
+
+    memset(cin,0,sizeof(cin));
+#if 1
+    for (i=1;i< nfft/2;++i) {
+        //cin[i].r = (kiss_fft_scalar)(rand()-RAND_MAX/2);
+        cin[i].r = rand_scalar();
+        cin[i].i = rand_scalar();
+    }
+#else
+    cin[0].r = 12000;
+    cin[3].r = 12000;
+    cin[nfft/2].r = 12000;
+#endif
+
+    // conjugate symmetry of real signal 
+    for (i=1;i< nfft/2;++i) {
+        cin[nfft-i].r = cin[i].r;
+        cin[nfft-i].i = - cin[i].i;
+    }
+
+    kiss_fft(kiss_fft_state,cin,cout);
+    kiss_fftri(kiss_fftr_state,cin,rout);
+    /*
+    printf(" results from inverse kiss_fft : (%f,%f), (%f,%f), (%f,%f), (%f,%f), (%f,%f) ...\n "
+            , (float)cout[0].r , (float)cout[0].i , (float)cout[1].r , (float)cout[1].i , (float)cout[2].r , (float)cout[2].i , (float)cout[3].r , (float)cout[3].i , (float)cout[4].r , (float)cout[4].i
+            ); 
+
+    printf(" results from inverse kiss_fftr: %f,%f,%f,%f,%f ... \n"
+            ,(float)rout[0] ,(float)rout[1] ,(float)rout[2] ,(float)rout[3] ,(float)rout[4]);
+*/
+    for (i=0;i<nfft;++i) {
+        sout[i].r = rout[i];
+        sout[i].i = zero;
+    }
+
+    printf( "nfft=%d, inverse=%d, snr=%g\n",
+            nfft,1, snr_compare(cout,sout,nfft/2) );
+    free(kiss_fft_state);
+    free(kiss_fftr_state);
+
+    return 0;
+}
diff --git a/src/3rd/kissfft/test/test_vs_dft.c b/src/3rd/kissfft/test/test_vs_dft.c
new file mode 100644
index 0000000..9a44129
--- /dev/null
+++ b/src/3rd/kissfft/test/test_vs_dft.c
@@ -0,0 +1,81 @@
+/*
+ *  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.h"
+
+
+void check(kiss_fft_cpx  * in,kiss_fft_cpx  * out,int nfft,int isinverse)
+{
+    int bin,k;
+    double errpow=0,sigpow=0;
+    
+    for (bin=0;bin<nfft;++bin) {
+        double ansr = 0;
+        double ansi = 0;
+        double difr;
+        double difi;
+
+        for (k=0;k<nfft;++k) {
+            double phase = -2*M_PI*bin*k/nfft;
+            double re = cos(phase);
+            double im = sin(phase);
+            if (isinverse)
+                im = -im;
+
+#ifdef FIXED_POINT
+            re /= nfft;
+            im /= nfft;
+#endif            
+
+            ansr += in[k].r * re - in[k].i * im;
+            ansi += in[k].r * im + in[k].i * re;
+        }
+        difr = ansr - out[bin].r;
+        difi = ansi - out[bin].i;
+        errpow += difr*difr + difi*difi;
+        sigpow += ansr*ansr+ansi*ansi;
+    }
+    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,10*log10(sigpow/errpow) );
+}
+
+void test1d(int nfft,int isinverse)
+{
+    size_t buflen = sizeof(kiss_fft_cpx)*nfft;
+
+    kiss_fft_cpx  * in = (kiss_fft_cpx*)malloc(buflen);
+    kiss_fft_cpx  * out= (kiss_fft_cpx*)malloc(buflen);
+    kiss_fft_cfg  cfg = kiss_fft_alloc(nfft,isinverse,0,0);
+    int k;
+
+    for (k=0;k<nfft;++k) {
+        in[k].r = (rand() % 65536) - 32768;
+        in[k].i = (rand() % 65536) - 32768;
+    }
+
+    kiss_fft(cfg,in,out);
+
+    check(in,out,nfft,isinverse);
+
+    free(in);
+    free(out);
+    free(cfg);
+}
+
+int main(int argc,char ** argv)
+{
+    if (argc>1) {
+        int k;
+        for (k=1;k<argc;++k) {
+            test1d(atoi(argv[k]),0);
+            test1d(atoi(argv[k]),1);
+        }
+    }else{
+        test1d(32,0);
+        test1d(32,1);
+    }
+    return 0;
+}
diff --git a/src/3rd/kissfft/test/testcpp.cc b/src/3rd/kissfft/test/testcpp.cc
new file mode 100644
index 0000000..a62f6e0
--- /dev/null
+++ b/src/3rd/kissfft/test/testcpp.cc
@@ -0,0 +1,80 @@
+/*
+ *  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 "kissfft.hh"
+#include <iostream>
+#include <cstdlib>
+#include <typeinfo>
+
+#include <sys/time.h>
+static inline
+double curtime(void)
+{
+    struct timeval tv;
+    gettimeofday(&tv, NULL);
+    return (double)tv.tv_sec + (double)tv.tv_usec*.000001;
+}
+
+using namespace std;
+
+template <class T>
+void dotest(int nfft)
+{
+    typedef kissfft<T> FFT;
+    typedef std::complex<T> cpx_type;
+
+    cout << "type:" << typeid(T).name() << " nfft:" << nfft;
+
+    FFT fft(nfft,false);
+
+    vector<cpx_type> inbuf(nfft);
+    vector<cpx_type> outbuf(nfft);
+    for (int k=0;k<nfft;++k)
+        inbuf[k]= cpx_type( 
+                (T)(rand()/(double)RAND_MAX - .5),
+                (T)(rand()/(double)RAND_MAX - .5) );
+    fft.transform( &inbuf[0] , &outbuf[0] );
+
+    long double totalpower=0;
+    long double difpower=0;
+    for (int k0=0;k0<nfft;++k0) {
+        complex<long double> acc = 0;
+        long double phinc = 2*k0* M_PIl / nfft;
+        for (int k1=0;k1<nfft;++k1) {
+            complex<long double> x(inbuf[k1].real(),inbuf[k1].imag()); 
+            acc += x * exp( complex<long double>(0,-k1*phinc) );
+        }
+        totalpower += norm(acc);
+        complex<long double> x(outbuf[k0].real(),outbuf[k0].imag()); 
+        complex<long double> dif = acc - x;
+        difpower += norm(dif);
+    }
+    cout << " RMSE:" << sqrt(difpower/totalpower) << "\t";
+
+    double t0 = curtime();
+    int nits=20e6/nfft;
+    for (int k=0;k<nits;++k) {
+        fft.transform( &inbuf[0] , &outbuf[0] );
+    }
+    double t1 = curtime();
+    cout << " MSPS:" << ( (nits*nfft)*1e-6/ (t1-t0) ) << endl;
+}
+
+int main(int argc,char ** argv)
+{
+    if (argc>1) {
+        for (int k=1;k<argc;++k) {
+            int nfft = atoi(argv[k]);
+            dotest<float>(nfft); dotest<double>(nfft); dotest<long double>(nfft);
+        }
+    }else{
+        dotest<float>(32); dotest<double>(32); dotest<long double>(32);
+        dotest<float>(1024); dotest<double>(1024); dotest<long double>(1024);
+        dotest<float>(840); dotest<double>(840); dotest<long double>(840);
+    }
+    return 0;
+}
diff --git a/src/3rd/kissfft/test/testkiss.py b/src/3rd/kissfft/test/testkiss.py
new file mode 100755
index 0000000..98320d7
--- /dev/null
+++ b/src/3rd/kissfft/test/testkiss.py
@@ -0,0 +1,167 @@
+#!/usr/bin/env python
+#  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.
+
+import math
+import sys
+import os
+import random
+import struct
+import popen2
+import getopt
+import numpy
+
+pi=math.pi
+e=math.e
+j=complex(0,1)
+
+doreal=0
+
+datatype = os.environ.get('DATATYPE','float')
+
+util = '../tools/fft_' + datatype
+minsnr=90
+if datatype == 'double':
+    fmt='d'
+elif datatype=='int16_t':
+    fmt='h'
+    minsnr=10
+elif datatype=='int32_t':
+    fmt='i'
+elif datatype=='simd':
+    fmt='4f'
+    sys.stderr.write('testkiss.py does not yet test simd')
+    sys.exit(0)
+elif datatype=='float':
+    fmt='f'
+else:
+    sys.stderr.write('unrecognized datatype %s\n' % datatype)
+    sys.exit(1)
+ 
+
+def dopack(x,cpx=1):
+    x = numpy.reshape( x, ( numpy.size(x),) )
+    
+    if cpx:
+        s = ''.join( [ struct.pack(fmt*2,c.real,c.imag) for c in x ] )
+    else:
+        s = ''.join( [ struct.pack(fmt,c.real) for c in x ] )
+    return s
+
+def dounpack(x,cpx):
+    uf = fmt * ( len(x) / struct.calcsize(fmt) )
+    s = struct.unpack(uf,x)
+    if cpx:
+        return numpy.array(s[::2]) + numpy.array( s[1::2] )*j
+    else:
+        return numpy.array(s )
+
+def make_random(dims=[1]):
+    res = []
+    for i in range(dims[0]):
+        if len(dims)==1:
+            r=random.uniform(-1,1)
+            if doreal:
+                res.append( r )
+            else:
+                i=random.uniform(-1,1)
+                res.append( complex(r,i) )
+        else:
+            res.append( make_random( dims[1:] ) )
+    return numpy.array(res)
+
+def flatten(x):
+    ntotal = numpy.size(x)
+    return numpy.reshape(x,(ntotal,))
+
+def randmat( ndims ):
+    dims=[]
+    for i in range( ndims ):
+        curdim = int( random.uniform(2,5) )
+        if doreal and i==(ndims-1):
+            curdim = int(curdim/2)*2 # force even last dimension if real
+        dims.append( curdim )
+    return make_random(dims )
+
+def test_fft(ndims):
+    x=randmat( ndims )
+
+
+    if doreal:
+        xver = numpy.fft.rfftn(x)
+    else:
+        xver = numpy.fft.fftn(x)
+    
+    open('/tmp/fftexp.dat','w').write(dopack( flatten(xver) , True ) )
+
+    x2=dofft(x,doreal)
+    err = xver - x2
+    errf = flatten(err)
+    xverf = flatten(xver)
+    errpow = numpy.vdot(errf,errf)+1e-10
+    sigpow = numpy.vdot(xverf,xverf)+1e-10
+    snr = 10*math.log10(abs(sigpow/errpow) )
+    print 'SNR (compared to NumPy) : %.1fdB' % float(snr)
+
+    if snr<minsnr:
+        print 'xver=',xver
+        print 'x2=',x2
+        print 'err',err
+        sys.exit(1)
+ 
+def dofft(x,isreal):
+    dims=list( numpy.shape(x) )
+    x = flatten(x)
+
+    scale=1
+    if datatype=='int16_t':
+        x = 32767 * x
+        scale = len(x) / 32767.0
+    elif datatype=='int32_t':
+        x = 2147483647.0 * x
+        scale = len(x) / 2147483647.0
+
+    cmd='%s -n ' % util
+    cmd += ','.join([str(d) for d in dims])
+    if doreal:
+        cmd += ' -R '
+
+    print cmd
+    p = popen2.Popen3(cmd )
+
+    open('/tmp/fftin.dat','w').write(dopack( x , isreal==False ) )
+
+    p.tochild.write( dopack( x , isreal==False ) )
+    p.tochild.close()
+
+    res = dounpack( p.fromchild.read() , 1 )
+    open('/tmp/fftout.dat','w').write(dopack( flatten(res) , True ) )
+    if doreal:
+        dims[-1] = int( dims[-1]/2 ) + 1
+
+    res = scale * res
+
+    p.wait()
+    return numpy.reshape(res,dims)
+
+def main():
+    opts,args = getopt.getopt(sys.argv[1:],'r')
+    opts=dict(opts)
+
+    global doreal
+    doreal = opts.has_key('-r')
+
+    if doreal:
+        print 'Testing multi-dimensional real FFTs'
+    else:
+        print 'Testing multi-dimensional FFTs'
+
+    for dim in range(1,4):
+        test_fft( dim )
+
+if __name__ == "__main__":
+    main()
+
diff --git a/src/3rd/kissfft/test/twotonetest.c b/src/3rd/kissfft/test/twotonetest.c
new file mode 100644
index 0000000..5f08daf
--- /dev/null
+++ b/src/3rd/kissfft/test/twotonetest.c
@@ -0,0 +1,101 @@
+/*
+ *  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 <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include "kiss_fft.h"
+#include "kiss_fftr.h"
+#include <limits.h>
+
+
+static
+double two_tone_test( int nfft, int bin1,int bin2)
+{
+    kiss_fftr_cfg cfg	= NULL;
+    kiss_fft_cpx *kout	= NULL;
+    kiss_fft_scalar *tbuf = NULL;
+
+    int i;
+    double f1 = bin1*2*M_PI/nfft;
+    double f2 = bin2*2*M_PI/nfft;
+    double sigpow=0;
+    double noisepow=0;
+#if FIXED_POINT==32
+    long maxrange = LONG_MAX;
+#else
+    long maxrange = SHRT_MAX;/* works fine for float too*/
+#endif
+
+    cfg = kiss_fftr_alloc(nfft , 0, NULL, NULL);
+    tbuf	= KISS_FFT_MALLOC(nfft * sizeof(kiss_fft_scalar));
+    kout	= KISS_FFT_MALLOC(nfft * sizeof(kiss_fft_cpx));
+
+    /* generate a signal with two tones*/
+    for (i = 0; i < nfft; i++) {
+#ifdef USE_SIMD
+        tbuf[i] = _mm_set1_ps( (maxrange>>1)*cos(f1*i)
+                             + (maxrange>>1)*cos(f2*i) );
+#else        
+        tbuf[i] =  (maxrange>>1)*cos(f1*i)
+                 + (maxrange>>1)*cos(f2*i);
+#endif        
+    }
+
+    kiss_fftr(cfg, tbuf, kout);
+
+    for (i=0;i < (nfft/2+1);++i) {
+#ifdef USE_SIMD        
+        double tmpr = (double)*(float*)&kout[i].r / (double)maxrange;
+        double tmpi = (double)*(float*)&kout[i].i / (double)maxrange;
+#else
+        double tmpr = (double)kout[i].r / (double)maxrange;
+        double tmpi = (double)kout[i].i / (double)maxrange;
+#endif
+        double mag2 = tmpr*tmpr + tmpi*tmpi;
+        if (i!=0 && i!= nfft/2)
+            mag2 *= 2; /* all bins except DC and Nyquist have symmetric counterparts implied*/
+
+        /* if there is power in one of the expected bins, it is signal, otherwise noise*/
+        if ( i!=bin1 && i != bin2 ) 
+            noisepow += mag2;
+        else
+            sigpow += mag2;
+    }
+    kiss_fft_cleanup();
+    /*printf("TEST %d,%d,%d noise @ %fdB\n",nfft,bin1,bin2,10*log10(noisepow/sigpow +1e-30) );*/
+    return 10*log10(sigpow/(noisepow+1e-50) );
+}
+
+int main(int argc,char ** argv)
+{
+    int nfft = 4*2*2*3*5;
+    if (argc>1) nfft = atoi(argv[1]);
+
+    int i,j;
+    double minsnr = 500;
+    double maxsnr = -500;
+    double snr;
+    for (i=0;i<nfft/2;i+= (nfft>>4)+1) {
+        for (j=i;j<nfft/2;j+=(nfft>>4)+7) {
+            snr = two_tone_test(nfft,i,j);
+            if (snr<minsnr) {
+                minsnr=snr;
+            }
+            if (snr>maxsnr) {
+                maxsnr=snr;
+            }
+        }
+    }
+    snr = two_tone_test(nfft,nfft/2,nfft/2);
+    if (snr<minsnr) minsnr=snr;
+    if (snr>maxsnr) maxsnr=snr;
+
+    printf("TwoToneTest: snr ranges from %ddB to %ddB\n",(int)minsnr,(int)maxsnr);
+    printf("sizeof(kiss_fft_scalar) = %d\n",(int)sizeof(kiss_fft_scalar) );
+    return 0;
+}