mirror of
https://github.com/mborgerding/kissfft.git
synced 2025-05-27 21:20:27 -04:00
fixed a bug with nfft==1
and added utility for more testing
This commit is contained in:
parent
95a7b856d1
commit
f4961ed74b
@ -18,6 +18,12 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
|
||||
fixed or floating point complex numbers. It also delares the kf_ internal functions.
|
||||
*/
|
||||
|
||||
/* the max length of the sequence of factors of nfft ,
|
||||
MAXFACTORS==32 this should be good up to more 1e15 samples
|
||||
( 2 * 3**31 )
|
||||
*/
|
||||
#define MAXFACTORS 32
|
||||
|
||||
static void kf_bfly2(
|
||||
kiss_fft_cpx * Fout,
|
||||
int fstride,
|
||||
@ -303,7 +309,7 @@ void * kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem )
|
||||
size_t memneeded = sizeof(kiss_fft_state)
|
||||
+ sizeof(kiss_fft_cpx)*nfft /* twiddle factors*/
|
||||
+ sizeof(kiss_fft_cpx)*nfft /* tmpbuf*/
|
||||
+ sizeof(int)*nfft /* factors*/
|
||||
+ sizeof(int)*2*MAXFACTORS /* factors*/
|
||||
+ sizeof(kiss_fft_cpx)*nfft; /* scratch*/
|
||||
|
||||
if ( lenmem==NULL ) {
|
||||
|
@ -13,6 +13,7 @@ BENCHKISS=bm_kiss_$(DATATYPE)
|
||||
BENCHFFTW=bm_fftw_$(DATATYPE)
|
||||
SELFTEST=st_$(DATATYPE)
|
||||
TESTREAL=tr_$(DATATYPE)
|
||||
FFTUTIL=kf_$(DATATYPE)
|
||||
|
||||
ifeq "$(DATATYPE)" "short"
|
||||
TYPEFLAGS=-DFIXED_POINT -Dkiss_fft_scalar=short
|
||||
@ -30,12 +31,15 @@ else
|
||||
endif
|
||||
|
||||
|
||||
all: $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL)
|
||||
all: $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL) $(FFTUTIL)
|
||||
|
||||
#CFLAGS=-Wall -O3 -ansi -pedantic -march=pentiumpro -ffast-math -fomit-frame-pointer
|
||||
# If the above flags do not work, try the following
|
||||
CFLAGS=-Wall -O3
|
||||
|
||||
$(FFTUTIL): ../kiss_fft.c fftutil.c kiss_fft2d.c kiss_fftr.c
|
||||
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
|
||||
|
||||
$(SELFTEST): ../kiss_fft.c $(SELFTESTSRC)
|
||||
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
|
||||
|
||||
@ -68,4 +72,4 @@ selftest_short.c:
|
||||
./mk_test.py -s 10 12 14 > selftest_short.c
|
||||
|
||||
clean:
|
||||
rm -f *~ bm_* st_* tr_*
|
||||
rm -f *~ bm_* st_* tr_* kf_*
|
||||
|
83
test/compfft.py
Executable file
83
test/compfft.py
Executable file
@ -0,0 +1,83 @@
|
||||
#!/usr/local/bin/python2.3
|
||||
|
||||
# use FFTPACK as a baseline
|
||||
import FFT
|
||||
from Numeric import *
|
||||
import math
|
||||
import random
|
||||
import sys
|
||||
import struct
|
||||
|
||||
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)
|
||||
|
||||
util = opts.get('-u','./kf_float')
|
||||
|
||||
try:
|
||||
n = int(opts['-n'])
|
||||
cpx = opts.get('-R') is None
|
||||
fmt=opts.get('-t','f')
|
||||
except KeyError:
|
||||
sys.stderr.write("""
|
||||
usage: compfft.py
|
||||
-n nfft
|
||||
-u utilname : see sample_code/fftutil.c
|
||||
-R : real-optimized version
|
||||
""")
|
||||
|
||||
x = randbuf(n,cpx)
|
||||
|
||||
cmd = '%s -n %d ' % ( util, n )
|
||||
if cpx:
|
||||
xout = FFT.fft(x)
|
||||
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 )
|
||||
|
||||
sigpow = sum( xout * conjugate(xout) )
|
||||
print len(xout)
|
||||
print len(xoutcomp)
|
||||
|
||||
diff = xout-xoutcomp
|
||||
noisepow = sum( diff * conjugate(diff) )
|
||||
|
||||
snr = 10 * math.log10(abs( sigpow / noisepow ) )
|
||||
print 'NFFT=%d,SNR = %f dB' % (n,snr)
|
||||
|
||||
def dopack(x,fmt,cpx):
|
||||
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()
|
@ -13,6 +13,7 @@ BENCHKISS=bm_kiss_$(DATATYPE)
|
||||
BENCHFFTW=bm_fftw_$(DATATYPE)
|
||||
SELFTEST=st_$(DATATYPE)
|
||||
TESTREAL=tr_$(DATATYPE)
|
||||
FFTUTIL=kf_$(DATATYPE)
|
||||
|
||||
ifeq "$(DATATYPE)" "short"
|
||||
TYPEFLAGS=-DFIXED_POINT -Dkiss_fft_scalar=short
|
||||
@ -30,12 +31,15 @@ else
|
||||
endif
|
||||
|
||||
|
||||
all: $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL)
|
||||
all: $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL) $(FFTUTIL)
|
||||
|
||||
#CFLAGS=-Wall -O3 -ansi -pedantic -march=pentiumpro -ffast-math -fomit-frame-pointer
|
||||
# If the above flags do not work, try the following
|
||||
CFLAGS=-Wall -O3
|
||||
|
||||
$(FFTUTIL): ../kiss_fft.c fftutil.c kiss_fft2d.c kiss_fftr.c
|
||||
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
|
||||
|
||||
$(SELFTEST): ../kiss_fft.c $(SELFTESTSRC)
|
||||
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
|
||||
|
||||
@ -68,4 +72,4 @@ selftest_short.c:
|
||||
./mk_test.py -s 10 12 14 > selftest_short.c
|
||||
|
||||
clean:
|
||||
rm -f *~ bm_* st_* tr_*
|
||||
rm -f *~ bm_* st_* tr_* kf_*
|
||||
|
@ -19,61 +19,90 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
|
||||
#include <unistd.h>
|
||||
|
||||
#include "kiss_fft.h"
|
||||
#include "kiss_fft2d.h"
|
||||
#include "kiss_fftr.h"
|
||||
|
||||
void fft_file(FILE * fin,FILE * fout,int nfft,int nrows,int isinverse,int useascii,int times)
|
||||
void fft_file(FILE * fin,FILE * fout,int nfft,int isinverse)
|
||||
{
|
||||
int i;
|
||||
void *st;
|
||||
kiss_fft_cpx * buf;
|
||||
kiss_fft_cpx * bufout;
|
||||
|
||||
|
||||
buf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft *nrows );
|
||||
bufout = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft *nrows);
|
||||
if (nrows!=1)
|
||||
st = kiss_fft2d_alloc( nrows,nfft ,isinverse ,0,0);
|
||||
else
|
||||
st = kiss_fft_alloc( nfft ,isinverse ,0,0);
|
||||
buf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft );
|
||||
bufout = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft );
|
||||
st = kiss_fft_alloc( nfft ,isinverse ,0,0);
|
||||
|
||||
while ( fread( buf , sizeof(kiss_fft_cpx) * nfft * nrows ,1, fin ) > 0 ) {
|
||||
for (i=0;i<times;++i)
|
||||
if (nrows!=1)
|
||||
kiss_fft2d( st , buf ,bufout);
|
||||
else
|
||||
kiss_fft( st , buf ,bufout);
|
||||
|
||||
if (useascii) {
|
||||
int i;
|
||||
for (i=0;i<nfft*nrows;++i)
|
||||
fprintf(fout, "(%g,%g) ", (double)bufout[i].r,(double)bufout[i].i);
|
||||
}else{
|
||||
fwrite( bufout , sizeof(kiss_fft_cpx) , nfft*nrows , fout );
|
||||
}
|
||||
while ( fread( buf , sizeof(kiss_fft_cpx) * nfft ,1, fin ) > 0 ) {
|
||||
kiss_fft( st , buf ,bufout);
|
||||
fwrite( bufout , sizeof(kiss_fft_cpx) , nfft , fout );
|
||||
}
|
||||
free(st);
|
||||
free(buf);
|
||||
free(bufout);
|
||||
}
|
||||
|
||||
void fft_file2d(FILE * fin,FILE * fout,int nfft,int nrows,int isinverse)
|
||||
{
|
||||
void *st;
|
||||
kiss_fft_cpx *buf;
|
||||
kiss_fft_cpx *bufout;
|
||||
|
||||
buf = (kiss_fft_cpx *) malloc (sizeof (kiss_fft_cpx) * nfft * nrows);
|
||||
bufout = (kiss_fft_cpx *) malloc (sizeof (kiss_fft_cpx) * nfft * nrows);
|
||||
st = kiss_fft2d_alloc (nrows, nfft, isinverse, 0, 0);
|
||||
|
||||
while (fread (buf, sizeof (kiss_fft_cpx) * nfft * nrows, 1, fin) > 0) {
|
||||
kiss_fft2d (st, buf, bufout);
|
||||
fwrite (bufout, sizeof (kiss_fft_cpx), nfft * nrows, fout);
|
||||
}
|
||||
free (st);
|
||||
free (buf);
|
||||
free (bufout);
|
||||
}
|
||||
|
||||
void fft_file_real(FILE * fin,FILE * fout,int nfft,int isinverse)
|
||||
{
|
||||
void *st;
|
||||
kiss_fft_scalar * rbuf;
|
||||
kiss_fft_cpx * cbuf;
|
||||
|
||||
rbuf = (kiss_fft_scalar*)malloc(sizeof(kiss_fft_scalar) * nfft );
|
||||
cbuf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * (nfft/2+1) );
|
||||
st = kiss_fftr_alloc( nfft ,isinverse ,0,0);
|
||||
|
||||
if (isinverse==0) {
|
||||
while ( fread( rbuf , sizeof(kiss_fft_scalar) * nfft ,1, fin ) > 0 ) {
|
||||
kiss_fftr( st , rbuf ,cbuf);
|
||||
fwrite( cbuf , sizeof(kiss_fft_cpx) , (nfft/2 + 1) , fout );
|
||||
}
|
||||
}else{
|
||||
while ( fread( cbuf , sizeof(kiss_fft_cpx) * (nfft/2+1) ,1, fin ) > 0 ) {
|
||||
kiss_fftri( st , cbuf ,rbuf);
|
||||
fwrite( rbuf , sizeof(kiss_fft_scalar) , nfft , fout );
|
||||
}
|
||||
}
|
||||
free(st);
|
||||
free(rbuf);
|
||||
free(cbuf);
|
||||
}
|
||||
|
||||
int main(int argc,char ** argv)
|
||||
{
|
||||
int nfft=1024;
|
||||
int isinverse=0;
|
||||
int isreal=0;
|
||||
FILE *fin=stdin;
|
||||
FILE *fout=stdout;
|
||||
int useascii=0;
|
||||
int times=1;
|
||||
int nrows=1;
|
||||
|
||||
while (1) {
|
||||
int c=getopt(argc,argv,"n:iax:r:");
|
||||
int c=getopt(argc,argv,"n:ir:R");
|
||||
if (c==-1) break;
|
||||
switch (c) {
|
||||
case 'a':useascii=1;break;
|
||||
case 'n':nfft = atoi(optarg);break;
|
||||
case 'r':nrows = atoi(optarg);break;
|
||||
case 'i':isinverse=1;break;
|
||||
case 'x':times=atoi(optarg);break;
|
||||
case 'R':isreal=1;break;
|
||||
}
|
||||
}
|
||||
|
||||
@ -89,7 +118,12 @@ int main(int argc,char ** argv)
|
||||
++optind;
|
||||
}
|
||||
|
||||
fft_file(fin,fout,nfft,nrows,isinverse,useascii,times);
|
||||
if (nrows>1 && !isreal)
|
||||
fft_file2d(fin,fout,nfft,nrows,isinverse);
|
||||
else if (nrows==1 && !isreal)
|
||||
fft_file(fin,fout,nfft,isinverse);
|
||||
else if (nrows==1 && isreal)
|
||||
fft_file_real(fin,fout,nfft,isinverse);
|
||||
|
||||
if (fout!=stdout) fclose(fout);
|
||||
if (fin!=stdin) fclose(fin);
|
||||
|
Loading…
Reference in New Issue
Block a user