47241


pythonでFFT（高速フーリエ変換）を実装しようと思っています  
コードはご覧の通りです

(FFT_sort.py)

[code]

    import numpy as np
    
    def sort(N):
    
        flag = ~(N & (N - 1))
        if flag != -1:
            return None
    
        result = np.zeros(N, dtype=np.int64)
    
        result[0] = 0
        result[1] = N / 2
    
        result[N / 2] = 1
        result[N / 2 + 1] = N / 2 + 1
    
        T = n  = N / 2
        oldN = n
        oldCat = 0
        cat = 1
    
        count = 1
        judge = 0
        next = 2
    
        linspace = np.array([i for i in range(N)])
    
        while(n != 2):
    
            n /= 2
            cat *= 2
    
            judge += cat
    
            normal = 2 ** (count + 1))
            for i in range(cat / 2):
    
                j = 1 + i * 2
                index = j * N / normal
    
                result[index] = next
                result[index + 1] = next + T
    
                result[index + T] = next + 1
                result[index + T + 1] = next + T + 1
    
                next += 2
    
            count += 1
            oldN = n 
    
        return result
    
[/code]

(FFT.py)

[code]

    import FFT_sort as fs
    import numpy as np
    from numpy import pi
    
    import matplotlib.pyplot as plt
    
    N = input()
    
    x = np.linspace(0, 2 * np.pi, N)
    y = np.sin(x)
    
    data = fs.sort(N)
    
    if data is None:
        print("error")
        quit()
    else:
        print(data)
    
    result = np.zeros(N, dtype=np.complex)
    result[:] = y[:]
    
    num = 1
    
    while(num != N):
    
        count = 0
    
        T = num
        F = N / (num * 2)
    
        print("F:", F)
    
        for i in range(N):
    
            judge = i / num
    
            if judge % 2 == 0:
                k = i + T
    
            else:
                k = i - T
    
            iF = i * F % N
    
            w = np.exp(float(i) * iF * -1j * 2 * pi / float(N))
    
            print("iF:", iF)
            print("w:", w)
    
            n = data[k]
    
            result[i] += w * result[n]
    
        count += 1
        num *= 2
    
    print(result)
    answer = np.fft.fft(y)
    print(answer)
    
[/code]

コードの入力はFFTのサンプリングの数で  
アウトプットは  
私のFFTの結果とnumpyでやった正規のFFTの結果です

ここで質問ですが  
私のコードでやったFFTと  
本家numpyのFFTの結果が全然違います

[code]

    [ 6.23785596 -4.42077404j  5.0713722  -3.1314097 j
     -3.8474327  -1.99237623j  2.05949337 -0.61485772j
      5.84009223 -4.84151225j -2.02845113 -0.82589246j
     -1.36232081 -8.6730753 j  1.26764366 -7.27397125j
     -0.01871131 -4.59613684j  1.0188491  -4.96565877j
      3.58757676-10.1782615 j  2.92161011 -6.39504849j
     -2.59079252 -3.55000541j  2.78417111-13.26650134j
     -6.0297906  -1.17150871j -4.92477154 -4.98328949j]
    [-2.49800181e-16+0.        j  1.49800460e+00-7.53097769j
     -2.88537029e-01+0.69659001j -2.36488255e-01+0.35392969j
     -2.22614343e-01+0.22261434j -2.16932373e-01+0.14494958j
     -2.14217111e-01+0.08873163j -2.12937210e-01+0.04235584j
     -2.12556562e-01+0.        j -2.12937210e-01-0.04235584j
     -2.14217111e-01-0.08873163j -2.16932373e-01-0.14494958j
     -2.22614343e-01-0.22261434j -2.36488255e-01-0.35392969j
     -2.88537029e-01-0.69659001j  1.49800460e+00+7.53097769j]
    
[/code]

このような結果になるのは何故でしょうか  
コードのどの部分に問題があるのでしょうか？

上の括弧がサンプル点が8つと仮定した時の結果が私のFFTの結果で  
下の括弧が本家のnp.fft.fftのアルゴリズムの結果です

ご覧の通り、結果は全然違います

自分では何回も検査しましたが、なかなかみつかりません  
バタフライ演算の方向もあっていますし  
wの値のとり方に問題がありそうですが  
自力では見つかりません

理由を教えて貰えないでしょうか？


CC BY-SA 4.0
2018-08-03T12:33:38.327
47221
2018-08-04T10:05:16.767
2018-08-03T13:59:30.497
3060
24226
question
python アルゴリズム 数学
pythonによるFFTの実装


すみません  
先ほど自力で解決しましつ

バタフライ演算の肝をしらなかったからですね

正解を下にコピーします

（FFT_sort.pyの所は同じなので省略 FFT.pyだけの修正）

[code]

    import FFT_sort as fs
    import numpy as np
    from numpy import pi
    
    import matplotlib.pyplot as plt
    
    
    
    N = input()
    
    
    x = np.linspace(0, 2 * np.pi, N)
    y = np.sin(x)
    
    #y = np.ones(N)
    
    data = fs.sort(N)
    
    if data is None:
        print("error")
        quit()
    else:
        print(data)
    
    
    
    result = np.zeros(N, dtype=np.complex)
    temp = np.zeros(N, dtype=np.complex)
    result[:] = y[data]
    
    
    num = 1
    
    print(y)
    print(result)
    
    while(num != N):
    
        count = 0
    
        T = num
        F = N / (num * 2)
    
        print("********************************************")
    
        print("F:", F)
    
    
        temp[:] = result[:]
        print("result:", result)
        print("temp:", temp)
    
        for i in range(N):
    
            judge = i / num
    
            f = iF = i * F % N
    
            if iF - N / 2 >= 0:
                iF -= N / 2 
                flag = -1.0
            else:
                flag = 1.0
    
            w = flag * np.exp(float(iF) * -1j * 2.0 * pi / float(N))
    
    
            if judge % 2 == 0:
                k = i + T
                #n = data[k]
                result[i] += w * temp[k]
            else:
                k = i - T
                #n = data[k]
                result[i] *= w
                result[i] += temp[k]
    
            print(judge % 2 == 0)       
    
    
            print("i:", i)
            print("f:", f)
            print("iF:", iF)
            print("flag:", flag)
            print("k:", k)
            #print("n:", n)
            print("w:", w)
    
            #result[i] += w * temp[n]
    
            #print("temp[n]:", temp[n])
            print("former result[i]:", temp[i])
            print("result[i]:", result[i])
    
        count += 1
        num *= 2
    
    
    
    
    print(result)
    answer = np.fft.fft(y)
    print(answer)
    
    print(result[:] - answer[:])
    
[/code]


CC BY-SA 4.0
2018-08-04T10:05:16.767
47241
2018-08-04T10:05:16.767
24226
47221
answer