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

Discussion in 'Python' started by Inten, 23 Jan 2020.

  1. Inten

    Inten New Member

    Joined:
    5 May 2011
    Messages:
    13
    Likes Received:
    2
    Reputations:
    0
    Помогите пож. найти ошибку в коде уже пару дней мучаюсь, на бумаге пару раз решил пример при n=8, всё сходится а при n=16 и т.д. не хочет (
    Ошибка где-то в индексах в циклах
    (работает только при n=8 илм 4)


    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:

    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:

    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)))

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

    Inten New Member

    Joined:
    5 May 2011
    Messages:
    13
    Likes Received:
    2
    Reputations:
    0
    для дибага решения с бумаги

    для любого n:
    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
    
     
    #2 Inten, 23 Jan 2020
    Last edited: 24 Jan 2020
    people2people likes this.