this repo has no description
removed unused and rotted code (closes #25)
-522
-10
.hgignore
-10
.hgignore
-97
test/compfft.py
-97
test/compfft.py
···
1
-
#!/usr/bin/env python
2
-
# Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
3
-
# This file is part of KISS FFT - https://github.com/mborgerding/kissfft
4
-
#
5
-
# SPDX-License-Identifier: BSD-3-Clause
6
-
# See COPYING file for more information.
7
-
8
-
# use FFTPACK as a baseline
9
-
import FFT
10
-
from Numeric import *
11
-
import math
12
-
import random
13
-
import sys
14
-
import struct
15
-
import fft
16
-
17
-
pi=math.pi
18
-
e=math.e
19
-
j=complex(0,1)
20
-
lims=(-32768,32767)
21
-
22
-
def randbuf(n,cpx=1):
23
-
res = array( [ random.uniform( lims[0],lims[1] ) for i in range(n) ] )
24
-
if cpx:
25
-
res = res + j*randbuf(n,0)
26
-
return res
27
-
28
-
def main():
29
-
from getopt import getopt
30
-
import popen2
31
-
opts,args = getopt( sys.argv[1:],'u:n:Rt:' )
32
-
opts=dict(opts)
33
-
exitcode=0
34
-
35
-
util = opts.get('-u','./kf_float')
36
-
37
-
try:
38
-
dims = [ int(d) for d in opts['-n'].split(',')]
39
-
cpx = opts.get('-R') is None
40
-
fmt=opts.get('-t','f')
41
-
except KeyError:
42
-
sys.stderr.write("""
43
-
usage: compfft.py
44
-
-n d1[,d2,d3...] : FFT dimension(s)
45
-
-u utilname : see sample_code/fftutil.c, default = ./kf_float
46
-
-R : real-optimized version\n""")
47
-
sys.exit(1)
48
-
49
-
x = fft.make_random( dims )
50
-
51
-
cmd = '%s -n %s ' % ( util, ','.join([ str(d) for d in dims]) )
52
-
if cpx:
53
-
xout = FFT.fftnd(x)
54
-
xout = reshape(xout,(size(xout),))
55
-
else:
56
-
cmd += '-R '
57
-
xout = FFT.real_fft(x)
58
-
59
-
proc = popen2.Popen3( cmd , bufsize=len(x) )
60
-
61
-
proc.tochild.write( dopack( x , fmt ,cpx ) )
62
-
proc.tochild.close()
63
-
xoutcomp = dounpack( proc.fromchild.read( ) , fmt ,1 )
64
-
#xoutcomp = reshape( xoutcomp , dims )
65
-
66
-
sig = xout * conjugate(xout)
67
-
sigpow = sum( sig )
68
-
69
-
diff = xout-xoutcomp
70
-
noisepow = sum( diff * conjugate(diff) )
71
-
72
-
snr = 10 * math.log10(abs( sigpow / noisepow ) )
73
-
if snr<100:
74
-
print xout
75
-
print xoutcomp
76
-
exitcode=1
77
-
print 'NFFT=%s,SNR = %f dB' % (str(dims),snr)
78
-
sys.exit(exitcode)
79
-
80
-
def dopack(x,fmt,cpx):
81
-
x = reshape( x, ( size(x),) )
82
-
if cpx:
83
-
s = ''.join( [ struct.pack('ff',c.real,c.imag) for c in x ] )
84
-
else:
85
-
s = ''.join( [ struct.pack('f',c) for c in x ] )
86
-
return s
87
-
88
-
def dounpack(x,fmt,cpx):
89
-
uf = fmt * ( len(x) / 4 )
90
-
s = struct.unpack(uf,x)
91
-
if cpx:
92
-
return array(s[::2]) + array( s[1::2] )*j
93
-
else:
94
-
return array(s )
95
-
96
-
if __name__ == "__main__":
97
-
main()
-107
test/fastfir.py
-107
test/fastfir.py
···
1
-
#!/usr/bin/env python
2
-
# Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
3
-
# This file is part of KISS FFT - https://github.com/mborgerding/kissfft
4
-
#
5
-
# SPDX-License-Identifier: BSD-3-Clause
6
-
# See COPYING file for more information.
7
-
8
-
from Numeric import *
9
-
from FFT import *
10
-
11
-
def make_random(len):
12
-
import random
13
-
res=[]
14
-
for i in range(int(len)):
15
-
r=random.uniform(-1,1)
16
-
i=random.uniform(-1,1)
17
-
res.append( complex(r,i) )
18
-
return res
19
-
20
-
def slowfilter(sig,h):
21
-
translen = len(h)-1
22
-
return convolve(sig,h)[translen:-translen]
23
-
24
-
def nextpow2(x):
25
-
return 2 ** math.ceil(math.log(x)/math.log(2))
26
-
27
-
def fastfilter(sig,h,nfft=None):
28
-
if nfft is None:
29
-
nfft = int( nextpow2( 2*len(h) ) )
30
-
H = fft( h , nfft )
31
-
scraplen = len(h)-1
32
-
keeplen = nfft-scraplen
33
-
res=[]
34
-
isdone = 0
35
-
lastidx = nfft
36
-
idx0 = 0
37
-
while not isdone:
38
-
idx1 = idx0 + nfft
39
-
if idx1 >= len(sig):
40
-
idx1 = len(sig)
41
-
lastidx = idx1-idx0
42
-
if lastidx <= scraplen:
43
-
break
44
-
isdone = 1
45
-
Fss = fft(sig[idx0:idx1],nfft)
46
-
fm = Fss * H
47
-
m = inverse_fft(fm)
48
-
res.append( m[scraplen:lastidx] )
49
-
idx0 += keeplen
50
-
return concatenate( res )
51
-
52
-
def main():
53
-
import sys
54
-
from getopt import getopt
55
-
opts,args = getopt(sys.argv[1:],'rn:l:')
56
-
opts=dict(opts)
57
-
58
-
siglen = int(opts.get('-l',1e4 ) )
59
-
hlen =50
60
-
61
-
nfft = int(opts.get('-n',128) )
62
-
usereal = opts.has_key('-r')
63
-
64
-
print 'nfft=%d'%nfft
65
-
# make a signal
66
-
sig = make_random( siglen )
67
-
# make an impulse response
68
-
h = make_random( hlen )
69
-
#h=[1]*2+[0]*3
70
-
if usereal:
71
-
sig=[c.real for c in sig]
72
-
h=[c.real for c in h]
73
-
74
-
# perform MAC filtering
75
-
yslow = slowfilter(sig,h)
76
-
#print '<YSLOW>',yslow,'</YSLOW>'
77
-
#yfast = fastfilter(sig,h,nfft)
78
-
yfast = utilfastfilter(sig,h,nfft,usereal)
79
-
#print yfast
80
-
print 'len(yslow)=%d'%len(yslow)
81
-
print 'len(yfast)=%d'%len(yfast)
82
-
diff = yslow-yfast
83
-
snr = 10*log10( abs( vdot(yslow,yslow) / vdot(diff,diff) ) )
84
-
print 'snr=%s' % snr
85
-
if snr < 10.0:
86
-
print 'h=',h
87
-
print 'sig=',sig[:5],'...'
88
-
print 'yslow=',yslow[:5],'...'
89
-
print 'yfast=',yfast[:5],'...'
90
-
91
-
def utilfastfilter(sig,h,nfft,usereal):
92
-
import compfft
93
-
import os
94
-
open( 'sig.dat','w').write( compfft.dopack(sig,'f',not usereal) )
95
-
open( 'h.dat','w').write( compfft.dopack(h,'f',not usereal) )
96
-
if usereal:
97
-
util = './fastconvr'
98
-
else:
99
-
util = './fastconv'
100
-
cmd = 'time %s -n %d -i sig.dat -h h.dat -o out.dat' % (util, nfft)
101
-
print cmd
102
-
ec = os.system(cmd)
103
-
print 'exited->',ec
104
-
return compfft.dounpack(open('out.dat').read(),'f',not usereal)
105
-
106
-
if __name__ == "__main__":
107
-
main()
-201
test/fft.py
-201
test/fft.py
···
1
-
#!/usr/bin/env python
2
-
# Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
3
-
# This file is part of KISS FFT - https://github.com/mborgerding/kissfft
4
-
#
5
-
# SPDX-License-Identifier: BSD-3-Clause
6
-
# See COPYING file for more information.
7
-
8
-
import math
9
-
import sys
10
-
import random
11
-
12
-
pi=math.pi
13
-
e=math.e
14
-
j=complex(0,1)
15
-
16
-
def fft(f,inv):
17
-
n=len(f)
18
-
if n==1:
19
-
return f
20
-
21
-
for p in 2,3,5:
22
-
if n%p==0:
23
-
break
24
-
else:
25
-
raise Exception('%s not factorable ' % n)
26
-
27
-
m = n/p
28
-
Fout=[]
29
-
for q in range(p): # 0,1
30
-
fp = f[q::p] # every p'th time sample
31
-
Fp = fft( fp ,inv)
32
-
Fout.extend( Fp )
33
-
34
-
for u in range(m):
35
-
scratch = Fout[u::m] # u to end in strides of m
36
-
for q1 in range(p):
37
-
k = q1*m + u # indices to Fout above that became scratch
38
-
Fout[ k ] = scratch[0] # cuz e**0==1 in loop below
39
-
for q in range(1,p):
40
-
if inv:
41
-
t = e ** ( j*2*pi*k*q/n )
42
-
else:
43
-
t = e ** ( -j*2*pi*k*q/n )
44
-
Fout[ k ] += scratch[q] * t
45
-
46
-
return Fout
47
-
48
-
def rifft(F):
49
-
N = len(F) - 1
50
-
Z = [0] * (N)
51
-
for k in range(N):
52
-
Fek = ( F[k] + F[-k-1].conjugate() )
53
-
Fok = ( F[k] - F[-k-1].conjugate() ) * e ** (j*pi*k/N)
54
-
Z[k] = Fek + j*Fok
55
-
56
-
fp = fft(Z , 1)
57
-
58
-
f = []
59
-
for c in fp:
60
-
f.append(c.real)
61
-
f.append(c.imag)
62
-
return f
63
-
64
-
def real_fft( f,inv ):
65
-
if inv:
66
-
return rifft(f)
67
-
68
-
N = len(f) / 2
69
-
70
-
res = f[::2]
71
-
ims = f[1::2]
72
-
73
-
fp = [ complex(r,i) for r,i in zip(res,ims) ]
74
-
print 'fft input ', fp
75
-
Fp = fft( fp ,0 )
76
-
print 'fft output ', Fp
77
-
78
-
F = [ complex(0,0) ] * ( N+1 )
79
-
80
-
F[0] = complex( Fp[0].real + Fp[0].imag , 0 )
81
-
82
-
for k in range(1,N/2+1):
83
-
tw = e ** ( -j*pi*(.5+float(k)/N ) )
84
-
85
-
F1k = Fp[k] + Fp[N-k].conjugate()
86
-
F2k = Fp[k] - Fp[N-k].conjugate()
87
-
F2k *= tw
88
-
F[k] = ( F1k + F2k ) * .5
89
-
F[N-k] = ( F1k - F2k ).conjugate() * .5
90
-
#F[N-k] = ( F1kp + e ** ( -j*pi*(.5+float(N-k)/N ) ) * F2kp ) * .5
91
-
#F[N-k] = ( F1k.conjugate() - tw.conjugate() * F2k.conjugate() ) * .5
92
-
93
-
F[N] = complex( Fp[0].real - Fp[0].imag , 0 )
94
-
return F
95
-
96
-
def main():
97
-
#fft_func = fft
98
-
fft_func = real_fft
99
-
100
-
tvec = [0.309655,0.815653,0.768570,0.591841,0.404767,0.637617,0.007803,0.012665]
101
-
Ftvec = [ complex(r,i) for r,i in zip(
102
-
[3.548571,-0.378761,-0.061950,0.188537,-0.566981,0.188537,-0.061950,-0.378761],
103
-
[0.000000,-1.296198,-0.848764,0.225337,0.000000,-0.225337,0.848764,1.296198] ) ]
104
-
105
-
F = fft_func( tvec,0 )
106
-
107
-
nerrs= 0
108
-
for i in range(len(Ftvec)/2 + 1):
109
-
if abs( F[i] - Ftvec[i] )> 1e-5:
110
-
print 'F[%d]: %s != %s' % (i,F[i],Ftvec[i])
111
-
nerrs += 1
112
-
113
-
print '%d errors in forward fft' % nerrs
114
-
if nerrs:
115
-
return
116
-
117
-
trec = fft_func( F , 1 )
118
-
119
-
for i in range(len(trec) ):
120
-
trec[i] /= len(trec)
121
-
122
-
for i in range(len(tvec) ):
123
-
if abs( trec[i] - tvec[i] )> 1e-5:
124
-
print 't[%d]: %s != %s' % (i,tvec[i],trec[i])
125
-
nerrs += 1
126
-
127
-
print '%d errors in reverse fft' % nerrs
128
-
129
-
130
-
def make_random(dims=[1]):
131
-
import Numeric
132
-
res = []
133
-
for i in range(dims[0]):
134
-
if len(dims)==1:
135
-
r=random.uniform(-1,1)
136
-
i=random.uniform(-1,1)
137
-
res.append( complex(r,i) )
138
-
else:
139
-
res.append( make_random( dims[1:] ) )
140
-
return Numeric.array(res)
141
-
142
-
def flatten(x):
143
-
import Numeric
144
-
ntotal = Numeric.product(Numeric.shape(x))
145
-
return Numeric.reshape(x,(ntotal,))
146
-
147
-
def randmat( ndims ):
148
-
dims=[]
149
-
for i in range( ndims ):
150
-
curdim = int( random.uniform(2,4) )
151
-
dims.append( curdim )
152
-
return make_random(dims )
153
-
154
-
def test_fftnd(ndims=3):
155
-
import FFT
156
-
import Numeric
157
-
158
-
x=randmat( ndims )
159
-
print 'dimensions=%s' % str( Numeric.shape(x) )
160
-
#print 'x=%s' %str(x)
161
-
xver = FFT.fftnd(x)
162
-
x2=myfftnd(x)
163
-
err = xver - x2
164
-
errf = flatten(err)
165
-
xverf = flatten(xver)
166
-
errpow = Numeric.vdot(errf,errf)+1e-10
167
-
sigpow = Numeric.vdot(xverf,xverf)+1e-10
168
-
snr = 10*math.log10(abs(sigpow/errpow) )
169
-
if snr<80:
170
-
print xver
171
-
print x2
172
-
print 'SNR=%sdB' % str( snr )
173
-
174
-
def myfftnd(x):
175
-
import Numeric
176
-
xf = flatten(x)
177
-
Xf = fftndwork( xf , Numeric.shape(x) )
178
-
return Numeric.reshape(Xf,Numeric.shape(x) )
179
-
180
-
def fftndwork(x,dims):
181
-
import Numeric
182
-
dimprod=Numeric.product( dims )
183
-
184
-
for k in range( len(dims) ):
185
-
cur_dim=dims[ k ]
186
-
stride=dimprod/cur_dim
187
-
next_x = [complex(0,0)]*len(x)
188
-
for i in range(stride):
189
-
next_x[i*cur_dim:(i+1)*cur_dim] = fft(x[i:(i+cur_dim)*stride:stride],0)
190
-
x = next_x
191
-
return x
192
-
193
-
if __name__ == "__main__":
194
-
try:
195
-
nd = int(sys.argv[1])
196
-
except:
197
-
nd=None
198
-
if nd:
199
-
test_fftnd( nd )
200
-
else:
201
-
sys.exit(0)
-26
test/tailscrap.m
-26
test/tailscrap.m
···
1
-
function maxabsdiff=tailscrap()
2
-
% test code for circular convolution with the scrapped portion
3
-
% at the tail of the buffer, rather than the front
4
-
%
5
-
% The idea is to rotate the zero-padded h (impulse response) buffer
6
-
% to the left nh-1 samples, rotating the junk samples as well.
7
-
% This could be very handy in avoiding buffer copies during fast filtering.
8
-
nh=10;
9
-
nfft=256;
10
-
11
-
h=rand(1,nh);
12
-
x=rand(1,nfft);
13
-
14
-
hpad=[ h(nh) zeros(1,nfft-nh) h(1:nh-1) ];
15
-
16
-
% baseline comparison
17
-
y1 = filter(h,1,x);
18
-
y1_notrans = y1(nh:nfft);
19
-
20
-
% fast convolution
21
-
y2 = ifft( fft(hpad) .* fft(x) );
22
-
y2_notrans=y2(1:nfft-nh+1);
23
-
24
-
maxabsdiff = max(abs(y2_notrans - y1_notrans))
25
-
26
-
end
-81
test/test_vs_dft.c
-81
test/test_vs_dft.c
···
1
-
/*
2
-
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
3
-
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
4
-
*
5
-
* SPDX-License-Identifier: BSD-3-Clause
6
-
* See COPYING file for more information.
7
-
*/
8
-
#include "kiss_fft.h"
9
-
10
-
11
-
void check(kiss_fft_cpx * in,kiss_fft_cpx * out,int nfft,int isinverse)
12
-
{
13
-
int bin,k;
14
-
double errpow=0,sigpow=0;
15
-
16
-
for (bin=0;bin<nfft;++bin) {
17
-
double ansr = 0;
18
-
double ansi = 0;
19
-
double difr;
20
-
double difi;
21
-
22
-
for (k=0;k<nfft;++k) {
23
-
double phase = -2*M_PI*bin*k/nfft;
24
-
double re = cos(phase);
25
-
double im = sin(phase);
26
-
if (isinverse)
27
-
im = -im;
28
-
29
-
#ifdef FIXED_POINT
30
-
re /= nfft;
31
-
im /= nfft;
32
-
#endif
33
-
34
-
ansr += in[k].r * re - in[k].i * im;
35
-
ansi += in[k].r * im + in[k].i * re;
36
-
}
37
-
difr = ansr - out[bin].r;
38
-
difi = ansi - out[bin].i;
39
-
errpow += difr*difr + difi*difi;
40
-
sigpow += ansr*ansr+ansi*ansi;
41
-
}
42
-
printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,10*log10(sigpow/errpow) );
43
-
}
44
-
45
-
void test1d(int nfft,int isinverse)
46
-
{
47
-
size_t buflen = sizeof(kiss_fft_cpx)*nfft;
48
-
49
-
kiss_fft_cpx * in = (kiss_fft_cpx*)malloc(buflen);
50
-
kiss_fft_cpx * out= (kiss_fft_cpx*)malloc(buflen);
51
-
kiss_fft_cfg cfg = kiss_fft_alloc(nfft,isinverse,0,0);
52
-
int k;
53
-
54
-
for (k=0;k<nfft;++k) {
55
-
in[k].r = (rand() % 65536) - 32768;
56
-
in[k].i = (rand() % 65536) - 32768;
57
-
}
58
-
59
-
kiss_fft(cfg,in,out);
60
-
61
-
check(in,out,nfft,isinverse);
62
-
63
-
free(in);
64
-
free(out);
65
-
free(cfg);
66
-
}
67
-
68
-
int main(int argc,char ** argv)
69
-
{
70
-
if (argc>1) {
71
-
int k;
72
-
for (k=1;k<argc;++k) {
73
-
test1d(atoi(argv[k]),0);
74
-
test1d(atoi(argv[k]),1);
75
-
}
76
-
}else{
77
-
test1d(32,0);
78
-
test1d(32,1);
79
-
}
80
-
return 0;
81
-
}