正直遅い
高速フーリエ変換より遅い。
でも、離散フーリエ変換より早い。
だから、「中速フーリエ変換」と名付けとく。
かかる時間
高速フーリエ変換<<<中速フーリエ変換<離散フーリエ変換
こんな感じ。
原理
離散フーリエ変換では、各周波数の三角関数と、入力信号の積の総和計算している。
三角関数の周期性を考慮して、入力信号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]")
N=2の場合のexp(-ikn/N)
左が実数で右が虚数。
青丸がサンプル点。
N=4の場合のexp(-ikn/N)
N=8の場合のexp(-ikn/N)
N=16の場合のexp(-ikn/N)
おれはいったい、何をやっているのだ。
【pythonの最新記事】
- 断面二次モーメントを、座標点の配列から計..
- 断面二次モーメントを、座標点の配列から計..
- fontファイルの文字データ(グリフ)を..
- matplotlibのpyplot.pl..
- 計算力学技術者試験の問題集 自炊(裁断→..
- pythonで、ホワイトノイズやピンクノ..
- 脳ドッグに行ってきた。→MRIの画像デー..
- matplotlibのimshowで円を..
- matplotlibの、cmapを、徐々..
- matplotlibのmake_axes..
- matplotlib floatinga..
- matplotlib plotの色を、値..
- Pythonで、「二次元フーリエ変換した..
- matplotlibのlinestyle..
- どちらが正しいRGBか。(matplot..
- matplotlibのannotateの..
- matplotlibで、x軸とy軸の数字..
- VBAで、pythonのrangeとか、..
- matplotlibのaxes3Dで、a..
- matplotlibのlatexで、行列..
https://wakelet.com/wake/ZdNBOrf7XFfswNliLP6qi
https://wakelet.com/wake/UBu_63OT6I6R7edPdpbf4
https://wakelet.com/wake/fXWQRKYjvvMHTC0gy3f2d
https://wakelet.com/wake/OXd3Yjky00HltvRoiMr7J
https://wakelet.com/wake/a7Sv1qVlRUsjHL4XAlMDs
https://wakelet.com/wake/DtVVpqPKtTV_opy5Z0jzy
https://wakelet.com/wake/63HqQeSULCUOyM5TvjY72
https://wakelet.com/wake/cCLBl8dRmoBSl_aqtVTcN
https://wakelet.com/wake/bZAKf68SLsmyH8rnFjGyY
https://wakelet.com/wake/OuBFhJLdbIMCs9REV0JSP
https://wakelet.com/wake/JBv8LRjUvzYwM7D-4Bbux
https://wakelet.com/wake/BJed_l7YObA61r2AlbmCa
https://wakelet.com/wake/U1K6_LScPmEyyL5Ynt60l
https://wakelet.com/wake/7G_3eQ19W2sa-qwBiBdbN
https://wakelet.com/wake/g8QgnJkFFYesOV_dSMXo6
https://wakelet.com/wake/sKFo1O0LE7M57iwMjK7m3
https://wakelet.com/wake/qvq8FF4HFF1qCRGCHyeEE
https://wakelet.com/wake/U4ySse10la9T4ZQrE7PXq
https://wakelet.com/wake/PH_9hC92J1nU89YP9PQGi
https://imgur.com/kDAWw4B
https://imgur.com/UDNikao
https://imgur.com/cbxdjpe
https://imgur.com/fQohWYg
https://imgur.com/EX4gGIA
https://imgur.com/zSsOyEJ
https://imgur.com/Jc5sLNB
https://imgur.com/su0hlL0
https://imgur.com/ZR4KI8L
https://imgur.com/OnPRBDh
https://imgur.com/tQvfvJ5
https://imgur.com/D0VSMEq
https://imgur.com/B7vmmyI
https://imgur.com/mF8KNqG
https://imgur.com/L4M9yng
https://imgur.com/MPFsE51
https://imgur.com/WUrD2Ml
https://imgur.com/QcRSOgS
https://imgur.com/fnTMbay
https://imgur.com/qa4NAJ3
https://imgur.com/fC6emza
https://imgur.com/Aj7EydX
https://imgur.com/AKsuabg
https://imgur.com/V7m7pga
https://imgur.com/4UWHGth
https://imgur.com/LAig5w1
https://imgur.com/VJFw2L6
https://imgur.com/t1l7NPl
https://imgur.com/4FVPhsL
https://imgur.com/9qT7A3O
https://imgur.com/x3OOYZ1
https://imgur.com/787pqOx
https://imgur.com/z9EzKUW
https://imgur.com/ThL9zIu
https://imgur.com/SMj4D8b
https://imgur.com/uqOJ6JF
https://imgur.com/X2Dfo63
https://imgur.com/iGPrWN2
https://imgur.com/WZ9JsdV
https://imgur.com/AfncFGX
https://imgur.com/6vHbSeb
https://imgur.com/reOf0fv
https://imgur.com/jxxK5fG
https://imgur.com/LP40ECO
https://imgur.com/7FBEZ48
https://imgur.com/mBJvMLl
https://imgur.com/VxhDjlA
https://imgur.com/yxwOtvn
https://imgur.com/xFWO7lL
https://imgur.com/BvzvbkR
https://imgur.com/2rGMEdH
https://imgur.com/sW0LSIc
https://imgur.com/entFfV5
https://imgur.com/v6r4rXx
https://imgur.com/i0miNG7
https://imgur.com/3SbsSsN
https://imgur.com/Wdq039S
https://imgur.com/xkhyOTC
https://imgur.com/I8Gi9FM
https://imgur.com/ooCHJei
https://imgur.com/JtxPMwE
https://imgur.com/eN3hHZO