はじめに
離散フーリエ変換をpythonで0から実装する備忘録です。自分のメモ用なので中身の丁寧さは補償しません。
1つ目参考資料の内容を自分なりに噛み砕いただけなので、もっと簡潔に実装していきたい方はそちらをご覧ください。
参考資料
- https://negligible.hatenablog.com/entry/2021/10/19/182702#%E5%8F%82%E8%80%83%E8%B3%87%E6%96%99
- https://cognicull.com/ja/woy58c08
- https://zenn.dev/yuto_mo/articles/1b2fdde0bf9d79
- https://youtu.be/fGos3wrKeHY?si=WByQu4bTTUciuGzp
- https://qiita.com/white_tiger/items/8cc1031dfe1860795073
DFTへの理解
この章ではDFTについて噛み砕いて理解していきます。
用語
| 用語 | 説明 |
|---|---|
| フーリエ級数展開（フーリエ変換の逆） | 周期的な関数 f(x) を、より簡単な三角関数（サイン波やコサイン波）の和として表現すること |
| フーリエ変換（フーリエ級数展開の逆） | 非周期的な関数 f(x) を、より細かな周波数成分に分解して表現すること |
| フーリエ係数 | フーリエ級数展開やフーリエ変換で使われる、関数を周波数成分に分解する際の重みや係数 |
| 離散 | 連続的な物理量や信号を、一定の間隔でサンプルし、有限個のデータ点で表現すること |
| DFT（離散フーリエ変換） | 離散的な関数（サンプルデータ）を、異なる周波数成分の和として表現すること |
やりたいこと：DFTを使用して時間領域の離散的な関数（サンプルデータ）を、異なる周波数成分の和として表現したい！&&どの時間領域の値がどの周波数成分に寄与しているのかを調査したい！
DFT（離散フーリエ変換）の導出
※この部分は1つ目参考資料を噛み砕いただけです
連続時間での信号をf(t)とします。この時f(t)は具体的に複素数をとります。
f(t)が1周期をTとする周期性を持つとすれば，f(t)をフーリエ級数展開した場合のフーリエ係数F[k]は次式で書くことができます。
F[k] = \frac{1}{T}\int_{0}^{T}f(t)e^{−jkωt}dt
これに対して周期Tを1とみなします。
T=1
またωは1周期がTとなる角周波数、すなわち基本波角周波数であり下記の式が成り立ちます。
ω=\frac{2π}{T} → T=\frac{2π}{ω}=1 また 2π=ω
下記のTを代入して、時間軸を規格化すると
F[k] = \frac{1}{T}\int_{0}^{T}f'(t')e^{−jkωt}dt
= \int_{0}^{1}f'(t')e^{−j2πkt}dt
関数が変わっているためf'(t')と表現しているのですが、もともと使用していたf(t)は今後現れないので、下記のようにf(t)で表現してしまいます。
F[k] = \int_{0}^{1}f(t)e^{−j2πkt}dt
F[k]の0→1の範囲をN分割して、f'(t')をiの関数として時間軸上に離散化します。つまりここでは、現在は積分区間で連続的な値になっているのを離散化して、有限個のデータ点で表現する作業をします。
F[k]の0→1の範囲をN分割すると、tは下記のように表せます。
t[i]=\frac{i}{N}(i=0,1,...,N-1) ※Nはデータの長さ
すると、F[k]も下記のように書き直せます。
t=\frac{i}{N} より
F[k] = \int_{0}^{1}f(t)e^{−j2πkt}dt
= \sum_{i=0}^{N-1}f[i]e^{−j\frac{2πki}{N}}dt
この時に下記のように定義でき、これを回転因子と言います。
e^{−j\frac{2πki}{N}} = W_{N}^{k}[i]
これを用いて表した式が 離散フーリエ変換(Discrete Fourier Transform, DFT) です。
F[k]= \sum_{i=0}^{N-1}f[i]W_{N}^{k}[i]
DFT（離散フーリエ変換）の解釈
離散フーリエ変換(Discrete Fourier Transform, DFT) の式は下記のように示せます。その解釈の仕方を説明します。
F[k]= \sum_{i=0}^{N-1}f[i]W_{N}^{k}[i]
問題となるのは回転因子でしょう。回転因子はこのサイトがとてもわかりやすいので、そちらを見ることをお勧めします。
W_{N}^{k}[i]
回転因子は簡単にいうと、複素平面上の単位円を角度でN分割したとき、2π/Nを1ステップとして、k*iステップ回転させた複素数のことを指します。
では長さN=8の時系列データに対して、k=0,1,...,N-1の範囲で離散フーリエ変換を行うときに具体的にどのような操作になるのかを考えます。まず基本的に複素平面は横軸を実数軸、縦軸を虚数軸とした平面で、そこに半径1となる円を描くことをイメージしてください。N=8なのでその円を8等分します。
k=0のとき、F[k]はiが0からN-1まで変動し、それを足し合わせていくので、
N=8 k＝0 より
F[0]=(f[0]W_{8}^{0}[0])+(f[1]W_{8}^{0}[1])+...+(f[N-1]W_{8}^{0}[N-1])
となります。k=0の時は、複素数平面上の円を2π/N=2π/8より8等分し、その上を
ki=0*0=0ステップ
ki=0*1=0ステップ
ki=0*2=0ステップ...
していったときの複素数の値を全て足し合わせます。k=0の時は全て0ステップになるので、すべての値を足し合わせるだけです。
※このステップの話を聞くとこの動画の解釈がわかってきます。
では、kが0以外の時を考えます。
N=8 k＝1 より
F[1]=(f[0]W_{8}^{1}[0])+(f[1]W_{8}^{1}[1])+...+(f[N-1]W_{8}^{1}[N-1])
となり、k=1の時は、複素数平面上の円を8等分し、その上を
ki=1*0=0ステップ
ki=1*1=1ステップ
ki=1*2=2ステップ...
となり、k=N-1になるまで、少しずつ円の上を進みながら値を足し合わせたものが最終的にフーリエ変換で得られる値になります。
このようにして得られたフーリエ変換の返り値（複素数）は、k番目の周波数として扱われ、下記のような値が得られます。
振幅A = \sqrt{(実部)^2+(虚部)^2}
位相ϕ = \tan^{-1}(\frac{虚部}{実部})
フーリエ変換でよくみられる返り値を使用した横軸が周波数、縦軸が振幅となっているグラフは、通常、フーリエ変換のパワースペクトルや振幅スペクトルを示します。
- 横軸（周波数）：各周波数成分を示し、若い値（左側）の方が低周波を表します
- 縦軸（振幅）：各周波数成分の振幅を示し、その周波数成分の強さを表します
低周波成分は通常、信号の全体的なゆっくりした変動を表し、信号全体に共通するような大きな傾向を捉えます。一方、高周波成分は信号の急速な変動を表します。
実装
いよいよ実装していきます。
F[k]= \sum_{i=0}^{N-1}f[i]e^{−j\frac{2πki}{N}}dt
import numpy as np
from scipy.fft import fft
def original_fft(data):
N = len(data) # データの個数
F = [] # DFTの結果を保存するリスト
# DFTの計算
for k in range(N):
F_k = 0
for i in range(N):
F_k += data[i] * np.exp(-2j * np.pi * k * i / N)
F.append(F_k)
return F
# 入力データ（例として4つのデータ）
f = [1, 2, 3, 4]
fft_res = original_fft(f)
fft_pac = fft(f)
# 結果の表示
print("original fft")
for k in range(len(fft_res)):
print(f"F[{k}] = {fft_res[k]}")
print("package fft")
for k in range(len(fft_pac)):
print(f"F[{k}] = {fft_pac[k]}")
出力は
original fft
F[0] = (10+0j)
F[1] = (-2.0000000000000004+1.9999999999999996j)
F[2] = (-2-9.797174393178826e-16j)
F[3] = (-1.9999999999999982-2.000000000000001j)
package fft
F[0] = (10-0j)
F[1] = (-2+2j)
F[2] = (-2-0j)
F[3] = (-2-2j)
高速フーリエ変換（fft）と比較しており、誤差のせいで値に差が出ています。
original_fft 関数の結果は理論的に正しいですが、手動計算のために浮動小数点の誤差がわずかに現れます
SciPy の fft 関数は、より最適化されたアルゴリズムを使用しており、誤差が異なる形で現れるかもしれませんが、通常は高精度です
（gptより）
大きい誤差ではないので、今回は四捨五入をして誤差の部分をカット。
rounded_fft_res = [np.round(x, 10) for x in fft_res]
rounded_fft_pac = [np.round(x, 10) for x in fft_pac]
print("rounded original fft")
for k in range(len(rounded_fft_res)):
print(f"F[{k}] = {rounded_fft_res[k]}")
print("rounded package fft")
for k in range(len(rounded_fft_pac)):
print(f"F[{k}] = {rounded_fft_pac[k]}")
結果はいい感じに整えられました。
rounded original fft
F[0] = (10+0j)
F[1] = (-2+2j)
F[2] = (-2-0j)
F[3] = (-2-2j)
rounded package fft
F[0] = (10-0j)
F[1] = (-2+2j)
F[2] = (-2-0j)
F[3] = (-2-2j)
他の誤差の対処法も後々実施してみようと思います。