mirror/kissfft
1
0
Fork 0
mirror of https://github.com/mborgerding/kissfft.git synced 2025-07-18 21:14:24 -04:00
Code Issues Packages Projects Releases Wiki Activity

initial commit of kiss_fft130.tar.gz contents

Browse Source
This commit is contained in:
Ian
2013-07-23 21:57:43 -04:00
commit 7d00183660
38 changed files with 4748 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

108
test/Makefile Normal file
View File

@ -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

94
test/benchfftw.c Normal file
View File

@ -0,0 +1,94 @@
#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

122
test/benchkiss.c Normal file
View File

@ -0,0 +1,122 @@
#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;
}

92
test/compfft.py Executable file
View File

@ -0,0 +1,92 @@
#!/usr/bin/env python
# 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()

129
test/doit.c Normal file
View File

@ -0,0 +1,129 @@
/* this program is in the public domain
A program that helps the authors of the fine fftw library benchmark kiss
*/
#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);
}

102
test/fastfir.py Executable file
View File

@ -0,0 +1,102 @@
#!/usr/bin/env python
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()

196
test/fft.py Executable file
View File

@ -0,0 +1,196 @@
#!/usr/bin/env python
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)

117
test/mk_test.py Executable file
View File

@ -0,0 +1,117 @@
#!/usr/bin/env python
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()

49
test/pstats.c Normal file
View File

@ -0,0 +1,49 @@
#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();
}

7
test/pstats.h Normal file
View File

@ -0,0 +1,7 @@
#ifndef PSTATS_H
#define PSTATS_H
void pstats_init(void);
void pstats_report(void);
#endif

26
test/tailscrap.m Normal file
View File

@ -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

172
test/test_real.c Normal file
View File

@ -0,0 +1,172 @@
#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;
}

74
test/test_vs_dft.c Normal file
View File

@ -0,0 +1,74 @@
#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;
}

73
test/testcpp.cc Normal file
View File

@ -0,0 +1,73 @@
#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;
}

162
test/testkiss.py Executable file
View File

@ -0,0 +1,162 @@
#!/usr/bin/env python
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()

94
test/twotonetest.c Normal file
View File

@ -0,0 +1,94 @@
#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;
}