From b24d80769b00987eba65997413282a5f613a7cd3 Mon Sep 17 00:00:00 2001 From: Mark Borgerding Date: Fri, 19 Jul 2019 10:43:49 -0400 Subject: [PATCH] removed unused and rotted code (closes #25) --- .hgignore | 10 --- test/compfft.py | 97 ---------------------- test/fastfir.py | 107 ------------------------ test/fft.py | 201 --------------------------------------------- test/tailscrap.m | 26 ------ test/test_vs_dft.c | 81 ------------------ 6 files changed, 522 deletions(-) delete mode 100644 .hgignore delete mode 100755 test/compfft.py delete mode 100755 test/fastfir.py delete mode 100755 test/fft.py delete mode 100644 test/tailscrap.m delete mode 100644 test/test_vs_dft.c diff --git a/.hgignore b/.hgignore deleted file mode 100644 index 978fa7f..0000000 --- a/.hgignore +++ /dev/null @@ -1,10 +0,0 @@ -syntax:glob -test/bm_* -test/st_* -test/tkfc_* -test/tr_* -tools/fastconv_* -tools/fastconvr_* -tools/fft_* -*.swp -*~ diff --git a/test/compfft.py b/test/compfft.py deleted file mode 100755 index d2671c1..0000000 --- a/test/compfft.py +++ /dev/null @@ -1,97 +0,0 @@ -#!/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/test/fastfir.py b/test/fastfir.py deleted file mode 100755 index 18662d4..0000000 --- a/test/fastfir.py +++ /dev/null @@ -1,107 +0,0 @@ -#!/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,'' - #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/test/fft.py b/test/fft.py deleted file mode 100755 index 4208a20..0000000 --- a/test/fft.py +++ /dev/null @@ -1,201 +0,0 @@ -#!/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/test/tailscrap.m b/test/tailscrap.m deleted file mode 100644 index abf9046..0000000 --- a/test/tailscrap.m +++ /dev/null @@ -1,26 +0,0 @@ -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/test/test_vs_dft.c b/test/test_vs_dft.c deleted file mode 100644 index 9a44129..0000000 --- a/test/test_vs_dft.c +++ /dev/null @@ -1,81 +0,0 @@ -/* - * 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;bin1) { - int k; - for (k=1;k