decimation in time split-radix FFT

sfft.py

import numpy as np
def sfft(x):
    n = len(x)
    if n==1:
        return x
    elif n==2:
        return np.r_[x[0]+x[1], x[0]-x[1]]
    else:
        w  = np.exp(-2j*np.pi/n*np.arange(n//4))
        w3 = np.exp(-2j*np.pi/n*3*np.arange(n//4))
        x02 = sfft(x[0::2])
        x1  = sfft(x[1::4])*w
        x3  = sfft(x[3::4])*w3
        x1_ = x1+x3
        x3_ = x1-x3
        return np.r_[x02[:n//4] + x1_,
                     x02[n//4:] - 1j*x3_,
                     x02[:n//4] - x1_,
                     x02[n//4:] + 1j*x3_]

if __name__=="__main__":
    print(sfft(np.arange(32)))