高速フーリエ変換

高速フーリエ変換という単語に遭遇したので、軽くまとめて、その上でお絵描きする。

はじめましょうか。

  • f:id:medical-science:20220108235851p:plain

    こちらのサイトをまとめる。
  • 多項式乗算問題
    • 入力が多項式2つ。（係数の配列2つ、長さm,n）
    • 出力が、f(x)g(x)。長さm+n-1
      • 1次式（係数は2つ）を2つかけても、2次式（係数は3つ）
    • for loop で解くと、O(nm)
    • 高速フーリエ変換で、O(n logn)
  • 複素数
  • 1の原始N乗根
    • N乗して初めて1になる数
    • \zeta_N = exp^{2 \pi i/N}
    • \Sigma_{i=0}^{N-1}\zeta_N^{i(j-k)}
      • j=kなら、N
      • それ以外なら、0
  • 離散フーリエ変換
    • \hat{f}(t) = \Sigma_{i=0}^{N-1} f(\zeta_N^i) t^i
  • 離散フーリエ逆変換
    • f(t) =\frac{1}{N} \Sigma_{i=0}^{N-1} \hat{f}(\zeta_N^{-i}) t^i
  • これを踏まえて、全体の流れを抑える
    • まず、m+nより、大きな2のべき乗の数を用意
    • f, gに、ζに関して離散フーリエ変換する。
    • そのDFT後の、係数の積を計算する。
      • これは、 \hat{f \cdot g}を求めている。
    • それを離散フーリエ逆変換すれば、求めたい配列が得られる。
  • ここから、高速に離散フーリエ変換を行う方法が出てくる。
    • 2のべき乗を探した理由と繋がる。
    • 基本的に、問題のサイズを半分にして、再帰的に解きたい。
      • この考え方は、DPとも繋がるし、二分探索とも、ダブリングとも繋がる。
      • 分割統治法
        • 分けて分けて、分けては、解く。すると、全体を解いたことになる。
    • DFTを、するときに、半分ずつに分ける方法がある。
    • 詳しい説明はこちら。とてもよかった。
  • ここからは、お絵描きでもする。
 
import scipy.fft
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import *
seed(seed = 21)

alpha = 0.1
N = 5000
T = 1.0/600.0
w = rand(N) * alpha
x = np.linspace(0.0, N*T, N)
y = 0.1*np.sin(60.0 * 2.0*np.pi*x)+0.4*np.sin(90.0*2.0*np.pi*x)
z = 0.2*np.sin(65.0 * 2.0*np.pi*x)+0.3*np.sin(30*2.0*np.pi*x)
y_f = scipy.fft.fft(w+ (1-alpha)*(y+z))
x_f = np.linspace(0.0, 1/(2*T), N//2)

Amp = 2.0/N * np.abs(y_f)
plt.plot(x_f, Amp[:N//2])
plt.show()

  • f:id:medical-science:20220108234819p:plain

    乱数を足してみると、超低周波数の振幅が大きくなった。
  • この性質を利用すれば、低周波数をカットすれば、高周波数の部分はノイズとは考えにくいので、意味ありげな情報を抽出できる（かも）。
  • 最後に、ちょっと乱数の割合をいじってみる。
  • いじった結果が上の画像。
  • 完全に余談だが、FFTした後、半分の長さだけ表示している。
    • これは、もう半分が対称的になっている。
    • 全体の長さに対して、半分。なのは、これは、データが有限だから仕方ない。
    • 10個のデータ、2個ごとの周期と見えるか、8個ごとの周期と見えるか、という問題。

バイバイ！