中速フーリエ変換 ~離散フーリエ変換よりは速い~: 長続き目指しブログ

<< Calculix 20本ノック:11本目: exampleの実行 Contact/Hertz2D | TOP | YouTubeのBlender公式チュートリアル（英語）で出てきた機能について、ひたすらメモする >>

20200419

中速フーリエ変換 ~離散フーリエ変換よりは速い~

正直遅い

高速フーリエ変換より遅い。
でも、離散フーリエ変換より早い。
だから、「中速フーリエ変換」と名付けとく。

かかる時間
高速フーリエ変換<<<中速フーリエ変換<離散フーリエ変換
こんな感じ。

原理

離散フーリエ変換では、各周波数の三角関数と、入力信号の積の総和計算している。

三角関数の周期性を考慮して、入力信号f(x)を加法定理でまとめてから計算する。
→かけ算の数が減り早くなる。はず。

最後のほうに、exp(-ikn/N)をNを何個か振った状態で載せておくので、にらめっこして原理を察してほしい。
考え方は、途中まで高速フーリエ変換と一緒なはず。

低速フーリエ変換（普通の離散フーリエ変換）

import time
import numpy as np
arr = np.random.rand(1024)

start = time.time()

import math

N = len(arr)
W = [[math.e**(-1j*2*math.pi*k*i/float(N)) for i in range(N) ] for k in range(N) ]
out = [sum([arr[i]*W[k][i] for i in range(N)]) for k in range(N)]

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

中速フーリエ変換

import time
import numpy as np
arr = np.random.rand(1024)

start = time.time()

import math

def order(Ns):
    if len(Ns) == 1:
        return Ns
    N1=Ns[1::2]
    N2=Ns[0::2]
    N2 = order(N2)
    return N1+N2

N = len(arr)
out = [0]*len(arr)#np.zeros(10)

N0 = range(N)
order1 = order(N0)

N1 = arr[:len(arr)/2]
N2 = arr[len(arr)/2:]
N3 = [N1[i]-N2[i] for i in range(len(N1))]
N4 = [N1[i]+N2[i] for i in range(len(N1))]
N5=[]
i1 = 0
i2 = 1

while i1 < len(arr):
    W = [math.e**(-1j*2*math.pi*order1[i1]*i/float(len(arr))) for i in range(len(arr)/2*i2)]
    N5.append(sum([N3[i]*W[i] for i in range(len(N3))]))
    if i1==len(arr)-2:
        N5.append(sum(N4))
        break
    if order1[i1] > order1[i1+1]:
        N1 = N4[:len(N4)/2]
        N2 = N4[len(N4)/2:]
        N3 = [N1[i]-N2[i] for i in range(len(N1))]
        N4 = [N1[i]+N2[i] for i in range(len(N1))]
        i2+=1
    
    i1+=1

for i in order1:
    out[order1[i]] = N5[i]

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

高速フーリエ変換

import time
import numpy as np
arr = np.random.rand(1024)
start = time.time()

import math
def FFT(f):
    N = len(f)
    if N == 1:
        return f[0]
    f_even = f[0:N:2]
    f_odd = f[1:N:2]
    F_even = FFT(f_even)
    F_odd = FFT(f_odd)
    W_N = [math.e**(-1j * (2 * math.pi * i) / float(N)) for i in range(0, N // 2)]
    F = [0]*N
    try:
        F[0:N//2] = [F_even[j] + [W_N[i] * F_odd[i] for i in range(len(W_N))][j] for j in range(len(W_N))]
    except IndexError:
        F = [F_even + W_N[0] * F_odd]
    try:
        F[N//2:N] = [F_even[j] - [W_N[i] * F_odd[i] for i in range(len(W_N))][j] for j in range(len(W_N))]
    except  IndexError:
        F.append(F_even - W_N[0] * F_odd)
    return F

a1 = FFT(arr)

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")