ブログ名

スペクトル解析補足記事【付録②】

はじめに

本記事は、「音声認識技術の実践:スペクトル解析」で説明しなかったプログラムに関する詳しい解説をします。

スペクトル解析記事のプログラムはコメントアウトによるメモを多く記載してあります。 本記事は、メモだけではプログラムがよく理解できない!という方に向けたものです。

Pythonによるプログラム構築をしたことがない初心者でも理解できるような解説に努めます。 よろしくおねがいします。

<今回の記事の内容>

  • リファレンス一覧
    実践記事で使用したライブラリのリファレンスを載せておきます
  • プログラムの解説1
    前回の記事で書いたプログラムについて、Pythonの技術に触れながら詳しく解説します。
    • コメントアウトについて
    • [必要なライブラリをインポート]import xx as xx というおまじないについて
    • [超簡単な正弦波を作る]変数の定義について
    • [超簡単な正弦波を作る]ライブラリの関数について
    • [グラフに表示①]matplotlibについて
    • [フーリエ変換]フーリエ変換をする関数
    • [グラフに表示②]茎グラフの表示
  • プログラムの解説2
    • [必要なライブラリをインポート]from xx import xxについて
    • [waveファイルから必要な情報を読み込む]
    • [一旦グラフにプロット]
    • [1024(N)個のサンプル値で構成され、1つ前のスライスとサンプル100(step)個分ずれるようなスライスに分割する]
    • [列ごとに1つのスライスがあるほうが便利なのでデータを転置する]
    • [...グラフにプロットする]
  • プログラムの解説3
    • [waveファイルから必要な情報を読み込む]
    • [グラフに表示する]

リファレンス一覧

英語ですが、関数について解説されているページを参照にすれば問題なくソースコードを理解できるという人は、下のページを参考にしてみてください。

プログラム解説1

コメントアウトについて

まず初めにコメントアウトに関する簡単な説明をします。

コメント文とは、プログラムとして読み込まないようにしてある文章です。 下の表記方法を用います。メモや記録として使われます。 '''これがコメントです'''(「'''」で挟む方法。行数が伸びても大丈夫) #これもコメントです(「#」の後ろに続ける。改行するとコメントアウトは無効になる)

今回はプログラムに全く触れたことがない人も想定して構築しました。 そのため、プログラムの1行目にコメントアウトで全体の処理について説明してあります。そしてソースコード内の処理ごとに段落分けを行い、先頭に処理内容をコメントアウトとして記録までしています。 ソースコード解説も、この段落をひとまとめと考え、説明をしていきたいと思います。

プログラムは自分で利用するだけなら、自分にとって読みやすいものがよいです。 先に説明した通り、コメントアウトはプログラムとして読み込まないものなので、プログラムからすべて消してしまっても問題ありません。

では、1つ目のプログラムを見てみます。スペクトル解析記事では「fft_numpy_sinwave.py」という名称です。

[必要なライブラリをインポート]import xx as xx というおまじないについて

環境構築の時に音声認識に必要なライブラリをインストールしましたね。 Pythonのプログラムにおいてライブラリ内の変数を使用する際は、まず上のおまじないを唱える必要があります。
たとえばimport numpyは 「このソースコードではnumpyというライブラリを利用するよ~」という宣言といったところです。
そして、as xxというのはこのソースコードファイルの中でライブラリをこのようなあだ名で呼びますという宣言です。 [fft_numpy_sinwave.py]において import numpy as np という呪文を唱えた後、np.sin()といった文章が現れますね。
これはNumpyというライブラリを今後はnpという名称で呼ぶと宣言したので、このライブラリ内にあるサイン波を生成する変数を呼び出す際、np.sin()という呼び出し方をしたのです。

[超簡単な正弦波を作る]変数の定義について

ここではプログラムにおける変数の定義について少し説明します。

プログラムにおいて変数は「箱」というイメージで表現されることが多いです。 プログラムを構築する際に自分で設定することができる記憶領域のようなものです。

たとえば、 f = 10 # 周波数f(frequency),単位はHz(サイクル毎秒) このような記述があります。
この時点で「fは周波数のことだよ!」と定義しているわけではありません。
あくまでfという箱を自分で勝手に宣言し、その中に10という値を入れたと考えてください。
変数は英数字で表現します。自分で名称づけるので自分がわかりやすい名称にするとよいでしょう。

t = np.linspace(0, 2, 2 * f_s, endpoint = False)# 時間を0から2秒に指定しサンプリング分割を行う(終点は要素に含まない)
x = np.sin(f * 2 * np.pi * t) # x=sin(2πf)の波形作成

これらの文も同様のことをしています。 「tという箱にnp.linespaceという計算処理を行った結果を入れる」 「xという箱にnp.sin()で計算した結果を入れる」 下にイメージとしての図を載せますが、変数には1つ以上の値を入れることもできます。 f:id:iTD_GRP:20200603234331j:plain

[超簡単な正弦波を作る]ライブラリの関数について

np.linspace()np.sin()はnumpyというライブラリに入っている関数です。 ()の中に計算してほしい数値や関数の設定などの入力ができます。 (中には入力をしなくても動いてくれる関数もあります)

入力や設定・詳しい関数の説明などはリファレンスにて説明されています。

[グラフに表示①]matplotlibについて

ここではmatplotlib(このソースコードではplt)を用いたグラフ表示を行っています。 まず最初の文章についても、matplotlibにおけるお約束の部分です。グラフプロット機能において、plt.subplots()という呪文は、「これからグラフの設定を行います」という宣言のようなものです。

私もこのプログラム(fft_numpy)を作成しているときはよくわからないままにライブラリを使用していました… きちんと調べてみたところ、subplotという関数を使うことで、1つの画面に複数のグラフを表示できるようなんですね!すごいですね。 こうして調べなおすと私が書いたプログラムも実はかなり機能を使いこなせていないかもしれません。

matplotlibに関する機能も掘り下げていくと非常に長く、脱線してしまうので機会があれば別の記事にしたいと思います。

ax.plot(t,x)について解説します これはx軸にtの値、y軸にxの値をplotするという意味です。 (このx, tは[超簡単な正弦波を作る]で定義した変数のことです) set_xlabelはx軸の名称を設定する変数です。 この辺りは図を見たほうがわかりやすいかもしれません

f:id:iTD_GRP:20200603234553j:plain

これらは()内で設定する関数なので、プログラム初心者はこう言った部分を自分で変更してみると動きがわかってよいと思います。

ここまでの内容が理解できましたか? スペクトル解析記事で使用したプログラムは、 * 変数に計算処理した結果を入れる * グラフに表示する

のみを行っているといえます。 従って、上の処理がなんとなくイメージできればプログラム内容がわかると思います。

[フーリエ変換]フーリエ変換をする関数

フーリエ変換を行う関数はnp.fftですが、

などフーリエ変換に関する様々な計算処理があります (詳しくは下のサイトにあります)
https://docs.scipy.org/doc/numpy/reference/routines.fft.html#module-numpy.fft

fftという機能の中から離散フーリエ変換を行うという設定のために、np.fft.fftという書き方をしているということですね。

また、 freqs = np.fft.fftfreq(len(x)) * f_s # サンプリング周波数を返すlen(x)ですが、これは変数のサイズを調べるための関数です。 Pythonにもともとある関数です。

[グラフに表示②]茎グラフの表示

①のグラフ表示で説明した部分と異なる点は下の関数です。 ax2.stem(freqs, np.abs(X), use_line_collection=True) ここでは、茎グラフを作成しているので、matplotlib内のstemという変数を利用しています。 np.abs()は絶対値を返す関数です。 use_line_collectionという設定は、茎部分の高速表示に必要な文章ですが、必須ではありません。

また、

ax2.set_xlim(- f_s / 2, f_s / 2)
ax2.set_ylim(-5, 110)

これらは、x軸、y軸の表示する範囲を指定する関数です。

plt.show()で作成したグラフを表示しています。 念のためグラフ結果②も載せておきますね。

f:id:iTD_GRP:20200603235009p:plain

これで一通りの解説をしました。一度プログラムを載せましたので、全体の流れをもう一度つかんでみてください。

'''簡単な正弦波を作ってnp.fft.fft()を試す'''

'''必要なライブラリをインポート'''
import matplotlib.pyplot as plt
import numpy as np

'''超簡単な正弦波を作る'''
f = 10 # 周波数f(frequency),単位はHz(サイクル毎秒)
f_s = 100 # サンプリングレート。1秒あたりの測定数

t = np.linspace(0, 2, 2 * f_s, endpoint = False)# 時間を0から2秒に指定しサンプリング分割を行う(終点は要素に含まない)
x = np.sin(f * 2 * np.pi * t) # x=sin(2πf)の波形作成

'''グラフに表示①'''
fig, ax= plt.subplots()
ax.plot(t, x) # 横軸にt, 縦軸にxの値をプロットする
ax.set_xlabel('Time[s]')
ax.set_ylabel('Signal amplitude')


'''フーリエ変換'''
X = np.fft.fft(x) # 波形のフーリエ変換
freqs = np.fft.fftfreq(len(x)) * f_s # サンプリング周波数を返す

'''グラフに表示②'''
fig2, ax2= plt.subplots()
ax2.stem(freqs, np.abs(X), use_line_collection=True) # フーリエ変換の結果をステム(茎)プロットする
ax2.set_xlabel('Frequency in Hertz[Hz]')
ax2.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax2.set_xlim(- f_s / 2, f_s / 2)
ax2.set_ylim(-5, 110)

plt.show()

プログラム解説2

おそらくここまでの説明を理解していればこの先はおさらいのように見ることができると思います。

プログラムが理解できたと感じた人は、変更できる値などをいじることで、動き方の違いを確かめてみてください。 次のプログラムを見てみましょう。「fft_voice_spectrogram.py」です。

オーバーラップしてフーリエ変換を行うことでスペクトログラムを作成するプログラムですね。

[必要なライブラリをインポート]from xx import xxについて

ライブラリのインポートで異なる表現がありましたね。 from xx import xxは、ライブラリの中から特定の関数のみ呼び出しています。 例えば from scipy.io import wavfile は、「scipy.ioからwavfileだけを使います」ということになります。

[waveファイルから必要な情報を読み込む]

早速scipyの機能であるwavfile.readを使って音声データを読み込んでいますね。 rate, audio = wavfile.read(r"音声データの置いてあるパス/音声ファイル名.wav", "r") ここでは=の左にある変数が2種類ありますね。
リファレンスも下に載せますが、wavefile.readという関数は「wavファイルからサンプルレートデータを出力する」関数です。そのため、出力された値を格納する変数を2つ用意しています。

また、=の右側、wavefile.readの後ろの()のはじめのrについて少し捕捉します。 執筆者はプログラム構築時はwindows7の環境でした。
wavefile.readの()内は音声データのパスと読み込みを意味する後ろの"r"のみを書いて実行したところ、エラーとなってしまったのです。
どういうことかはわかりませんが、時たま実行環境によって読み込みであったり文章を書いた文字の種類などでエラーが起きてしまうことがあります。
この場合は調べたところ、始めにrを入れると直るとあるサイトで見かけ、実際に入力して実行できました。

原因不明のエラーが起きてしまうこともしばしばあるので、処理ごとに実行してみて動きを確認することも小さなことではありますが、重要だと思います。(こういう作業はデバッグと呼ばれたりします。)

また、np.mean()は平均を出力する関数です。この関数でステレオ音声(左右で異なる音声データ。2値)を1つのデータに変換しています。

[一旦グラフにプロット]

これは上で説明したので、特に必要ないと思います。
実行すると下のようなグラフが表示されます。

f:id:iTD_GRP:20200603235423p:plain

[1024(N)個のサンプル値で構成され、1つ前のスライスとサンプル100(step)個分ずれるようなスライスに分割する]

これは窓関数を作成する部分ですね。
slicesでは、音声データを1024のスライスに分割しています。これは100ステップごとにずらしてスライスしています。
そして、スライスに窓を掛け合わせることでオーバーラップを表現しています。 以前表示したこの図の作業を行っているということですね。

f:id:iTD_GRP:20200603235502p:plain

ここで、slices=slices*windowのように左の変数を右でも表示しています。少し混乱するかもしれません。
slicesという箱の中に入っていた数値を取り出してwindowと掛け合わせた計算結果をslicesの中に入れなおしていると考えてみるとよいかもしれません。 変数は一時的に数値や計算結果を記憶している箱なので、更新すれば新しい値で上書きされ、以前の数値はなくなります。

[列ごとに1つのスライスがあるほうが便利なのでデータを転置する]

これは、numpyの計算処理技術の話になってしまうので、理解できなくても大丈夫です。
numpyは関数を呼び出せば簡単に、かつ高速に計算処理をするという機能を備えたライブラリです。
高速な計算処理を可能にする技術として、計算を行列によって行っています(プログラムではn次元配列と呼ばれたりもします)。

f:id:iTD_GRP:20200603235526p:plain

slices.Tでは行列の行・列を逆転する転置を行っています。
スライスが列にあったほうがイメージしやすいし、スペクトログラムにするときもやりやすいようです。

[...グラフにプロットする]

文章が長かったので省略しました。
ここでは、対数値を取っています。極端に大きい値と小さい値の両方が含まれるときは対数値を取ることで範囲を圧縮することが多いです。多分このくらいの理解でも十分だと思います。

また、ここではグラフのプロット方法が違いますね。
ax.imshow()は、画像または配列を表示する際のプロット方法です。
先ほど少し説明した通り、slicesをフーリエ変換し、それを圧縮した結果であるSは配列値として計算結果が格納されています。
配列を画像としてプロットするために使用した関数でした。
(表示されている結果からもわかる通り、x軸、y軸の他に色としてのz軸と考えられる情報がプロットされています)

ax.axis()はグラフの座標軸の設定です。今回は、最大値と最小値を見えるように表示する設定を行っています。

では、最後に全体のプログラムを載せます。

'''fftを使って短時間フーリエ変換(stft)'''

'''必要なライブラリをインポート'''
from scipy.io import wavfile
import numpy as np
import matplotlib.pyplot as plt
from skimage import util

'''waveファイルから必要な情報を読み込む'''
rate, audio = wavfile.read(r"音声データの置いてあるパス/音声ファイル名.wav", "r") # waveファイルを読み込み用ファイルとして開く

audio = np.mean(audio, axis=1) # ステレオをモノラルに変換
N = audio.shape[0] # 要素数を返す
L = N/rate # 長さの算出

'''一旦グラフにプロット'''
f, ax = plt.subplots()
ax.plot(np.arange(N)/rate, audio)
ax.set_xlabel('Time[s]')
ax.set_ylabel('Amplitude[unknown]')
plt.show()

'''
1024(N)個のサンプル値で構成され、1つ前のスライスとサンプル100(step)個分ずれるようなスライスに分割する
'''
M = 1024 # スライスのサンプル数

slices = util.view_as_windows(audio, window_shape=(M,), step=100) # スライス分割

window = np.hamming(M+1)[:-1] # 窓関数の生成
slices = slices * window # スライスに窓を掛け合わせる

'''列ごとに1つのスライスがあるほうが便利なのでデータを転置する'''
slices = slices.T

spectrum = np.fft.fft(slices, axis=0)[:N] # スライスをフーリエ変換
spectrum = np.abs(spectrum) # スペクトラムの絶対値を取っている

'''
スペクトルには、非常に大きい値と小さい値の両方が含まれるため、対数値を取って範囲を圧縮する
そしてグラフにプロットする
'''
f, ax = plt.subplots()

S = np.abs(spectrum) # スペクトラムの絶対値を取る
S = 20 * np.log10(S / np.max(S)) # スペクトラムの対数値を取る

ax.imshow(S, origin='lower', cmap='viridis', extent=(0, L, 0, rate/2/1000))
ax.axis('tight')
ax.set_ylabel('Frequency[kHz]')
ax.set_xlabel('Time[s]')

plt.show()

プログラム解説3

これで最後です! 「scipy_stft_spectrogram.py

これはscipyのスペクトログラムを作る関数を使ったものでした。

[waveファイルから必要な情報を読み込む]

ここまで理解できた人は、ファイルの読み込みや変数の定義、関数を呼び出した計算処理などがコメントアウトだけでも理解できるようになったのではないかと思います。

少しだけ最後の文について説明しますね。

freqs, times, Sx = signal.spectrogram(y, fs=framerate, window='hanning',
                                      nperseg=1024, noverlap=N-100,
                                      detrend=False, scaling='spectrum') # スペクトログラム変数

実は上のプログラム2つでも少しだけ出ていたのですが、この関数では計算したいデータを入力するだけでなくパラメータの設定も行っています。
例えば、計算するときの窓関数をハニング窓に設定していますね(window='hanning')。 これらの詳しい説明もリファレンスに書かれています。 https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.spectrogram.html 入力しない場合も大体デフォルトの設定で計算してくれるので便利ですね。

[グラフに表示する]

ここで異なるグラフ表示の方法はこれですね。 ax.pcolormesh(times, freqs/1000, 10* np.log10(Sx), cmap='viridis') カラーメッシュという連続変化する色で配列を可視化するためのグラフ表示方法です。 x,y,z軸の情報を軸と色で表現していますね。

imshowとの違いについて調べました。 どうやらpcolormeshは色と座標を結びつけることができる様です。 ただし、imshowのほうが計算処理が速く、データを一定の間隔で配置することができるとのことでした。 pcolormeshのほうが何かと柔軟性があるという記述が多くみられました。

では最後にプログラム全文です。

'''
短時間フーリエ変換(stft)を関数で実行
'''

import wave
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

'''waveファイルから必要な情報を読み込む'''
wavefile = wave.open(r"C:\...\VoiceRecognition\sample\sample_voice01.wav", "r") # waveファイルを読み込み用ファイルとして開く

nframes = wavefile.getnframes() # フレーム総数を調べる
framerate = wavefile.getframerate() # フレームレート(サンプリング周波数)を調べる

y = wavefile.readframes(nframes) # フレームの読み込み
y = np.frombuffer(y, dtype="int16") # 波形を変換出来る様に変形する
t = np.arange(0, len(y))/float(framerate) # 音声データの長さ(x軸)

wavefile.close() # waveファイルを閉じる

N=1024

freqs, times, Sx = signal.spectrogram(y, fs=framerate, window='hanning',
                                      nperseg=1024, noverlap=N-100,
                                      detrend=False, scaling='spectrum') # スペクトログラム変数
```グラフに表示する```
f, ax = plt.subplots()
ax.pcolormesh(times, freqs/1000, 10* np.log10(Sx), cmap='viridis')
ax.set_ylabel('Frequency[kHz]')
ax.set_xlabel('Time[s]')
plt.show()

まとめ

今回は、スペクトル解析の実践で作ったプログラムについて解説しました。 個人的にPythonは慣れてくるとわかりやすい言語だと思います。 ただ、関数内でどのような計算処理が行われてるか見えないので、複雑なモデルを作る際は注意が必要かもしれません。
(それでもリファレンスがあるのでわかりやすいかと思いますが)

プログラムや関数理解の助けになれたら幸いです。

参照したページ

https://numpy.org/doc/stable/ https://docs.scipy.org/doc/scipy/reference/# https://matplotlib.org/#


前の記事へ 目次へ戻る