import numpy as np
def p(q,r):
print "%s=" % q
print (r)
print
return None
def swap(x,i,j):
temp=x[i]
x[i]=x[j]
x[j]=temp
return None
def bitreversal(x):
N = int(x.size)
r = np.zeros(N, dtype=int)
M = 1
while M < N:
k = 0
while k < M:
T = int(2*r[k])
r[k] = T
r[k+M] = T+1
k = k+1
M = int(2*M)
n=1
while n < N-1:
if r[n]>n:
swap(x,n,r[n])
n=n+1
return x
def fft0 (z) :
Zo=bitreversal(z)
p=int(np.log2(int(Zo.size)))
x=Zo.astype(np.complex_)
twiddle=np.exp(-2j*np.pi*np.arange(0,0.5,1./x.size,dtype=np.complex_))
N=x.size
Bp=N/2;
Np=2;
twiddle_step_size=N/2
for P in range(0, p):
Npp= Np/2
baseT=0
for b in range(0, Bp):
baseB=baseT+Npp
for n in range(0, Npp):
twiddle_factor=n*twiddle_step_size
top= x[baseT+n]
bot=x[baseB+n]*twiddle[twiddle_factor]
x[baseT+n]=top+bot
x[baseB+n]=top-bot
baseT=baseT+Np
Bp=Bp/2
Np=Np*2
twiddle_step_size=twiddle_step_size/2
return x
def steps(g):
np.set_printoptions(precision=3,suppress=1)
N1=len(g)
N2=N1/2
x1=g[::2]
x2=g[1::2]
x2j=[1j*m for m in x2]
x=[a+b for a,b in zip(x1,x2j)]
x0=np.array(x)
x_=x0.astype(np.complex_)
X=fft0(x_)
Xc=np.conjugate(X)
X1=[0]*N2
X2=[0]*N2
for k in range(0,N2):
N2_k_mod=(N2-k)%N2
X1[k]=0.5*(X[k]+Xc[N2_k_mod])
X2[k]=-0.5*1j*(X[k]-Xc[N2_k_mod])
G=[0]*N1
Gn=[0]*N1
np.vectorize(complex)(G)
np.vectorize(complex)(X1)
np.vectorize(complex)(X2)
for k in range(0,N2):
G[k]=X1[k]+np.exp(-np.pi*1j*k/N2)*X2[k]
G=np.array(G)
G[N2]=0.5*((X[0]+Xc[0])+1j*(X[0]-Xc[0]))
for k in range (1,N2):
G[N1-k]= np.conjugate(G[k])
return G
def test():
flag=0
for r in range (3,8):
s=np.power(2,r)
g=np.random.rand(s)
G=steps(g)
Gpy=np.fft.fft(g)
p("G",G)
p("Gpy",Gpy)
tf=np.allclose(G,Gpy)
if tf==0: flag=1
assert(tf)
print tf
if flag==0: print ("All results were correct.")
print
test()