Pythonでお手軽音声処理(5) フーリエ変換・スペクトログラム表示

シェアする

  • このエントリーをはてなブックマークに追加

お久しぶりです。ふるやん(@furuya1223)です。

今回もPython音声処理シリーズです。

前回scipy を使って色々な波を生成して聞いてみました。

今回はついに離散フーリエ変換を扱います。また、スペクトログラムとよばれる画像を表示します。

今回の記事では離散フーリエ変換の理論について詳しくは解説しませんが、以下の書籍に理論と実装について詳しく書かれています。おすすめです。

詳しい理論の解説記事も後日投稿するかもしれません。私にやる気があれば……

できるようになること

  • 離散フーリエ変換をふんわり理解する
  • 正弦波の重ね合わせでノコギリ波を作成する
  • 対数パワースペクトログラムを描画する

離散フーリエ変換(DFT)とは

離散フーリエ変換は「離散的な信号を正弦波の足し算に分解すること」を意味します。

ざっくり解説していきましょう。

正弦波について

突然ですが、正弦波と呼ばれる種類の波があります。

正弦波は、その名の通り sin(x) の関数を用いて表現される波です。

正弦波は、「振幅」「周波数」「初期位相」の 3 つの特徴を決めると 1 つに決まります。

振幅は、振れ幅の大きさです。

周波数は、波の細かさ(1秒間に何往復するか)です。

初期位相は、時刻 0 のときにどの位置から始まるかです。上の図では「0から上がる方」で始まっています。初期位相は「横方向(時間方向)のずれ」とも解釈できます。

初期位相は単に「位相」と呼ばれることもあります。厳密には別物ですが。

位相はとりあえずあまり深く考えなくて大丈夫です。

フーリエ変換の主張

フーリエ(Fourier)さんは「どんな波でも色々な周波数・振幅・位相の正弦波の足し算だけで表現できる」ということを発見しました。

さらに、「どんな周波数・振幅・位相の正弦波を足せば元の波を再現できるか」も計算することができます。この計算方法をフーリエ変換と呼びます。

しかも、元の波が周期的な波(一定間隔で同じ形を繰り返す波)であれば、「繰り返しの周波数(周期の逆数)の定数倍の周波数の正弦波のみで構成できる」のです。つまり離散的な(飛び飛びの)周波数のみ考えればよいのです。

(このときはフーリエ級数展開と呼ばれます。まあ名前はあまり気にしなくていいです)

ノコギリ波を例に確認してみましょう。周波数は 440 Hz です。

(フーリエ級数展開しやすいように、前回作ったものを横にずらしています)

上のような波形を持つノコギリ波は、直線と角で構成されているので、正弦波では作れなさそうです。

このノコギリ波に対してフーリエ級数展開を行うと、以下のような結果が得られます。

この結果をの位相と振幅を「フーリエ係数」と呼びます。

周波数 位相 振幅
440 Hz 0 0.5*2/π
880 Hz 0 -0.5*2/2π
1320 Hz 0 0.5*2/3π
... ... ...
440*k Hz 0 (-1)^(n-1)*0.5*2/kπ
... ... ...

ノコギリ波の周波数と同じ周波数(今回は 440 Hz)を「基本周波数」と呼びます。今回は、基本周波数成分の振幅は、ノコギリ波の振幅に 2/π を掛けた値になります。

基本周波数の整数倍の周波数の成分を「倍音(成分)」と呼びます。2倍なら 2倍音です。

ノコギリ波では、n倍音の成分の振幅は、ノコギリ波の振幅を 1/n 倍して、符号を交互に反転させた値になるそうです。

厳密にノコギリ波を再現するには無限個の正弦波を足さないといけませんが、だんだん振幅が小さくなるので、適当なところで打ち切ってもノコギリ波に近い波が得られそうです。

本当にノコギリ波になるのか、確かめてみましょう。

今回は位相が全て 0 なので無視するとして、基本周波数と、各倍音の振幅の列が与えられたときに、それらの振幅で正弦波を重ね合わせる関数を作成します。

最初 data を空配列とし、前回作成した正弦波作成関数を用いて作成した正弦波を足していきます。

振幅の値は、入力された振幅列 fourier_series[i] とします。周波数は i 倍音なので frequency の i 倍とします。

main関数でノコギリ波のフーリエ係数を作成します。

dim 変数で、正弦波をいくつ足すかを決めます。

まず、1つの正弦波を表示してみます。

普通に正弦波ですね。

次に dim=2

あーなんかノコギリ波っぽくなりましたね。

続いて dim=3

順調にノコギリ波に近づいています。

次は一気に dim=10 にしてみましょう

もうかなりノコギリ波ですね。

最後に dim=50 でやってみましょう。

角の部分がちょっと変ですけづ、まあほとんどノコギリ波です。

音を聞き比べてみましょう。

ノコギリ波

フーリエ係数から作成したノコギリ波もどき(dim=50)

ほとんど同じですね。フーリエ係数から元の波を再現できています。

ちなみに、dim=10 ではこんな感じ。高音が足りないですねぇ。

フーリエ変換について、さらに嬉しい事実があります。

今扱ったノコギリ波は、連続時間信号として扱いましたが、コンピュータで扱う信号は必ずサンプリングされており、飛び飛びの時間にしか値を持ちません。(第2回参照)

このような離散時間信号は、有限の範囲の周波数の正弦波で表せるのです。

第2回で「サンプリング定理」のお話をしました。サンプリング周波数の半分の周波数までしか表現できないという話です。あれのことです。

そのため、コンピュータで扱う音声に関しては、無限の倍音を考える必要が無いのです。

つまり、「離散時間で周期的な信号」を正弦波の重ね合わせで再現するとき、「有限範囲かつ離散的な周波数の正弦波」で完全に再現できるのです。

ということは、データ数が有限であり、フーリエ係数(各周波数成分の振幅と位相のデータ)もコンピュータで扱えるデータになります!

これが離散フーリエ変換(Discrete Fourier Transform: DCT)です。

離散フーリエ変換をやってみる

適当な音声に対して、離散フーリエ変換をやってみましょう。

これをやると、「どの周波数の音がどれくらいの強さ(振幅)で含まれているか」が分かります。

やってみたいのですが、離散フーリエ変換では「周期的な信号」しか扱えないのでした。

実際の音声(人間の話し声など)はほとんど周期的ではありません。さてどうしましょう。

答えは簡単で、「対象の音声をリピート再生させて周期的な信号にする」です。

音声の長さが10秒なら、リピート再生させることで周期10秒(周波数0.1Hz)の周期信号になります!

ただし、実際に知りたいのは「ある瞬間にどの高さの音がどれくらい含まれていて、それがどのように時間変化するか」であることが多いので、音声を短く切り分けて、各区間ごとに離散フーリエ変換を行います。

この短く切った区間を「フレーム」と呼びます。

音声を短いフレームに分割し、各フレームについて「その区間の音声が繰り返されていると考えて離散フーリエ変換を行う」という処理を行います。

短い時間に対してフーリエ変換を行うので、短時間フーリエ変換(Short-Time Fourier Transform: STFT)と呼ばれます。

STFTでは、フレームの長さの他に、フレームをずらす間隔(フレームシフトサイズ)窓関数を指定する必要がありますが、詳しくはググってください。

1フレーム分の信号をリピートすると、つなぎ目で値が不連続に変化してしまいます。それを防ぐため、端で0になるような滑らかな関数を掛け算して、強制的につなぎ目を0にします。このときに掛ける関数を窓関数と呼びます。

シンプルなハニング窓(ハン窓とも呼ばれる)がよく用いられます。librosa.core.stft() のデフォルトもハニング窓です。

ただし、窓関数を掛けるとフレーム区切りの情報が失われてしまうので、それを防ぐために、フレーム分割の仕方を工夫します。1フレーム目と2フレーム目が重なるように、少しずつずらしながらフレームを取るのです。このずらし幅がフレームシフトサイズです。

フレームシフトサイズがフレームサイズの半分のとき、フレームが半分ずつ重なるので「ハーフオーバーラップ」と呼びます。フレームサイズの 1/4 ずつずらすときは「3/4オーバーラップ」と呼ぶこともあります。

フーリエ変換は librosa という音声処理用パッケージを用いるのが便利です。

pip install librosa

もしくは

conda install -c conda-forge librosa

でインストールできるはずです。

フーリエ変換は librosa.core.stft() で実行できます。引数で信号データ (data) とフレームサイズ (fft_size)、フレームシフトサイズ (hop_length) などを指定します。

FFTは高速フーリエ変換の略です。フレームサイズが 2 の累乗(256, 512, 1024 など)の場合、離散フーリエ変換は高速に計算することができ、これを高速フーリエ変換(FFT)と呼びます。普通はFFTを行うので、fft_size は 2 の累乗にしましょう。

stft() の結果は、numpy の ndarray で帰ってきます。shape は (fft_size/2+1, フレーム数) です。

結果の [f, t] の要素には、第 t フレームにおける f 番目の周波数成分のフーリエ係数が記録されています。

周波数番号は 0 から fft_size/2 まであり、0 が 0 Hz、 fft_size/2 が rate/2 Hz(サンプリングレートの半分)を表します。

STFTの結果が 0 Hz から rate/2 Hz の範囲を表すのは常に同じですが、fft_size/2+1 段階で表現されるので、フレームサイズを大きくすると周波数分解能が上がります。

ただし、フレームサイズを大きくすると、「短い瞬間の様子」が得られなくなります。これってなんだか不確定性原理みたいですね。

STFTの結果は ndarray ですが、その要素は複素数です。1つの値で振幅と位相を表現しないといけないため、振幅は絶対値、位相は偏角で表現しています。

そのため、振幅が欲しい場合は numpy の abs() 関数で絶対値に変換します。

np.abs(librosa.core.stft(data, fft_size, hop_length)) で得られた振幅データ(2 次元の ndarray )について、2乗して対数を取ります。(2乗は無くてもいいですが)

人間の聴覚特性に基づいてのことですが、まあ難しいことは気にしなくていいです。気になる人は「ウェーバー・フェヒナーの法則」とかでググってください。

以上の処理で「対数パワースペクトログラム」が得られます。

あるフレームに対する周波数ごとの情報を「スペクトル」と呼び、それを複数のフレーム(時間)で並べたものを「スペクトログラム」と呼びます。

関数にまとめたのがこちら。

この結果を図示するために、librosa.display.specshow() を利用します。

librosa.display.specshow() はスペクトログラムのデータ、サンプリングレート、フレームシフトサイズなどを入力すると、図を描画してくれます。

色の付け方は cmap で変えます。色々変えたい場合は "matplotlib cmap" などでググってください。

main() 関数で音声を読み込んで以上の関数を実行します。

とりあえず、フレームサイズは 1024、シフトサイズはその 1/4 とします。

JSUT コーパスの音声 sample.wav に対して図を書くと、以下のようになります。

横軸が時間、縦軸が周波数です。

色が明るい部分が、その成分を強く含むことを表します。

人間の声には 5000 Hz を超える成分が少ないこと、低周波領域で縞模様を描くことが分かります。

縞模様は、人間の声が、声帯振動の基本周波数とその倍音構造で構成されていることを表します。

縞模様が斜め上向きのところでは、声が高くなっていっています。

これだけでは少し分かりづらいので、前回作成した正弦波とノコギリ波の対数パワースペクトログラムも見てみましょう。

正弦波

一番下の得に明るい線が 440 Hz を表します。

正弦波なので 1 本の線だけになってほしいところですが、フレームサイズで切ってリピートしたり、窓関数を掛けたりしているので、厳密に 1 本になってくれません。まあこういうものです。

ノコギリ波

なんかすごい。

一番下の 440 Hz 以外の細かい横線は、倍音(880 Hz, 1320 Hz, ...)を表します。多いですね。

倍音以外に明るくなっているのは、先ほどと同じくフレームと窓関数の影響でしょう。

ついでに、今回作成した 50 個の正弦波の重ね合わせによるノコギリ波でも表示してみましょう。

440*50=22000 なので、44.1kHz のサンプリングレートで表現できる倍音をすべて含んでいます。

むしろこちらの方がいらない部分が無くてきれいですね。

10個の正弦波の重ね合わせで作成した、物足りないノコギリ波でも見てみましょう。

綺麗に 10 本しか無いですね。正弦波よりきれいなのは何故。

ということで、かなり長い記事になりましたが、離散フーリエ変換をしてスペクトログラムを表示することができました。やったー。

今回使用したコードは GitHub リポジトリで公開しています。

次回は簡単なノイズ除去ができれば良いなと思っています。

では。

スポンサーリンク
レクタングル(大)
レクタングル(大)

シェアする

  • このエントリーをはてなブックマークに追加

フォローする