Numerical derivatives

def ifft_dif_recursive(f):
    f = np.asarray(f)
    N = f.shape[0]
    W = np.exp(2. * np.pi * 1.0j / N)
    if N == 1:
        return f
    else:
        indrange = np.arange(N // 2)
        f1,f2 = f[:N//2], f[N//2:]
        f_even = f1 + f2
        f_odd = (f1 - f2) * (W**(indrange))
        F1 = ifft_dif_recursive(f_even)
        F2 = ifft_dif_recursive(f_odd)

        F = np.zeros(N, dtype=complex)
        indrange = np.arange(N)
        F[indrange[::2]] = F1
        F[indrange[1::2]] = F2
        return F

def gauss(x):
    return np.exp(-x**2)
def dgauss(x):
    return -2*x*np.exp(-x**2)

x = np.linspace(-3, 3, 32)
y = gauss(x)
dy = dgauss(x)
N = x.shape[0]
dx = x[1] - x[0]

print(N)
plt.plot(x, y, color='red', lw=3)
plt.plot(x, dy, color='blue', lw=3)
plt.xlabel('x')
plt.ylabel('y(x)')

F = fft_dif_recursive(y)
freq = (2 * np.pi / N/ dx) * np.concatenate((np.arange(N // 2), np.arange(-N // 2, 0, 1)))
F = 1.0j * F * freq
dy2 = ifft_dif_recursive(F) / (N)
plt.plot(x, dy2.real, 'o')

Comparison with dft2

Ntest = 10
timings3 = []
dimension = []
for i in range(1, Ntest):
    f = [random.gauss(0, 1) for n in range(2**i)]
    timing = %timeit -o -n 1 fft_dif_recursive(f)
    dimension.append(2**i)
    timings3.append(timing.average)
plt.subplot(121)    
plt.plot(np.array(dimension), np.array(timings3), 'o', markersize=12, color='red')
plt.plot(np.array(dimension), np.array(timings3), color = 'red')
plt.plot(np.array(dimension2), np.array(timings2), 'o', markersize=12, color='blue')
plt.plot(np.array(dimension2), np.array(timings2), color='blue')


plt.xlabel("$N$")
plt.ylabel("Timing / s")
plt.show()