2019-04-10

教師なしのクラスタリング系アルゴリズム

機械学習クラスタリング

教師なしのクラスタリング系アルゴリズムまとめ。

・したいこと教師なしクラスタリングをしたい。いくつかのクラスタ、どのクラスタにも属さないクラスタというラベルをつけたい。さらにクラスタの数もあらかじめ指定したくない。

・クラスタ数の自動推定

qiita.com

blog.kzfmix.com

・クラスタに属さないデータをノイズとして捉える

qiita.com

2019-02-03

短時間フーリエ変換

x軸が時間、y軸が振幅を表す波形にフーリエ変換(離散フーリエ変換)を施し、その絶対値をプロットするとx軸が周波数、y軸が強度を表す。
波形の全時間に対してどれくらいの周波数がどの程度含まれているのかを知りたい場合はこれで十分。
しかし、どの時間にどれくらいの周波数が含まれているかを知りたい場合には使えない。
それを解決する手法が短時間フーリエ変換(Short -Term Fourier Transformation, STFT)である。
こちらの説明が分かりやすい。
https://www.jstage.jst.go.jp/article/jasj/72/12/72_764/_pdf

通常のフーリエ変換に対して一定の区間だけで大きさを持つ窓関数を元の波形にかけてフーリエ変換を行う。
一度フーリエ変換が終われば窓関数を一定の幅だけスライドさせて元の波形にかけてフーリエ変換を行う。
これを繰り返すことにより時間と周波数、強度の3つの関係を表す図(スペクトログラム)を作成できる。

窓関数を元の波形に対して作用させると

${ \begin{align} x_m(t - mS) = w_a(t-mS)x(t) \end{align} }$

${w_a}$ が窓関数でその長さを ${N}$ とする。つまり、 ${0 \le t\le N-1 }$ の範囲だけ値をとり、その他は0である。
このtは窓関数の時間座標から見たときの尺度なので、元の音源に対しては ${mS \le t \le t mS + N -1 }$ となる。
ここで ${S}$ は窓関数をスライドさせる幅である。一度窓関数をスライドさせるごとに ${m}$ が1ずつ増加していく。
${m}$ は窓関数を当てた後のスペクトログラムのフレーム数のインデックスに対応する。

これに対してフーリエ変換を施すので

${ \begin{align} X(m, k) = \sum^{mS + N-1}_{t=mS} x_m(t - mS) e^{- \frac{2\pi (t-mS)}{N} i } \end{align} }$

となる。
この式から窓関数のサイズ ${N}$ を大きく取るとより細かく周波数成分を調べることができる。
つまり周波数成分の分解能を高くすることができる。周波数と時間の分解能にはトレードオフの関係があることに注意する。 wikiのExampleの比較も分かりやすい。

https://en.wikipedia.org/wiki/Short-time_Fourier_transform

スペクトログラムはx軸に ${m}$ を、y軸に ${k}$ をとってz軸に ${|X(m, k) |}$ をプロットしたものになる。 3次元の図となるが基本的にはz軸を色の濃淡で表した2次元の図で表されることが通例である。

スペクトログラムもlibrosaで簡単に作成できる。

import librosa
import librosa.display
import matplotlib.pyplot as plt

sr = 44100
music, sr_ = librosa.load("fanfare.wav", sr=sr)
D = np.abs(librosa.stft(y))
plt.figure(figsize=(10,4))
librosa.display.specshow(librosa.amplitude_to_db(D, ref=np.max), y_axis='log', x_axis='time', sr=sr)
plt.title('Power spectrogram')
plt.colorbar(format='%+2.0f dB')
plt.tight_layout()

結果

f:id:butch416:20190203001522p:plain

(librosa.display.specshowではsrを指定しないと意図しない値になるので注意)

Dの次元は(周波数、時間、強度)になっている。 n_fftが窓関数のサイズ。上の数式で説明した長さ ${N}$ に対応する。デフォルト値はn_fft=2048。
hop_lengthがスライドさせる幅。上の数式で説明した長さ ${S}$ に対応する。デフォルト値はhop_length=n_fft/4。
win_lengthは窓関数の中でも実際に窓関数を適用させたい長さかな？基本的にはwin_length=n_fftで、win_length<=n_fftである。
大きさをn_fft以下に設定した場合は0でパディングされる。
windowは窓関数の種類。デフォルトではhannになっている。

n_fft( ${N}$ )をデフォルト値よりも大きくしてみる。
これは周波数成分の分解能を高くすることを意味している。
n_fft = 2048 * 2で同様の計算をすると

f:id:butch416:20190203002005p:plain

以前と比較するとより正確に周波数成分が抽出されていることが確認できます。(横縞が見えやすくなっている。)
窓関数のサイズを2倍にしたので時間は半分になっています。

もっと大きくして、n_fft = 2048 * 8でやると

f:id:butch416:20190203003036p:plain

より横縞がくっきり見えるようになりましたが、その分時間軸の境目の描写が曖昧になったような気がします。

この3つのスペクトルの配列を見てみると

n_fft=2048: (1025, 1007)
n_fft=2048 * 2: (2049, 504)
n_fft=2048 * 8: (8193, 126)

となり、時間軸の刻み幅が小さくなっていくことがわかる。時間の分解能が低く、周波数の分解能が高くなっていることを意味している。

2019-02-02

デジタル波形について

音声処理

仕事で音声処理をする機会があったので基本的な知識をまとめておく。

信号の種類

アナログ信号：手を叩いて実際に我々の耳に入ってくる音(連続値)
デジタル信号：手を叩いてマイクで集音しコンピュータ上で処理ができる音(離散値)

アナログ信号は連続値を取りますが、デジタル信号はアナログ信号をデジタル化(0と1で表す)したものであるので離散値を取ります。
波形のグラフはx軸に時間、縦軸にy幅をそれぞれ表します。
よく見るウニョウニョとした波の図です。

デジタル化

音のデジタル化は次の3つのステップを踏むようです。

標本化(サップリング)：時間軸の刻み幅(細かさ)を決める
量子化：振幅の刻み幅を決める
符号化

・参考
https://www.infraexpert.com/study/telephony2.html http://www.soraotona.net/weblog/?p=319

標本化(サンプリング、サンプリングレート)についてはこちらの説明が分かりやすいです。 http://www.aichi-c.ed.jp/contents/joho/contents/jissyuu/049/sound.files/oto.htm

以下は個人的な解釈、理解、要約です。
音の高さ(周波数)は1秒間の振動数で決まります。
振動数とは波形がx軸と交わる回数のことです。
よりx軸と多く交わっている波形の方が高い音になります。
サンプリングレートが低い場合は波形のx軸の刻み幅を大きくなってしまいます。
この場合、x軸と交わった回数を調べる精度が荒くなってしまうため、実際の信号の持つ周波数よりも低い周波数しか捉えることができなくなってしまいます。
ゆえにサンプリングレートを高くした方が高周波数の音を捉えることができ、実際の信号に近い音をデジタル化できるようです。

サンプリングレートを8kHzに設定すると1秒間に1/8000秒ごとの振幅データが得られます。
1秒の音源の場合のデータ数は8000個、2秒の場合は16000個になります。
逆にデータ数をサンプリングレートで割れば音源の再生時間が分かりますね。

サンプリング定理(Shannon-染谷の定理)によりサンプリングした音をアナログ信号に戻す場合の周波数は、サンプリングレートの0.5倍までしか復元できないようです。
8kHzでサンプリングした音をアナログ信号に復元すると4kHzになってしまいます。
なので最高周波数の2倍以上の周波数でサンプリングする必要があるようです。
(2点だけ波形の谷と山を観測できれば少なくとも1点はx軸と交わっているので、2倍以上のサンプリングをすれば良い的なロジックだったと思います。)

python上で音源の加工をする時はLibROSA(https://librosa.github.io/librosa/)というモジュールが便利です。
短時間フーリエ変換を行い簡単にスペクトログラムを得ることもできます。

音源はこちらのfanfare.wavを使用しました。 http://www.ne.jp/asahi/music/myuu/wave/wave.htm

%matplotlib inline
import matplotlib.pyplot as plt
import librosa
import librosa.display
sr = 44100
music, sr_ = librosa.load("fanfare.wav", sr=sr)
print("play time : {} [sec]".format(len(music) / sr))
librosa.display.waveplot(music, sr=sr)

結果

play time : 11.679637188208616 [sec]

f:id:butch416:20190202225534p:plain

librosa.loadやlibrosa.display.waveplotのsr(sampling rate)のデフォルト値が22050なので注意してください。
何も指定しない場合はsr=22050になるので波形のx軸の上限が実際の再生時間とは異なる値になってしまいます。

IPython上で音源を再生したい場合はIPython.displayを使用します。
https://musicinformationretrieval.com/ipython_audio.html

import IPython.display as ipd
ipd.Audio(music, rate=sr)

2018-03-19

無名関数(lambda)について

無名関数(lambda)について
sqr = lambda x : x**2
sqr(2)=4が出力される。関数の引数は,で指定すれば複数用いることが出来る。

myfunc = lambda x,y : x**2+y
出力結果
myfunc(2,3)=7

参考:

www.lifewithpython.com

普通はsortやmapの引数で整理する時になどによく用いられる。
注意点はpython3ではmapがiteratorで返されるのでlistで指定して返さねければならない。

hiroto1979.hatenablog.jp

2018-03-07

import sys,osについて

Python

ゼロから作るDeep Learning(p73)を読んでいてMNISTを読み込む時に使うsysとosについて。

osはファイルのパスを指定するモジュール

それぞれの使い方は↓を参考に。

hiro-itasuto7.hatenadiary.jp

os.pardirは現在のディレクトリより一つ上の階層のディレクトリの場所を指定する「".."」に相当する。

絶対パスを指定したい場合はos.listdirを用いる。

www.yukun.info

sys.path.appendはライブラリ読み込み時のパスを追加するために用いられる。

メモ置き場

個人的備忘録(物理、機械学習、プログラミング...etc)