алгоритм Кули-Тьюки

Помогите пож. найти ошибку в коде уже пару дней мучаюсь, на бумаге пару раз решил пример при n=8, всё сходится а при n=16 и т.д. не хочет (
Ошибка где-то в индексах в циклах
(работает только при n=8 илм 4)

Spoiler

Code:

def fft(data: np.ndarray) -> np.ndarray:

    """захардкоденый пример c n=8"""
    #data = np.zeros(8)
    #for i in range(len(data)):
    #    data[i] = (i+1)/10

    fdata = np.asarray(data, dtype='complex128')
    n = fdata.size

    # check if input length is power of two
    if not n > 0 or (n & (n - 1)) != 0:
        raise ValueError

    # shuffle data
    X = bit_rev(fdata)

    j = 1
    m = 1
    k = 0
    i = 0
    while m < int(np.log2(n)+2):
        for k in range(2**(m-1)):
            if k >= n/2:
                break
            w = np.exp(-2*1j*np.pi*k/(2*m))
            for i in range(k, n, 2**m):
                j = i+m
                temp = X[i]
                X[i] = temp + X[j]*w
                X[j] = temp - X[j]*w
        m *= 2
    # normalize fft signal
    X = 1 / np.sqrt(n) * X

    return X

bit_rev:

Spoiler

Code:

def bit_rev(data: np.ndarray) -> np.ndarray:
    v = np.ndarray(len(data), dtype='complex128')
    n = len(data)-1
    #num_bits = np.log(n, 2) + 1

    for i in range(n+1):
        if i != 0:
            bin_rep = '{:0{width}b}'.format(i, width=n.bit_length())
            bin_rev = bin_rep[::-1]
            bin_rev_int = int(bin_rev, 2)
            v[i] = data[bin_rev_int]
        else:
            v[i] = data[i]
    return v

test:

Spoiler

def test_5_fft(self):
data = np.random.randn(2)
data1 = ifft(fft(data))
self.assertTrue(np.allclose(data, data1))
self.assertTrue(np.allclose(fft(data), np.fft.fft(data)/np.sqrt(data.size)))

Спасибо люди добрые

 

для дибага решения с бумаги

для любого n:

Spoiler

Code:

def fft(data: np.ndarray) -> np.ndarray:

    """doc downt there"""

    fdata = np.asarray(data, dtype='complex128')
    n = fdata.size

    # check if input length is power of two
    if not n > 0 or (n & (n - 1)) != 0:
        raise ValueError

    #first step of FFT: shuffle data
    X = bit_rev(fdata)

    #recursively merge transforms
    m = 0
    h = 0
    while m < int(np.log2(n)):
        for k in range(0, (2**m)):
            w = np.exp(-2*1j*np.pi*k/(2**(m+1)))
            for i in range(k, n):
                i = k+h*2**(m+1)
                j = i+(2**m)
                if not i < n:
                    break
                temp = X[i]
                X[i] = X[i] + X[j]*w
                X[j] = X[j]*w*(-1) + temp
                h += 1
            h = 0
        m += 1
    #normalize fft signal
    X = 1 / np.sqrt(n) * X

    return X