mirror/kissfft
1
0
Fork 0
mirror of https://github.com/mborgerding/kissfft.git synced 2026-04-08 05:32:42 -04:00
Code Issues Packages Projects Releases Wiki Activity

fixed scaling for short real

Browse Source 
moved fft.py from ./ to sample_code
This commit is contained in:
Mark Borgerding
2003-12-14 05:13:11 +00:00
parent 559c14b49b
commit 573536f48f
14 changed files with 96 additions and 244 deletions
Download Patch File Download Diff File Expand all files Collapse all files

55
test/Makefile
View File

@ -10,58 +10,61 @@ ifeq "$(DATATYPE)" ""
DATATYPE=float
endif
UTIL=fftutil_$(DATATYPE)
TESTREAL=tr_$(DATATYPE)
BENCHKISS=bm_kiss_$(DATATYPE)
BENCHFFTW=bm_fftw_$(DATATYPE)
all: $(UTIL) $(BENCHKISS) $(TESTREAL)
SELFTEST=st_$(DATATYPE)
TESTREAL=tr_$(DATATYPE)
ifeq "$(DATATYPE)" "short"
TYPEFLAGS=-DFIXED_POINT -Dkiss_fft_scalar=short
SELFTESTSRC=selftest_short.c
else
TYPEFLAGS=-Dkiss_fft_scalar=$(DATATYPE)
SELFTESTSRC=selftest.c
endif
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 -ansi -pedantic
all: $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL)
$(UTIL): ../kiss_fft.c ../kiss_fft2d.c fftutil.c
gcc -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) ../kiss_fft.c ../kiss_fft2d.c fftutil.c -lm
#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
$(SELFTEST): ../kiss_fft.c $(SELFTESTSRC)
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
$(TESTREAL): ../kiss_fft.c ../kiss_fftr.c test_real.c
gcc -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
$(BENCHKISS): benchkiss.c ../kiss_fft.c pstats.c
gcc -o $@ $(CFLAGS) -I.. benchkiss.c $(TYPEFLAGS) ../kiss_fft.c pstats.c -lm
fftw: $(BENCHFFTW)
./$(BENCHFFTW) -x $(NUMFFTS) -n $(NFFT)
$(CC) -o $@ $(CFLAGS) -I.. benchkiss.c $(TYPEFLAGS) ../kiss_fft.c pstats.c -lm
$(BENCHFFTW): benchfftw.c pstats.c
gcc -o $@ $(CFLAGS) -DDATATYPE$(DATATYPE) benchfftw.c pstats.c -lm -lfftw3f -lfftw3 -L /usr/local/lib/
time: all
@./$(BENCHKISS) -x $(NUMFFTS) -n $(NFFT)
@$(CC) -o $@ $(CFLAGS) -DDATATYPE$(DATATYPE) benchfftw.c pstats.c -lm -lfftw3f -lfftw3 -L /usr/local/lib/ || echo "FFTW not available for comparison"
POW2=256 512 1024 2048 4096 8192
POW3=243 729 2187
POW5=25 125 625 3125
mtime: all $(BENCHFFTW)
for n in $(POW2) $(POW3) $(POW5) ;do \
mtime: all
@for n in $(POW2) $(POW3) $(POW5) ;do \
echo ============================;\
./$(BENCHKISS) -x $(NUMFFTS) -n $$n;\
[ -x ./$(BENCHFFTW) ] && ./$(BENCHFFTW) -x $(NUMFFTS) -n $$n || true ; \
done
snr: all $(TESTREAL)
@echo "testkiss( $(NFFT) , 1, '$(DATATYPE)' );" | octave -q
@echo "testkiss( $(NFFT) , $(NROWS), '$(DATATYPE)' );" | octave -q
./$(TESTREAL)
test: all $(TESTREAL)
@echo "======SELF TEST"
@./$(SELFTEST)
@echo "======REAL FFT TEST"
@./$(TESTREAL)
@echo "======TIMING TEST"
@./$(BENCHKISS) -x $(NUMFFTS) -n $(NFFT)
@[ -x ./$(BENCHFFTW) ] && ./$(BENCHFFTW) -x $(NUMFFTS) -n $(NFFT) ||true
test: snr time fftw
selftest.c:
./mk_test.py 10 12 14 > selftest.c
selftest_short.c:
./mk_test.py -s 10 12 14 > selftest_short.c
clean:
rm -f *~ fftutil_* bm_* tr_* *.dat
rm -f *~ bm_* st_* tr_* *.dat

2
test/benchfftw.c
View File

@ -1,7 +1,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>
#include <getopt.h>
#include <unistd.h>
#include "pstats.h"
#ifdef DATATYPEdouble

2
test/benchkiss.c
View File

@ -1,6 +1,6 @@
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <unistd.h>
#include <sys/times.h>
#include <unistd.h>
#include "kiss_fft.h"

124
test/fft.py Executable file
View File

@ -0,0 +1,124 @@
#!/usr/local/bin/python2.3
import math
import sys
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]
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
if __name__ == "__main__":
main()

19
test/fileload.m
View File

@ -1,19 +0,0 @@
function data = fileload( fname , prec , iscomplex ,nrows )
f = fopen(fname,"r", "native");
data = fread(f,Inf,prec);
fclose(f);
len = length(data);
if iscomplex,
data = (data(1:2:len) + j*data(2:2:len) );
len = len/2;
endif
tdata = zeros(len/nrows,nrows);
tdata(:) = data;
data = tdata .';
endfunction

20
test/filesave.m
View File

@ -1,20 +0,0 @@
function filesave( fname , prec , data )
# kiss fft stores data by rows
data = vec(data.');
len = length(data);
f = fopen(fname,"w", "native");
if is_complex(data),
flat=zeros( 1,2*len);
flat(1:2:2*len) = real(data);
flat(2:2:2*len) = imag(data);
data = flat;
endif
fwrite(f,data,prec);
fclose(f);
endfunction

17
test/mk_test.py
View File

@ -23,6 +23,10 @@ def c_format(v,round=0):
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)
@ -44,12 +48,13 @@ def test_cpx( n,inverse ,short):
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);""" % inverse + """
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("FFT n=%d, inverse=%d, snr = %g dB\\n",NFFT,""" + str(inverse) + """,snr);
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
@ -58,6 +63,8 @@ def test_cpx( n,inverse ,short):
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;
@ -96,13 +103,13 @@ def main():
fftsizes = [ 1800 ]
print '#include "kiss_fft.h"'
print compare_func()
print "int main() {"
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 0;
return exit_code;
}
"""

30
test/test_real.c
View File

@ -1,5 +1,6 @@
#include "kiss_fft.h"
#include "_kiss_fft_guts.h"
#include <sys/times.h>
#include <time.h>
#include <unistd.h>
static double cputime()
@ -9,8 +10,6 @@ static double cputime()
return (double)(t.tms_utime + t.tms_stime)/ sysconf(_SC_CLK_TCK) ;
}
double snr_compare( kiss_fft_cpx * vec1,kiss_fft_cpx * vec2, int n)
{
int k;
@ -18,15 +17,19 @@ double snr_compare( kiss_fft_cpx * vec1,kiss_fft_cpx * vec2, int n)
sigpow = noisepow = .00000000000000000001;
for (k=0;k<n;++k) {
sigpow += vec1[k].r * vec1[k].r +
vec1[k].i * vec1[k].i;
err = vec1[k].r - vec2[k].r;
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 = vec1[k].i - vec2[k].i;
err = (double)vec1[k].i - (double)vec2[k].i;
noisepow += err * err;
if (vec1[k].r)
scale += vec2[k].r / vec1[k].r;
scale +=(double) vec2[k].r / (double)vec1[k].r;
/*
fprintf(stderr,"vec1=");pcpx(vec1+k);
fprintf(stderr,"vec2=");pcpx(vec2+k);
*/
}
snr = 10*log10( sigpow / noisepow );
scale /= n;
@ -38,17 +41,13 @@ double snr_compare( kiss_fft_cpx * vec1,kiss_fft_cpx * vec2, int n)
#ifndef RANDOM
#define NFFT 8
#else
#define NFFT 2*3*5*7*11
#define NFFT 8*3*5
#endif
#ifndef NUMFFTS
#define NUMFFTS 1000
#endif
void pcpx(const char * msg, kiss_fft_cpx * c)
{
printf("%s: %g + %gi\n",msg,c->r,c->i);
}
int main()
{
@ -62,6 +61,7 @@ int main()
void * kiss_fft_state;
void * kiss_fftr_state;
srand(time(0));
for (i=0;i<NFFT;++i) {
#ifdef RANDOM
@ -69,7 +69,6 @@ int main()
#endif
cin[i].r = sin[i];
cin[i].i = 0;
/* printf("in[%d]",i);pcpx("",cin+i); */
}
kiss_fft_state = kiss_fft_alloc(NFFT,0,0,0);
@ -106,9 +105,6 @@ int main()
for (i=0;i<NFFT;++i) {
sout[i].r = sin[i];
sout[i].i = 0;
/* printf("sin[%d] = %f\t",i,sin[i]);
printf("cin[%d]",i);pcpx("",cin+i);
printf("sout[%d]",i);pcpx("",sout+i); */
}
printf( "nfft=%d, inverse=%d, snr=%g\n",

58
test/testkiss.m
View File

@ -1,58 +0,0 @@
function snr= testkiss( nfft , nrows, prec )
printf('### testing SNR for %d x %d point %s FFTs\n' , nfft,nrows,prec);
if strcmp( prec ,'short')
scale_t2f=nfft*nrows;
scale_f2t=nfft*nrows;
else
scale_t2f=1;
scale_f2t=1/(nfft*nrows);
endif
kfft= sprintf('./fftutil_%s',prec);
sig=floor(32767*rand(nrows,nfft)) + j*floor(32767*rand(nrows,nfft));
filesave('time.dat',prec,sig);
if nrows > 1
cmd = sprintf('%s -r %d -n %d time.dat freq.dat',kfft,nrows,nfft);
else
cmd = sprintf('%s -n %d time.dat freq.dat',kfft,nfft);
endif
system(cmd);
Fsigcomp=fileload('freq.dat',prec,1,nrows ) * scale_t2f;
if nrows == 1,
Fsig=fft(sig);
else
Fsig=fft2(sig);
endif
diff = Fsig - Fsigcomp;
noise_pow = sum( conj(diff).*diff );
sig_pow = sum( conj(Fsig).*Fsig );
snr_t2f = 10*log10( sig_pow / noise_pow )
avg_scale = mean( abs(Fsig) ./ abs(Fsigcomp) );
var_scale = var( abs(Fsig) ./ abs(Fsigcomp) );
if nrows > 1
cmd = sprintf('%s -r %d -i -n %d freq.dat time2.dat',kfft,nrows,nfft);
else
cmd = sprintf('%s -i -n %d freq.dat time2.dat',kfft,nfft);
endif
system(cmd);
sigcomp=fileload('time2.dat',prec,1,nrows) * scale_f2t;
diff = sig - sigcomp;
noise_pow = sum( conj(diff).*diff );
sig_pow = sum( conj(sig).*sig );
snr_f2t = 10*log10( sig_pow / noise_pow )
avg_scale = mean( abs(sig) ./ abs(sigcomp) );
var_scale = var( abs(sig) ./ abs(sigcomp) );
snr=[snr_t2f snr_f2t];
endfunction

63
test/testvec.c
View File

File diff suppressed because one or more lines are too long