ヒルベルト変換のデモスクリプトをPythonで書いた話

表題の通りのまさに備忘録。

音声波形をヒルベルト変換して包絡および瞬時位相を計算し、そこから元の音声波形を再構成するスクリプトPythonで書いた話。SPTKに付属のdata.shortをwavに変換して用いた。発話内容は「青い植木鉢」である。

処理の核となるのは以下の通りである。ヒルベルト変換を計算した結果は複素数値の「解析信号」であり、絶対値を取ることで音声波形の包絡線が得られる。位相は瞬時位相と呼ばれる。これら振幅と位相の情報から元の音声波形が再構成される。

scipyのsignalパッケージを使うと上記の処理が本質的に3行で済んでしまうので、scipy最高〜〜、というわけである。

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

# wave_dataは音声データ
envelop = np.abs(signal.hilbert(wave_data))  # 包絡
angle = np.unwrap(np.angle(signal.hilbert(wave_data)))  # 瞬時位相
reconst = envelop * np.cos(angle)  # 再構成

元の波形と再構成後の波形、そして再構成誤差を図1に示す。誤差は完全に0にはならないことがわかる。

f:id:tam5917:20200512023239p:plain
図1:元の波形(上段)、再構成後の波形(中段)、再構成誤差(下段)

使ったプログラムは以下の通り。

フーリエ級数展開のデモンストレーションをPythonで書いた話

はじめに

東京大学の小山先生が、フーリエ級数展開のデモンストレーションをMATLABでお書きになった。

この素晴らしいアニメーションをPythonで再現するスクリプトを書いても良いのではないかと思い、今回の表題に至るわけである。
ちなみに再現したアニメーションは以下の通りである。グラフの軸ラベルがずっと固定であったり、描画範囲が微妙に異なるので完全再現ではないが、それなりに再現できていると思われる。

スクリプトの解説(描画まわり)

小山先生のMATLABスクリプトを参考に、以下のPythonスクリプトを書いた。
https://gist.github.com/tam17aki/4dc145b3f7d128fdb15ebf0e9137cd9b

基本的にはオリジナルのMATLABスクリプトを逐次Pythonに翻訳していく方針である。配列操作にはnumpy、またグラフの描画とアニメーション作成にはmatplotlibというライブラリを用いた。

アニメーションを作成する以前に、まず複数のグラフを一枚のキャンバスに重ね描きしないといけない。matplotlibはデフォルトで重ね描きが有効になっているので、MATLABのようにhold onとhold offを切り替える必要はない。

キャンバスを用意するにはまず以下を書く。

import matplotlib.pyplot as plt
fig = plt.figure(figsize=[8.0, 8.0])  # 800 x 800

これで図の外枠ができ上がるので、add_axes関数により実際のグラフ描画枠を配置していく。描画位置を指定するsubplotをより細かく制御する働きをしていると考えればよい。

ax1 = fig.add_axes([0.04, 0.85, 0.14, 0.14])
ax2 = fig.add_axes([0.24, 0.85, 0.74, 0.14])
ax3 = fig.add_axes([0.04, 0.65, 0.14, 0.14])
ax4 = fig.add_axes([0.24, 0.65, 0.74, 0.14])
...

add_axesの引数は「x軸の開始位置, y軸の開始位置, x軸の長さ, y軸の長さ」を、全体に対する比率で書いて与える。詳細な仕様は公式マニュアルを参考にしたほうが良い(URL)。

グラフのプロットは各描画枠(axes)ごとに行うことができる。例えば

ax2.plot(time_axis, sig, linestyle="-", color="b", linewidth=1.5)
ax4.plot(time_axis, coef[0] * numpy.sin(phi[:, 0]),
         linestyle="-", color="b", linewidth=1.5)

のような形である。

重ね描きは各axesごとにplotを並べるだけで実現できる。例えば以下のようにする。

ax3.plot(coef[0] * numpy.cos(ANGLE), coef[0] * numpy.sin(ANGLE),
         color="k", linewidth=1.5)
ax3.plot([0, circ[0, 0]], [0, circ[1, 0]],
         linestyle="-", color="b", marker="o",
         markerfacecolor="b", markersize=4)
ax3.plot([circ[0, 0], 0.6], [circ[1, 0], circ[1, 0]],
         linestyle=":", color="b", marker="o", markersize=4,
         markerfacecolor="b", linewidth=1)

以上の要素を用いることで、アニメーションの1コマに相当するグラフを描くことができるので、あとはそれらをくっつければ動画にすることができる。アニメーションの作成にはmatplotlibのanimationモジュールを用いた。

import matplotlib.animation as animation

1コマに相当する各プロットのグラフを集めて、それをリストに追加していくのが基本的な使い方である。ちなみにplotの戻り値は描画内容を格納する特殊なリストになっているので、「+=」で次々と追加していくことができる(URL)。

images = []  # アニメーションの総プロットを格納するリスト
...
for t0 in range(TIME_NUM):
...
    im = ax1.plot(coef[0] * numpy.cos(ANGLE),
                  coef[0] * numpy.sin(ANGLE),
                  color="k", linewidth=1.5)
    im += ax1.plot(coef[1] * numpy.cos(ANGLE) + circ[0, 0],
                   coef[1] * numpy.sin(ANGLE) + circ[1, 0],
                   color="k", linewidth=1.5)
...
    im += ax2.plot(time_axis, sig, linestyle="-", color="b", linewidth=1.5)
...
    images.append(im)  # 1コマ分の全プロットを集めたものをリストに追加

最終的に、以下のArtistAnimation関数を用いてアニメを作成し、save関数でmp4に書き出すことができる。

ANIME = animation.ArtistAnimation(fig, images, interval=40)
ANIME.save("plot_sawtooth.mp4", writer="ffmpeg", dpi=300)

ArtistAnimation関数の引数intervalには各フレーム間の遅延をmsecで与える。save関数の引数dpiには解像度を与える。今回はあまり高解像度にする必要はないだろう。mp4作成にあたって、事前にffmpegをインストールしておく必要がある。

おわりに

今回はmatplotlibの優秀さを実感した。アニメーション作成がやや重たいのが課題である。FuncAnimation関数を用いたとき、軽くなるかを一度検証してみる必要がある。

フーリエ級数展開はいいぞ。

navi2chの設定

事前準備

1. 2chproxy.plをダウンロードして、適切な場所に置く
https://github.com/yama-natuki/2chproxy.pl

2. bash_profileに以下を追記

/path/to/2chproxy.pl -d

navi2ch起動前に2chproxy.plをバックグラウンドで起動しておくという主旨である。バックグラウンドでの実行ができれば、他の手段で構わないと思う。

3. navi2chをインストール
melpaからインストール可能。https://melpa.org/#/navi2ch

起動

4. 以下の設定を使う

(eval-when-compile (require 'navi2ch))
(autoload 'navi2ch "navi2ch" "Navigator for 2ch for Emacs" t)
(with-eval-after-load "navi2ch"
  (setq navi2ch-list-bbstable-url "http://menu.5ch.net/bbsmenu.html")
  (setq navi2ch-list-valid-host-regexp
        (concat "\\("
                (regexp-opt '(".5ch.net" ".2ch.net"))
                "\\)\\'"))
  (setq navi2ch-net-http-proxy "localhost:8080"))

5. Emacsを立ち上げる

6. M-x navi2ch で起動する

Waveデータに対してLSB置換法に基づくステガノグラフィをPythonでやってみた

はじめに

ステガノグラフィとは、秘密のメッセージを「ばれないように、こっそりと」隠す技術である。画像メディアに対するステガノグラフィPythonパッケージは見つかるのだが、音メディア系はちょっと見当たらなかったので、試しにPythonで書いてみたということ。

今回の実装はステガノグラフィ技術の中でも一番簡単なLSB置換法に基づくものである。LSBとは最下位ビットのことを指す。音の各サンプル点における振幅値(整数)の最下位ビットは相対的に重要度が低いと考えられるので、ここに情報を埋め込む、すなわち振幅値のLSBをメッセージのビットで置き換えるわけである。

準備

1. ステガノグラフィPythonパッケージsteganoをインストール

pip install stegano

2. soundfile, numpy, scipyをインストール

3. メッセージを隠したいwavを用意

Pythonスクリプト

今回は秘密のメッセージを文字列として与えることにした(『Hello, World!』)。以下のコードではwavの入力まわりのつくりがまだ甘いが。隠したメッセージはきちんと復元できるのがポイントである。

$ python stegano_lsb.py
Hidden message is "Hello, World!"

gist.github.com

結果

(サイトからのロードの関係で初回の再生は音の冒頭が欠けて聞こえてしまう…)

  • 元の音

soundcloud.com

  • メッセージが隠された音

soundcloud.com

先頭の数十サンプルに情報を埋め込んでいるので、ぱっと聞いただけでは情報が埋め込まれているかは分からないというわけだ。

今回は先頭部分のサンプルにメッセージ情報を埋め込んだが、ランダムに埋め込む位置を変えるアルゴリズムがsteganoパッケージには実装されているので、今後はそちらのアルゴリズム実装をトライする予定である。

Pyroomacousticsを使ってfastMNMF法に基づく音源分離を試してみた

概要

ブラインド音源分離手法の1つであるfastMNMF法が件のPythonパッケージに実装されているので、手元の音源で音源分離を試してみたということ。

結果

  • オリジナルの音源信号(ドラム)

soundcloud.com

  • オリジナルの音源信号(ピアノ)

soundcloud.com

  • 混合観測信号(ドラム+ピアノ)

soundcloud.com

  • 推定された音源信号(ドラム)

soundcloud.com

  • 推定された音源信号(ピアノ)

soundcloud.com

いい感じに分離できている。素晴らしい。

ヘルダーの不等式の証明をスケッチ

ヘルダーの不等式の証明はいつも忘れるので、備忘録として証明の「スケッチ」は残しておこうという主旨である(はてなブログにおける数式の練習も兼ねて)。

1. ヘルダーの不等式とは

$1 \leq p \leq \infty$,$q$は$p$の調和共役とする.すなわち, $$ \frac{1}{p} + \frac{1}{q} = 1 $$

が成り立つものとする.3つ組 $(X, \mathcal{B}, \mu)$を測度空間とする.$f$,$g$は$X$上の可測関数とする. このとき,以下が成り立つ.

$$ ||f \cdot g||_{1} \leq ||f|| _ {p} ||g|| _ {q} $$

2. 記号の準備
  • $L _ 1$ノルム $$ ||f|| _ 1 = \int _ X\; |f| \; d\mu $$

  • $L _ p$ノルム $$ ||f|| _ p = \int _ X\; |f| ^ p \; d\mu $$

3. 証明のスケッチ

まず以下のヤングの不等式を認めることにする. $$ ab \leq \frac{a ^ p}{p} + \frac{b ^ q}{q} \;\;\; (a,\; b > 0) $$

この不等式を$\theta > 0$を用いて少し変形する. $$ ab = (\theta a) (\theta ^ {-1} b) \leq \frac{\theta ^ {a} a ^ {p}}{p} + \frac{\theta ^ {-q} b ^ {q}}{q} $$

ここで$a = |f|$および$b = |g|$として $$ |f||g| \leq \frac{\theta ^ {p} |f| ^ {p} }{p} + \frac{\theta ^ {-q} |g| ^ {q}}{q} $$ となるので、両辺を積分することで

$$ ||f \cdot g|| _ {1} \leq \frac{\theta ^ {p} ||f|| _ {p} ^ {p} }{p} + \frac{\theta ^ {-q} ||g|| _{q} ^ {q}}{q} $$

という不等式が成り立つ。そこで $\theta > 0$を $$ \theta ^ {p} ||f|| _ {p} ^ {p} = \theta ^ {-q} ||g|| _{q} ^ {q} $$

となるように選ぶ.すると不等式の右辺は$\theta ^ {p} ||f|| _ {p} ^ {p}$でくくり出すことができて $$ ||f \cdot g|| _ {1} \leq \left ( \frac{1}{p} + \frac{1}{q} \right ) \cdot \theta ^ {p} ||f|| _ {p} ^ {p} = \theta ^ {p} ||f|| _ {p} ^ {p} $$

となる.さてここで $$ \left ( \frac{||g|| _ {q}}{\theta} \right ) ^ {q} = \left ( \theta ||f|| _ {p} \right ) ^ {p} $$

なのだから,両辺の「$q$乗根」を取って

$$ \frac{||g|| _ {q}}{\theta} = \left ( \theta||f|| _ {p} \right ) ^ {\frac{p}{q}} $$

である.$q/p$は

$$ 1 + \frac{p}{q} = p \Leftrightarrow \frac{p}{q} = p -1 $$

と計算できるから,結局

$$ \begin{align} \frac{||g|| _ {q}}{\theta} &= \left ( \theta ||f|| _ {p} \right ) ^ {p -1} \\ ||g|| _ {q} &= \theta ^ {p} ||f|| _ {p} ^ {p -1} \end{align} $$ となる.よって,

$$ ||f \cdot g|| _ {1} \leq ||f|| _ {p} \cdot \theta ^ {p} ||f|| _ {p} ^ {p-1} = ||f|| _ {p} ||g|| _ {q} $$

を得て、証明が終わる.

Pythonの音声区間検出ライブラリ inaSpeechSegmenterを試してみた話

Pythonでいい感じの音声区間検出してくれるライブラリはないかなと探していたら、inaSpeechSegmenterというものが見つかったので使ってみた。
github.com

デフォルトでは音声の区間、音楽の区間、ノイズの区間、無音の区間を検出し、その区間ラベルと時間情報(開始・終了時刻)の情報を返してくれる。面白いのは音声区間においてジェンダーの区別もデフォルトでしてくれるということ。

区間ラベルを列挙しておこう:

  • speech (male/female) 音声
  • music 音楽
  • noise ノイズ
  • noEnergy 無音

インストールは以下。

pip install inaSpeechSegmenter

こんなスクリプトを書いてみた。
gist.github.com

inaSpeechSegmenterは区間ラベルと時間情報のみを取得できる。同ライブラリのseg2csvをインポートすることで、検出結果をcsvに書き出すことができる 。

seg2csv(segmentation, 'myseg.csv')

実用上は複数の音声ファイルに分割したいので、上記のスクリプトではpydubのAudioSegmentの力を借りて音声区間ごとにwavに分割保存したわけである。