MATLAB版WORLD を GNU Octave で使う

音声分析合成システム WORLD を GNU Octave 上で動かしてみる。

WORLDにはC++版とMATLAB版がある。ふと思い立って試作するには、C++ よりも MATLAB のほうがはるかにラクだ。しかし MATLAB は高価なので、個人で手を出すのは難しい。そこで MATLAB 互換のオープンソースソフトウェア Octave を使うことにした。

この記事では MATLAB版WORLD を Octave で動かすときの注意点と、WORLD で声をこねくりまわして遊ぶためのメモを書き連ねていく。

注 : この記事のとおりにOctaveで動作させると、C++版や本物のMATLABで実行したときと結果が一致しない（原因はおそらくビルトイン関数の結果が微妙に異なるため）。一応それらしい結果になるし、単に「動いたら面白いな」という程度なので、ここでは問題ないと考えて割り切っている。正しい結果が必要なときは Julia の WORLD.jl パッケージ を利用したほうがよいとおもう。

WORLD の使い方は、おもに森勢先生の公演スライドとPDF資料をベースにしている（URLを記事の最後に載せてある）。とてもわかりやすい資料なので、目を通しておくべし。

環境


WORLD を使う準備

Matlab版WORLD のZIPファイルをダウンロードして解凍しておく。

基本的に Octave は MATLAB 互換なので、普通はこのまま動作できる。
ところが、Octave には実装されていないMATLABの関数が存在する。
今回使った WORLD v0.2.1_4 (Matlab) で使われている rng関数 と interp1q関数 は、現時点では Octave に実装されていない。

というわけで、埋めあわせ用の関数を用意した。
rng.m
function rng(x)
    randn('seed', x)
    rand('seed', x)
end
interp1q.m
function yi = interp1q(x, y, xi)
    yi = interp1(x, y, xi, 'linear');
end
これら2つのファイルを作成して、WORLDと同じディレクトリに置いておく。

次に Octave を開いて、カレントディレクトリを移動する。
cd 'C:\WORLDをダウンロードしたパス'
さらに、信号処理のパッケージを読み込む。以下のコマンドで MATLAB の Signal Processing Toolbox に実装されている関数を利用できるようになる。
pkg load signal

これでようやくWORLDを動かす準備ができた。


分析と合成の流れ

WORLD で音声波形を分析する際の基本的なフローは、次のようになる。
  1. 波形データから 基本周波数 (F0) を推定 → Dio
  2. 波形データとF0から スペクトル包絡 を推定 → CheapTrick
  3. 波形データとF0から 非周期性指標 を推定 → D4C
上記のほかにも、Dioで分析した基本周波数データを補正する StoneMask や、Dioよりも高精度に基本周波数を分析できる Harvest も利用できる。

逆に、分析データから波形を合成するときは Synthesis をつかって一発でできる。
（ちなみに、C++版WORLDにはリアルタイム合成できる Synthesis2 などもある。）

Octave から以下のコマンドを順番に実行していくと、一連の分析合成を試すことができる。
[x, fs] = audioread('vaiueo2d.wav');

r = Dio(x, fs);
f = CheapTrick(x, fs, r);
q = D4C(x, fs, r);
y = Synthesis(q, f);
audiowrite('resynth.wav', y, fs);

1. WORLD付属の音声ファイル vaiueo2d.wav を読み込む
[x, fs] = audioread('vaiueo2d.wav');
2. Dioで基本周波数(F0)を推定する
r = Dio(x, fs);
※基本周波数の範囲を決めることもできる (例: 70Hz から 200Hz まで)
　推定した基本周波数が不連続になるときは、分析する周波数の範囲を変えてみるとうまくいくことがある。
op.f0_floor = 70;
op.f0_ceil = 200;
r = Dio(x, fs, op);
3. 推定した基本周波数 r を補正する（任意）
r.f0 = StoneMask(x, fs, r.temporal_positions, r.f0);
4. スペクトル包絡 f を推定する
f = CheapTrick(x, fs, r);
5. 非周期性指標 q を推定する
q = D4C(x, fs, r);
6. 再合成して波形を生成する
y = Synthesis(q, f);
7. 再合成した音声波形を resynth.wav へ書き出す
audiowrite('resynth.wav', y, fs);

Tips: 上記であつかっている変数 r, f, q は構造体になっている。Octave では fieldnames 関数を使って、構造体のフィールド一覧を取得できる。
fieldnames(r)

プロットする

WORLD で分析して得られたパラメータを可視化してみる。
以下は、すでに上記の 1.～5. までのコマンドを実行した前提で進める。


基本周波数 (F0)

figure
plot(r.temporal_positions, r.f0)
xlabel('Time [s]'); ylabel('F0 [Hz]')
基本周波数 (F0)


スペクトル包絡

スペクトログラムとして表示する。
見やすくするために、振幅の対数をとっている。
T = f.temporal_positions;
F = 0:(fs/2/size(f.spectrogram,2)):(fs/2);
S = log(f.spectrogram);
figure; imagesc(T, F, S); set(gca, 'YDir', 'normal')
xlabel('Time [s]'); ylabel('Frequency [Hz]')
スペクトル包絡 (リニア周波数軸)
このスペクトル包絡のデータからは、フォルマント周波数や近似的な声道断面積などを求められる。


非周期性指標

これもスペクトログラム的に表示する。
T = q.temporal_positions;
F = 0:(fs/2/size(q.aperiodicity,2)):(fs/2);
S = q.aperiodicity;
figure; imagesc(T, F, S); set(gca, 'YDir', 'normal')
xlabel('Time [s]'); ylabel('Frequency [Hz]')
周期性指標の時間変化
非周期性指標は 0 (周期的, 紫) ～ 1 (非周期的, 黄) の値で、各周波数帯ごとに計算される。
すべての周波数帯で 非周期性指標 = 0 ならインパルス列、非周期性指標 = 1 ならホワイトノイズをあらわすことになる。


パラメータを加工して合成

分析パラメータを編集して再合成すると、ボイスチェンジャのようなことができる。


基本周波数 (F0) を1オクターブ上げる
Synthesis(q, f) に渡す q.f0 の各要素を2倍する (=1オクターブ上げる)。
q_new = q;
q_new.f0 = q.f0 * 2;
y = Synthesis(q_new, f);
audiowrite('resynth_edit_f0.wav', y, fs);
うわずったような声になった resynth_edit_f0.wav が出力される。


話速を 1/2 倍する
Synthesis(q, f) に渡す q.temporal_positions の各要素を2倍する(=再生時間を2倍にする=再生速度を1/2倍にする)。
q_new = q;
q_new.temporal_positions = q_new.temporal_positions * 2;
y = Synthesis(q_new, f);
audiowrite('resynth_edit_speed.wav', y, fs);
アァイウゥエォ とゆっくりとした声になった resynth_edit_speed.wav が出力される。

ただ、この方法は単純にすべての発音の長さを 1/2 倍しているだけなので、あまり自然な音声は期待できない。
もっと自然に話速を変えるには、母音と子音で再生速度を変えたり、調音（音と音の移り変わり）の速度を指定できるようにしたり……といった工夫が必要。


声道長を 0.8 倍に伸縮させる
いわゆるボイスチェンジャのように、声質を変える。
線形的に伸縮させたスペクトルと、声道の長さを伸縮したときのフォルマントの移動は、だいたい一致する。つまり、スペクトルを補間すれば、声道長を変化させたような声質にすることができる。
ただし、線形振幅のまま補間すると包絡が崩れやすいので、いったん対数振幅へ変換して処理している。
f_new = f;
S = f.spectrogram;
fftsize = (size(S, 1) - 1) * 2;
w = (0:fftsize-1) * fs / fftsize;
w2 = w * 0.8;
for i = 1:size(S, 2)
    tmp = log(S(:, i));
    tmp = [tmp(:); tmp(end-1:-1:2)];
    tmp2 = interp1(w, tmp, w2, 'linear');
    f_new.spectrogram(:, i) = exp(tmp2(1:fftsize/2+1));
end
y = Synthesis(q, f_new);
audiowrite('resynth_edit_vtlen.wav', y, fs);
声質の変わった resynth_edit_vtlen.wav が出力される。


雑感

とても簡単につかえるので快適。MATLABやOctaveはサクッとゴツい処理ができる点がよく、WORLDのAPIが簡潔なことも相まって、トライアルアンドエラー的に実験を繰り返しやすい。

再合成音声の品質もよい。再合成波形と元波形を比べると、見た目は違うのに、同じように聴こえるのが不思議。また、再合成波形のエンベロープが変わることがあるが、おそらく位相の情報が欠落した結果だと思う。

少しパラメータをいじって再合成すると、ボイスチェンジャーのようなことができておもしろい。ただ、パラメータに手を加えすぎると、あまり自然でなくなってしまう。パラメータの加工をうまくやる必要がありそう。

分析・合成にはすこし時間がかかるが、個人的には許容できる範囲だと思っている。
ただしF0推定法にHarvestをつかう場合だけは、めちゃめちゃ時間がかかる。
パフォーマンスが気になる場合は、C++版やJulia版 (WORLD.jl) を使ったほうがよいかもしれない。


おまけ: Windows版Octave でうまくプロットできないとき

GNU Octave はWindowsサポートがやや雑なようで、グラフをプロットするときに発狂してフリーズしやすい。
Octave を起動したら、まずは
plot(1:10)
を実行してみる。
初回の描画には、数分かかることがある。どうやらフォントを読み込んでいるようす。うまくいけば、2回目以降のプロットは問題なく表示される。

「応答なし」になってしまう場合は、GUIツールキットを FLTK や Gnuplot に変えると（やはり最初の描画に数分かかることがあるが）多くの場合、プロットできるようになる。
graphics_toolkit('fltk')
graphics_toolkit('gnuplot')


参考

めざせ音声分析合成マスター！～「よくわからない」から「ちょっとわかる」へのチュートリアル～
WORLDの使い方に関する資料。「発表資料」からダウンロードできるPowerPointファイルがとてもわかりやすい。

Plot window not responding - Stack Overflow
Windows版Octaveのプロットウィンドウに関する不具合の対処方法。

コメント

コメントを投稿

このブログの人気の投稿

基本波形の生成

イメージ
オーディオやマルチメディア関連のプログラムを書いていると、手早くノコギリ波や矩形波を生成したいことがある。 そこで、簡単にいろいろな波形を生成する方法をまとめる。 1. ノコギリ波 1サンプルごとに少しずつ値を足していくと、直線的に増加していく波形ができる。 値が 1.0 を超えた瞬間に 0.0 へ戻すことで、ノコギリ状の周期波形が生成できる。 これは、値の小数点以下のみを取り出すことで実装できる。 実装してみる。 以下のコードでは、オーディオ信号として扱いやすいよう出力値の範囲を 0〜1 から -1〜1 へ変換している。 double sawtooth(double 周波数, double サンプルレート) { static double phase = 0.0; // 初期値 phase += 周波数 / サンプルレート; phase -= floor(phase); // 整数部分を引き算 return 2.0 * (phase - 0.5); // -1.0～1.0 の範囲に広げる } 「小数点以下を取り出す」を「1.0で割ったときの剰余」として実装する方法もある。 double another_sawtooth(double 周波数, double サンプルレート) { static double phase = 0.0; // 初期値 phase += 周波数 / サンプルレート; phase = fmod(phase, 1.0); // 小数部分を取り出す return 2.0 * (phase - 0.5); // -1.0～1.0 の範囲に広げる } プログラミング言語によっては（C# など）、浮動小数点数（float や double）に対して剰余演算子 % が使えるので、この実装方法のほうがシンプルに書ける。 2. 矩形波 ノコギリ波が 0.5 を超えたら 1.0、超えなければ -1.0 を出力する。 条件を「0.75 を超えたら 1.0 を出力」に変えれば、デューティ比25%の矩形波になる。 double pulse(double 周波数, double サンプルレート) { static double phase =
続きを読む

(1)C言語で音声合成もどき ～WAVファイルを生成する～

この記事は素人が音声合成で遊んでいるだけの記事です。完全に行き当たりばったりなので、紹介している内容の保証はできません。 また、記事の内容を予告なく変更することがあります。 「C言語で音声合成もどき」ということで、このシリーズではフォルマント合成と呼ばれる手法を使った音声合成プログラムをC言語で作っていきます。 本題の音声合成に入る前に、まずは肩慣らしとして 音声ファイルを作るプログラム をつくります。 今回は  8bit, 44100Hz, 1ch  の WAV形式 のファイルを作成するプログラムを書いてみます。 ( この記事では、比較的実装が簡単なWAV形式を扱います。 ひょっとするとMP3形式等のほうが馴染み深いかもしれませんが、MP3ではデータ圧縮にやや複雑なアルゴリズムを使うことになります。) なお、WAVファイルの構造については WAVファイルフォーマット WAVE file format を参考にしました。
続きを読む

(2)C言語で音声合成もどき ～母音の生成～

イメージ
この記事は素人が音声合成で遊んでいるだけの記事です。完全に行き当たりばったりなので、紹介している内容の保証はできません。また、記事の内容を予告なく変更することがあります。 生成した「お」の波形 今回は本題の音声合成に移っていきます。 前回 は波形データをWAVファイルに書き出す部分を作ったので、ここからはさっそく母音の合成をしていきます。 ざっくりいうと人間の声は、 声帯 で基となる音がつくられ、 声道 で加工されたものです。 今回のプログラムは、この「声帯→声道」という構成を参考にして音声を合成していきます。 まずは、声帯にあたる部分から実装することにします。 声帯が振動すると、周期的な波が発生します。 つまり、なんらかの周期的な波を使えば模倣できそうです。 一般的なシンセサイザーのようにノコギリ波や三角波を使うこともできますが、ここでは声帯振動の波を模倣した Rosenberg波 とよばれる波形を使います。 以下がRosenberg波の式です。 出典： https://dspace.jaist.ac.jp/.../896paper.pdf のP14 (2.1)、一部式の変形・表記変更あり τ1は 声門開大期 *1 の長さ 、τ2は 声門閉小期 *2 の長さ を表します。 *1 声門開大期……声門が閉じた状態から開ききるまでの時間 *2 声門閉小期…… 声門が開いた状態から閉じきるまでの時間 これらの値は、τ1 + τ2 が 0〜1 の範囲におさまるように決めます。 今回は  τ1 = 0.90, τ2 = 0.05 としました。 Cで書くとこんな感じです。 double GenRosenberg(double freq) { /* Rosenberg波を生成 */ static double t = 0; double tau = 0.90; /* 声門開大期 */ double tau2 = 0.05; /* 声門閉小期 */ double sample = 0.0; t += freq / (double)サンプリング周波数; t -= floor(t); if (t <= tau) {
続きを読む