音響解析関連

Last Change: 19-Aug-2015.
author : qh73xe

とりあえず,julia に手を出した理由は 音響解析を手軽にできそうだったからです.

目標としては wav を読み込み,可視化し,特徴量の抽出を行います.

WAV

WAV ライブラリはその名の通り, wav ファイルを扱うためのライブラリです, このライブラリを導入するには以下のようにします.

install WAV
Pkg.add("WAV")

ライブラリの使用は python というよりも R に似ています. とりあえず,wav ファイルを読み込み,再生するコードを記述します.

wav ファイルを読み込み,再生
using WAV
x, fs = wavread("hogehoge.wav")
wavplay(x, fs)

読み込んだ波形をとりあえず可視化します. 可視化には PyPlot というライブラリを使用しました.

wav ファイルの可視化
# Pkg.add("PyPlot")  # ライブラリのインストール
using PyPlot
plot(0:1/fs:(length(x)-1)/fs, x)
xlabel("Time [s]")

WORLD

音響情報の基本的な値を求めるには WORLD というライブラリを 利用するのが便利そうです(というかこれが使いたくて julia 触った.だっていちいち コンパイルとか...).

WORLD はもともと C で書かれた音響解析 & 合成ツールでSTRAIGHTの高速版という認識です.

  • 音声分析合成システム 「WORLD」
    • 最新版(0.2.0)で実装したD4Cにより,全ての条件においてTANDEM-STRAIGHTを上回る上位 互換になったそうです.
    • このツールを使用すると,音声ファイルの基本周波数,スペクトル包絡,非周期性指標を分析,変換,合成が可能になります.

このライブラリはその WORLD の julia 用ラッパーで r9y9 さんによって作成されました.

  • r9y9/WORLD.jl
    • 私のように少し音響情報がほしいだけの人間には,すごくありがたいです.

詳しくは公式のデモを確認するのがよいと思いますが, とりあえず,簡単な使い方をメモしておきます.

デモ環境の準備

とりあえず,デモ用の音声を julia に読み込みます. 実際に手持ちのファイルを使用する場合は,filepath にファイルまでのパスを書けば 動くかと思います.

using PyCallmatplotlib = pyimport("matplotlib")
PyDict(matplotlib["rcParams"])["figure.figsize"] = (12, 5)
using PyPlot
using WAV
using WORLD

# データの読み込み(デモデータを使用します)
filepath = joinpath(Pkg.dir("WORLD"), "test" , "data", "test16k.wav")
x, fs = wavread(filepath)
wavplay(x, fs)  # どんな音なのか確認(ATR だった)

これで, x, fs に音声ファイルの情報が入っています.

F0 の抽出

とりあえず,読み込んだ wav から F0 を抽出します. WORLD に音声ファイルの情報を入れる場合, データ型の変換が必要みたいなので注意です.

F0
# データ型の調整
x = vec(x)
fs = int(fs)

# 各種設定
period = 5.0
opt = DioOption(
    f0floor=40.0,
    f0ceil=700.0,
    channels_in_octave=2.0,
    period=period,
    speed=4
)
f0, timeaxis = dio(x, fs, opt)  # F0 の抽出

# 結果の可視化
figure(figsize=(16, 6), dpi=80, facecolor="w", edgecolor="k")
plot(timeaxis, f0, label=" F0 trajectory estimated by DIO", linewidth="2")
xlabel("Time [sec]")
ylabel("Frequency  [Hz]")
legend()

注釈

DioOption の引数に関して

DioOption の引数に関していくつか不明なものがあるので 調査が必要.

  • f0floor : 解析する F0 の一番低い値
  • f0ceil : 解析する F0 の一番高い値
  • period : 解析フレームの大きさ(msec)
  • channels_in_octave: 不明
  • speed : 不明

スペクトラム包絡

F0 の値を利用してスペクトラム包絡を計算します.

Spectral envelope estimatinon
# Spectral envelope estimatino by CheapTrick
spectrogram = cheaptrick(x, fs, timeaxis, f0)
# for plot
imshow(10log10(spectrogram), origin="lower", aspect="auto")
colorbar()

非周期性の推定

最後は非周期性指標です.

Aperiodicity ratio estimation
aperiodicity = d4c(x, fs, timeaxis, f0)
imshow(20log10(aperiodicity), origin="lower", aspect="auto")
colorbar()

再合成

上記の3つが WORLD で算出できる基本的な音響情報になります. 続いては, 上記3つの値を利用して再合成を行います.

Sysnthesis from f0, spectral envelope and aperiodicity
y = synthesis(f0, spectrogram, aperiodicity, period, fs, length(x))
wavplay(y, fs)

これで,今まで解析を行った音声を再合成して音を鳴らすことができるかと思います. もちろん,このままでも面白いのですが,いくつか単純な変換例を追記します. とりしま,わかりやすくF0を対象に.

f0 の値を2倍にする
f0_b = f0 * 2
y_b = synthesis(f0_b, spectrogram, aperiodicity, period, fs, length(x))
wavplay(y_b, fs)
f0 の値を均一にする
f0_c = []
for n in [1:length(f0)]
    if f0[n] == 0
        f0_c = [f0_c, f0[n]]  # もともと0の部分は0のママにしないとおかしい
    else
        f0_c = [f0_c, mean(f0)]  # それ以外の部分はとりあえず平均値に
    end
end
y_c = synthesis(f0_c, spectrogram, aperiodicity, period, fs, length(x))
wavplay(y_c, fs)

ここまでで,とりあえず julia/WORLD の感覚はつかめるかと思います. 音声変換, 合成の正式な(ホットな??)変換はまた,後日.