Decimation in frequency FFT (recursive implementation)

def fft_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 = fft_dif_recursive(f_even)
        F2 = fft_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

f = data_piano[:2000]
f = np.pad(f, 24)
F = fft_dif_recursive(f)
N = len(f)
dx = 1. / sampling_rate
freq = 1.0 / (N * dx) * np.arange(N)
plt.plot(freq[:100], F.real[:100])

F_numpy = np.fft.fft(f)
frequency_numpy = np.fft.fftfreq(N, dx)
print(np.allclose(F, F_numpy))
plt.plot(frequency_numpy[:100], F_numpy[:100].real)

Our implementation fft_dif_recursive as well as fft from Numpy give us the same result: