Chris@366: #!/usr/bin/env python Chris@366: Chris@366: import math Chris@366: import sys Chris@366: import random Chris@366: Chris@366: pi=math.pi Chris@366: e=math.e Chris@366: j=complex(0,1) Chris@366: Chris@366: def fft(f,inv): Chris@366: n=len(f) Chris@366: if n==1: Chris@366: return f Chris@366: Chris@366: for p in 2,3,5: Chris@366: if n%p==0: Chris@366: break Chris@366: else: Chris@366: raise Exception('%s not factorable ' % n) Chris@366: Chris@366: m = n/p Chris@366: Fout=[] Chris@366: for q in range(p): # 0,1 Chris@366: fp = f[q::p] # every p'th time sample Chris@366: Fp = fft( fp ,inv) Chris@366: Fout.extend( Fp ) Chris@366: Chris@366: for u in range(m): Chris@366: scratch = Fout[u::m] # u to end in strides of m Chris@366: for q1 in range(p): Chris@366: k = q1*m + u # indices to Fout above that became scratch Chris@366: Fout[ k ] = scratch[0] # cuz e**0==1 in loop below Chris@366: for q in range(1,p): Chris@366: if inv: Chris@366: t = e ** ( j*2*pi*k*q/n ) Chris@366: else: Chris@366: t = e ** ( -j*2*pi*k*q/n ) Chris@366: Fout[ k ] += scratch[q] * t Chris@366: Chris@366: return Fout Chris@366: Chris@366: def rifft(F): Chris@366: N = len(F) - 1 Chris@366: Z = [0] * (N) Chris@366: for k in range(N): Chris@366: Fek = ( F[k] + F[-k-1].conjugate() ) Chris@366: Fok = ( F[k] - F[-k-1].conjugate() ) * e ** (j*pi*k/N) Chris@366: Z[k] = Fek + j*Fok Chris@366: Chris@366: fp = fft(Z , 1) Chris@366: Chris@366: f = [] Chris@366: for c in fp: Chris@366: f.append(c.real) Chris@366: f.append(c.imag) Chris@366: return f Chris@366: Chris@366: def real_fft( f,inv ): Chris@366: if inv: Chris@366: return rifft(f) Chris@366: Chris@366: N = len(f) / 2 Chris@366: Chris@366: res = f[::2] Chris@366: ims = f[1::2] Chris@366: Chris@366: fp = [ complex(r,i) for r,i in zip(res,ims) ] Chris@366: print 'fft input ', fp Chris@366: Fp = fft( fp ,0 ) Chris@366: print 'fft output ', Fp Chris@366: Chris@366: F = [ complex(0,0) ] * ( N+1 ) Chris@366: Chris@366: F[0] = complex( Fp[0].real + Fp[0].imag , 0 ) Chris@366: Chris@366: for k in range(1,N/2+1): Chris@366: tw = e ** ( -j*pi*(.5+float(k)/N ) ) Chris@366: Chris@366: F1k = Fp[k] + Fp[N-k].conjugate() Chris@366: F2k = Fp[k] - Fp[N-k].conjugate() Chris@366: F2k *= tw Chris@366: F[k] = ( F1k + F2k ) * .5 Chris@366: F[N-k] = ( F1k - F2k ).conjugate() * .5 Chris@366: #F[N-k] = ( F1kp + e ** ( -j*pi*(.5+float(N-k)/N ) ) * F2kp ) * .5 Chris@366: #F[N-k] = ( F1k.conjugate() - tw.conjugate() * F2k.conjugate() ) * .5 Chris@366: Chris@366: F[N] = complex( Fp[0].real - Fp[0].imag , 0 ) Chris@366: return F Chris@366: Chris@366: def main(): Chris@366: #fft_func = fft Chris@366: fft_func = real_fft Chris@366: Chris@366: tvec = [0.309655,0.815653,0.768570,0.591841,0.404767,0.637617,0.007803,0.012665] Chris@366: Ftvec = [ complex(r,i) for r,i in zip( Chris@366: [3.548571,-0.378761,-0.061950,0.188537,-0.566981,0.188537,-0.061950,-0.378761], Chris@366: [0.000000,-1.296198,-0.848764,0.225337,0.000000,-0.225337,0.848764,1.296198] ) ] Chris@366: Chris@366: F = fft_func( tvec,0 ) Chris@366: Chris@366: nerrs= 0 Chris@366: for i in range(len(Ftvec)/2 + 1): Chris@366: if abs( F[i] - Ftvec[i] )> 1e-5: Chris@366: print 'F[%d]: %s != %s' % (i,F[i],Ftvec[i]) Chris@366: nerrs += 1 Chris@366: Chris@366: print '%d errors in forward fft' % nerrs Chris@366: if nerrs: Chris@366: return Chris@366: Chris@366: trec = fft_func( F , 1 ) Chris@366: Chris@366: for i in range(len(trec) ): Chris@366: trec[i] /= len(trec) Chris@366: Chris@366: for i in range(len(tvec) ): Chris@366: if abs( trec[i] - tvec[i] )> 1e-5: Chris@366: print 't[%d]: %s != %s' % (i,tvec[i],trec[i]) Chris@366: nerrs += 1 Chris@366: Chris@366: print '%d errors in reverse fft' % nerrs Chris@366: Chris@366: Chris@366: def make_random(dims=[1]): Chris@366: import Numeric Chris@366: res = [] Chris@366: for i in range(dims[0]): Chris@366: if len(dims)==1: Chris@366: r=random.uniform(-1,1) Chris@366: i=random.uniform(-1,1) Chris@366: res.append( complex(r,i) ) Chris@366: else: Chris@366: res.append( make_random( dims[1:] ) ) Chris@366: return Numeric.array(res) Chris@366: Chris@366: def flatten(x): Chris@366: import Numeric Chris@366: ntotal = Numeric.product(Numeric.shape(x)) Chris@366: return Numeric.reshape(x,(ntotal,)) Chris@366: Chris@366: def randmat( ndims ): Chris@366: dims=[] Chris@366: for i in range( ndims ): Chris@366: curdim = int( random.uniform(2,4) ) Chris@366: dims.append( curdim ) Chris@366: return make_random(dims ) Chris@366: Chris@366: def test_fftnd(ndims=3): Chris@366: import FFT Chris@366: import Numeric Chris@366: Chris@366: x=randmat( ndims ) Chris@366: print 'dimensions=%s' % str( Numeric.shape(x) ) Chris@366: #print 'x=%s' %str(x) Chris@366: xver = FFT.fftnd(x) Chris@366: x2=myfftnd(x) Chris@366: err = xver - x2 Chris@366: errf = flatten(err) Chris@366: xverf = flatten(xver) Chris@366: errpow = Numeric.vdot(errf,errf)+1e-10 Chris@366: sigpow = Numeric.vdot(xverf,xverf)+1e-10 Chris@366: snr = 10*math.log10(abs(sigpow/errpow) ) Chris@366: if snr<80: Chris@366: print xver Chris@366: print x2 Chris@366: print 'SNR=%sdB' % str( snr ) Chris@366: Chris@366: def myfftnd(x): Chris@366: import Numeric Chris@366: xf = flatten(x) Chris@366: Xf = fftndwork( xf , Numeric.shape(x) ) Chris@366: return Numeric.reshape(Xf,Numeric.shape(x) ) Chris@366: Chris@366: def fftndwork(x,dims): Chris@366: import Numeric Chris@366: dimprod=Numeric.product( dims ) Chris@366: Chris@366: for k in range( len(dims) ): Chris@366: cur_dim=dims[ k ] Chris@366: stride=dimprod/cur_dim Chris@366: next_x = [complex(0,0)]*len(x) Chris@366: for i in range(stride): Chris@366: next_x[i*cur_dim:(i+1)*cur_dim] = fft(x[i:(i+cur_dim)*stride:stride],0) Chris@366: x = next_x Chris@366: return x Chris@366: Chris@366: if __name__ == "__main__": Chris@366: try: Chris@366: nd = int(sys.argv[1]) Chris@366: except: Chris@366: nd=None Chris@366: if nd: Chris@366: test_fftnd( nd ) Chris@366: else: Chris@366: sys.exit(0)