アカデミック

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

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

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

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

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

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

お題

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

流れ

今回対象とするデータは「あいうえお」ではなく,データ長が2をべき乗に制限した「ドレミファソラシド」にしようと思います。以下のファイルを対象とします。

データ長を調整したドレミファソラシド

短時間フーリエ変換の概要

STFTの簡単な概要は以下のページで説明しています。ここでは,実装方法をメインにお伝えしていこうと思います。

【超初心者向け】wavファイルからpythonでMCMS(時間/周波数分解能の異なるメル周波数スペクトログラム)を抽出する 時間周波数分解能の異なるメル」周波数スペクトログラムを抽出したいゾ。 時間周波数分解能ってどうやて変えるんだろう? ...

実装

前提として,窓関数は既成のライブラリを使用します。また,FFTの実装はこちらのページで説明しているものを利用します。また,wavファイルの読み込みと書き出しにはlibrosaを利用します。

必要なライブラリのインポート

import numpy as np
import matplotlib.pyplot as plt
import librosa
from scipy import hamming

FFTの定義

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

STFTの定義

librosaの公式ドキュメントにもある通り,librosaのSTFTでは窓関数を端っこに適用するときは反転パディングを行っています。今回もそれに習いましょう。また,実信号に対するFFTの結果が対象であることに注意して,後ろ半分の情報を捨ててしまうように実装しています。

def STFT(x, win_length, hop=0.5):
    '''
    今回はwin=Mとして実装。
    データサイズも簡単のため2のべき乗に制限。
    hopはデフォルトで窓幅の半分とする。
    '''
    hop_length = int(win_length * hop)
    
    # 窓関数をかける時に端点が問題になります。
    # 今回はlibrosaのデフォルトに習って反転パディングをしてみます。
    pad_first = x[:hop_length]
    pad_last = x[-hop_length:][::-1]
    x_pad = np.concatenate([pad_first, x, pad_last])

    # データのサンプル数
    N = x_pad.shape[0]

    # FFTの結果を半分にした長さ
    M = int(win_length//2) + 1

    # 窓関数を適用する回数
    T = int((N - hop_length)/hop_length)

    # 今回はハミング窓(ハン窓)を利用します
    han = hamming(win_length)
    
    # 結果を格納する箱です
    spec = np.zeros((M,T), dtype="complex")
    
    for t in range(T):
        # まず窓関数を適用します
        windowed_x = x_pad[t*hop_length:t*hop_length+win_length] * han
        
        # 次にFFTを実行します。
        spec[:,t] = FFT(windowed_x)[:int(win_length//2)+1]
    return spec

データの読み込み

# [path_to_wavfile]は皆さんの環境に合わせて変更してください。
y, sr = librosa.load("[path_to_wavfile]")

STFTの実行

スペクトログラムの描画にはlibrosaを使います。また,こちらのページでdb(デシベル)への変換は説明していますので,今回は振幅からdbへの変換はlibrosaに頼ってしまいます。

spec = STFT(y, 1024, 0.5)
spec_db = librosa.amplitude_to_db(np.abs(spec))
librosa.display.specshow(spec_db, y_axis="log")
「ドレミファソラシド」のスペクトログラム

逆変換

短時間フーリエ変換の逆変換を実装するのも,よい練習になります。今回私も実装してみたのですが,復元された「ドレミファソラシド」の音質が非常に悪くなってしまいました。あとでlibrosaの逆変換と比べてみます。まずは,実装からみてみます。

実装

高速フーリエ変換の逆変換はこちらのページで説明しています。そのままのコードを記しておきます。

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

こちらのIFFTを利用して,本当に文字通りSTFTとは反対の操作を施していきます。

def ISTFT(X, win_length, hop=0.5):
    
    # スペクトログラムの横(時間方向)の長さ
    T = X.shape[1]

    # 窓関数のホップ長
    hop_length = int(win_length * hop)

    # 今回もハン窓を利用します
    han = hamming(win_length)

    # STFTでFFTの結果の半分を捨てていましたので対称性に気をつけながら復元します。
    X_recover = np.concatenate([X[:-1], X[::-1, :][1:]])

    # 改めて復元するデータの長さを計算します
    # STFTでTを計算した式の逆を求めているだけです
    N = int((T * hop_length) + hop_length) - win_length
    
    # 改めて窓関数を適用する回数を計算します
    nm_hop = int((N - hop_length)/hop_length)

    # データを格納する箱です
    x = np.zeros(N, dtype="complex")

    # 端点だけ場合分けしてIFFTを施していきます
    for i in range(nm_hop):
        if i==0:
            x[i*hop_length:i*hop_length + int(win_length//2)] += IFFT(X_recover[:,i])*han[int(win_length//2):]
        if i==X_recover.shape[0]-1:
            x[i*hop_length:i*hop_length + int(win_length//2)] += IFFT(X_recover[:,i])*han[:int(win_length//2)+1]
        else:
            x[i*hop_length:i*hop_length + win_length] += IFFT(X_recover[:,i])*han
    return x

結果

まず,librosaで復元した音源を聞いてみましょう。

spec = librosa.stft(y)
y_recover_librosa = librosa.istft(spec)
librosa.output.write_wav("test_librosa.wav", y_recover_librosa, sr)
librosaによる復元結果

次に,自作関数による復元結果を聞いてみましょう。最後に「np.asfortranarray」でfortranarrayに変換しておかないとエラーを吐くので注意です。

x_recover = np.real(ISTFT(spec, 1024, 0.5))
librosa.output.write_wav('test.wav', np.asfortranarray(x_recover), sr)
自作関数による復元結果。なぜか1MBにもなってしまったので圧縮をかけました。圧縮をかける前後は音質にほとんど変わりはありませんでしたので,こちらが復元結果のほぼオリジナルの音質です。

同じディジタル信号処理のはずなのに,なぜここまで差が出るのかについては要検討です。

COMMENT

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