アカデミック

【超初心者向け】pythonで音声認識④「高速フーリエ変換を実装してみよう」

pythonで簡単な音声認識をやってみたいぞ。

そもそも何から始めればいいのかしら。

今回は,基本的な音声認識をpythonで行う方法をお伝えしていこうと思います。本記事はpython実践講座シリーズの内容になります。その他の記事は,こちらの「Python入門講座/実践講座まとめ」をご覧ください。

【超初心者向け】python入門講座/実践講座まとめ目次入門講座 1.実行環境 2.文字の出力 3.データ型 4.変数 5.更新と変換 6.比較演算子 7.論理演算子 8.条件分岐 9.リスト...

コーディングに関して未熟な部分がたくさんあると思いますので,もし何かお気づきの方は教えていただけると幸いです。また,誤りについてもご指摘していただけると非常に助かります。

お題

FFTをpythonで実装してみよう。

流れ

今回対象とするデータは「あいうえお」ではなく,データ長が1024である合成サイン波とします。なぜなら,DFTのデータ長が2のべき乗のときに,効率的な再帰計算により高速フーリエ変換が可能になるからです。

実際に使われているFFTには様々なアルゴリズムが存在し,データ長が2のべき乗でない場合にも適用可能になっています。

実際に合成サイン波を作ってみましょう。

import numpy as np
import matplotlib.pyplot as plt

# データ長
N = 128
x = np.arange(N)

# 周期
t1 = 10
t2 = 20
t3 = 30

# 3つの周期の重ね合わせ
y = np.sin(t1 * 2 * np.pi * (x/N)) + np.sin(t2 * 2 * np.pi * (x/N)) + np.sin(t3 * 2 * np.pi * (x/N))

# 可視化
plt.plot(x, y)
合成サイン波

FFT(高速フーリエ変換)

詳しい導出は省略させていただきますが,簡単に導入しておきます。まず,離散フーリエ変換において$e^{-i\frac{2\pi \omega n}{N}}=W_N^{n}$とおきます。すると,離散フーリエ変換は以下のように表されます。

\begin{eqnarray}
X[\omega] &=& \sum_{n=0}^{N-1} x[n] W_N^{n\omega}
\end{eqnarray}

そして,ここから少し天下り的なのですが,離散フーリエ変換の式を偶数番目のデータ点で行った場合と,奇数番目のデータ点で行った場合に分けます。データ長が2のべき乗であるから,半分ずつに分けられますね。

\begin{eqnarray}
X_{even}[\omega] &=& \sum_{n=0}^{\frac{N}{2}-1} x[2n] W_{\frac{N}{2}}^{n\omega}\\
X_{odd}[\omega] &=& \sum_{n=0}^{\frac{N}{2}-1} x[2n+1] W_{\frac{N}{2}}^{n\omega}
\end{eqnarray}

このとき,$X_{even}[\omega]$,$X_{odd}[\omega]$,$W_N^{n}$に規則性が現れます。この規則性こそが高速フーリエ変換を「高速」たらしめる所以です。

\begin{eqnarray}
X_{even}[\omega] &=& X_{even}[\omega + \frac{N}{2}]\\
X_{odd}[\omega] &=& X_{odd}[\omega + \frac{N}{2}]\\
W_N^{n} &=& -W_N^{n+\frac{N}{2}}
\end{eqnarray}

以上を踏まえると,離散フーリエ変換を再帰的に計算することができます。まず,離散フーリエ変換を奇数部分と偶数部分で表してみましょう。偶数部分はインデックス0から始まりますので,偶数部分に奇数部分の変形版を足し合わせることで元の離散フーリエ変換を実現するという方針です。

\begin{eqnarray}
X[\omega] &=& X_{even}[\omega] + X_{odd}[\omega]W_N^{n}
\end{eqnarray}

この式に先ほどの周期性を適用すると,高速フーリエ変換が完成します。前半部分と後半部分が同時に求まっていくイメージです。

\begin{eqnarray}
X[\omega] &=& X_{even}[\omega] + X_{odd}[\omega]W_N^{n}\;(\omega=0,\cdots,\frac{N}{2}-1)\\
X[\omega+\frac{N}{2}] &=& X_{even}[\omega] – X_{odd}[\omega]W_N^{n}\;(\omega=0,\cdots,\frac{N}{2}-1)
\end{eqnarray}

実装

import numpy as np

def FFT(x):
    N = x.shape[0]
    
    # 再起動作の一番最後
    if N==1:
        return x[0]
    
    x_even = x[0:N:2]
    x_odd = x[1:N:2]
    
    # ここで再帰動作をする
    X_even = FFT(x_even)
    X_odd = FFT(x_odd)
    
    # DFTと同じようにWを求める
    W = []
    for t in range(N//2):
        W.append(np.exp(-1j * ((2*np.pi*t) / N)))
    W = np.array(W)
    
    # ここで型をcomplexに指定しないとエラーを吐くので注意
    X = np.zeros(N, dtype="complex")
    X[0:N//2] = X_even + W*X_odd
    X[N//2:N] = X_even - W*X_odd
    
    return X

結果

# 実行
Y = FFT(y)

# 結果は複素数なので実部と虚部に分けて可視化
Y_real = np.abs(Y)
Y_imag = np.imag(Y)
plt.plot(np.arange(0, Y_real.shape[0]), Y_real)
実部の可視化

DFTと同じ結果になりましたね。

plt.plot(np.arange(0, Y_real.shape[0]), Y_imag)
虚部の可視化

DFTと同じ結果になりましたね。ちなみに,「np.fft.fft」で試しても全く同じ結果が得られます。

高速逆フーリエ変換

せっかくなので逆変換をしてみましょう。高速フーリエ変換の逆変換は非常に単純です。逆変換は,以下のように定義されます。

\begin{eqnarray}
x[n] &=& \frac{1}{N}\sum_{n=0}^{N-1} X[\omega] W_N^{n\omega}
\end{eqnarray}

この式を変形していきます。

\begin{eqnarray}
\frac{1}{N}\sum_{n=0}^{N-1} X[\omega] W_N^{n\omega} &=& \frac{1}{N}\overline{\sum_{n=0}^{N-1} \overline{X[\omega]} W_N^{-n\omega}}
\end{eqnarray}

つまり,「FFTの結果の複素共役を取り」「普通のFFTを行い」「その結果を再び複素共役にして」「Nで割る」ことによって逆変換が可能になります。実際に行ってみましょう。

実装

def IFFT(X):
    X = X.conjugate()
    x = FFT(X)
    return (1/N) * x.conjugate()

今回は見やすいようにデータを作り直してみます。

N = 1024
x = np.arange(N)
t = 10

y_origin = np.sin(t * 2 * np.pi * (x/N))
plt.plot(x, y_origin)
変換前の信号
Y = FFT(y_origin)
y_recover = IFFT(Y)
y_recover_real = np.real(y_recover)
plt.plot(np.arange(y_recover_real.shape[0]), y_recover_real)
復元後の波

時間比較

せっかくなので,DFT・FFT・np.fft.fftの速度を比較してみましょう。

t1 = time.time()
Y = DFT(y)
t2 = time.time()
elapsed_time = t2-t1
print("経過時間:{elapsed_time}")

t3 = time.time()
Y = FFT(y)
t4 = time.time()
elapsed_time = t4-t3
print("経過時間:{elapsed_time}")

t5 = time.time()
Y = np.fft.fft(y)
t6 = time.time()
elapsed_time = t6-t5
print("経過時間:{elapsed_time}")
Output:
経過時間:3.5613720417022705
経過時間:0.04334902763366699
経過時間:0.00020003318786621094

圧倒的にDFTの計算が遅いですね。今回のFFTの実装ではfor文を使ってしまいましたので,nmpyのfftよりも時間がかかってしまっています。

COMMENT

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です