[rogy Advent Calendar 2016] MAT少女育成計LAB #01

[rogy Advent Calendar 2016] MAT少女育成計LAB #01

信号処理音声処理
MATLAB, ケプストラム, ピッチ推定, フォルマント推定, ボイスチェンジャー

こんにちは！ところたんだポン！ この記事は，rogy Advent Calendar 2016 の18日目の記事に選ばれたポン！

(以降"ま○いく"要素はありません…)

前置き

みなさまごきげんよう．早く美少女になりたい．石油王でもいい．

来たるべき年末の諸々に向けて，みなさま進捗があったりなかったりすることだろうと思います． ちなみにわたしは任意のものに進捗がありません．

せめて，せめて一人寂しく大福にろうそくを立てるクリスマスではなく，できれば彩りがほしいものです． というわけで，MATLABで自分のおっさん声を美少女ボイスに変換するプログラムの一例を書いていきます．

注意

  • 私は信号処理が専門ではありませんし，まだ思いつきでやっている状態ので，変なことをやっている可能性が非常に高いです
  • Signal Processing Toolboxが必要です
  • 今回は使っていませんが，DSP System Toolboxを使うことでここに書かれている内容はかなり簡略化できます

なぜボイスチェンジャーか

FaceRigというものをご存知でしょうか？ 一言で言えば，自分の表情の変化をカメラに映すと，キャラクターがそれに連動して動くというソフトウェアです．

これを使うと，見た目はこんな風に美少女になれるのですが…

残念ながら声までは変換してくれません． したがって，がっかりしないためには，声の高さやその他の特徴を美少女風に変更してくれるボイスチェンジャーが必要になります．

ボイスチェンジャーはプロ用の高価なものから無償で配布されているシンプルなものまで，様々な用途・ターゲット向けのものが世に出回っています． 基本的な考え方はどれも似通っていて，かつ比較的シンプルに思えたので，昔学習した信号処理の知識を引っ張り出してきて実装してみようと考えました．

音声の分析

人の声の音，すなわち音声を処理して別の人が発したかのような音声信号を作ることになるわけですが，それにあたって音声を構成する要素について理解しておく必要があります． といっても専門でも何でもない私が解説しても無駄かと思いますので，参考資料を示しておきます．

東京大学『音響音声学』講義資料

全て目を通したわけではないのですが，こちらの資料は音・音声にまつわるあらゆるメカニズムを図入りで分かりやすく解説されている日本語資料として非常に参考になるのではないかと思います． 今回必要となる音声に関する情報は主に2つです．

  1. 基本周波数(ピッチ)
  2. フォルマント周波数

試しに，今回サンプルとする「にゃーん」という音声データを使ってこれらの情報を見ていこうと思います．

使用するデータは こちら に置いておきますので，皆さんも試してみてください．

まず，音の周波数成分の時系列変化を見るため，スペクトログラムというプロットを作成してみます．

% 音声データの読み込み
[Y, Fs] = audioread('nya-n.wav');

% ダウンサンプリング
Yd = resample(Y, 1, 5);
Fsd = Fs/5;

% スペクトログラムの描画
nfft = 2048; nwindow = 256; nshift = 128;
spectrogram(Yd, nwindow, nshift, nfft, Fsd, 'yaxis');

をMATLABで実行すれば同じプロットが得られます． こちらのスペクトログラムの横軸は時間，縦軸は周波数，パワーが強い周波数帯で明るい色が表示されています(パワーそのものではなくスペクトル密度ですが)．

見ていただくとわかるように，きれいな縞模様が出ていますね． これは，人間が感じる声の高さに相当する基本周波数(以降簡単に「ピッチ」と表現)の倍音が縦に並んでいるためです． 皆さんよくご存知の通り，物体が振動するときにいわゆる単振動をすることは稀で，通常はその整数倍の周波数での振動が起こります． 人間の声の場合，声帯で基本周波数の音を出すことになるのですが，そこで発生する「ブーブー」というラッパのような音に含まれる倍音の成分が縞模様としてスペクトログラムにあらわれているというわけです．

しかし，我々が聴く声というものはラッパの音とは程遠いものです．つまり，声は声帯だけで出しているものではないのです． ではどのようにしているかと言うと，骨格や口の内部の形状を変化させることによって，声を生み出しているわけです． 音が通過する空間の形が変わると何が起こるかというと，当然音の共鳴の仕方が変化します． ここでもう一度スペクトログラムを見てみましょう．

よく見ると縞模様は一様ではなくて，ある時間変化する周波数帯で強調されていることがわかります． これが，母音や声質の違いを生み出す鍵になります．このようにして現れる時間とともに変化するピークのことをフォルマントと呼びます． 男性と女性が同じ基本周波数で声を出したとしても，一般的にフォルマント周波数の値に差が出ます． この差を利用することで，男声を女声に，女声を男声に変換することができるようになります．

まじめに音声処理をしようと思うともっとたくさんのことを考慮する必要がありますが，今回は簡単にピッチとフォルマント(第1,第2のみ)を音声データから抽出し，最終的な声の合成に使用するということを考えます． この記事では，以下の図のようにピッチ(F0)とフォルマント(F1,F2)を推定して時系列上に可視化することを目標にします．

ピッチ推定

先ほどのスペクトログラムで見たように，基本周波数は倍音を伴ってスペクトルにあらわれるため，周波数に対して周期的な信号として捉えることができます． したがって，ピッチを推定するためにはその周期的信号の周期を求めてあげればいいことになりますよね！ という考え方を用いた方法が，ケプストラムを用いる手法です．

ケプストラムは，もとの音声信号をフーリエ変換して得られた周波数特性をさらにフーリエ変換(実用上は逆フーリエ変換を施す)したもののことを指します． 今回は，ケプストラムのピークを取ることによってピッチの推定を行います．

便利なことに，MATLABにはケプストラムを計算してくれるコマンドがあります(演算自体そんなに複雑ではないのでベタ書きしてもいいのですが…)．

%% 前処理
% てきとーなタイミングの信号を切り出す
ts = 8000;
y1 = y(ts:ts+nwindow-1);

% 窓関数をかける
w = window(@hamming, nwindow);
yw = y1.*w;

%% (実数)ケプストラム算出
c = rceps(yw);

plot(c(1:end/2)); % ケプストラムは折り返しがあるので半分だけプロット

プロットは以下のようになります．

ここで，ケプストラムにきれいにピークが立っているのが見えます．これがピッチに相当する部分です． 位置を見てみると，61サンプルごとにあらわれると出ていますので，サンプリング周波数が 44100/5 = 8820 Hz なので，8820/61 = 144.6 Hzと計算できます． スペクトログラムに戻って値を比較してみても正しそうだということが分かります．

これを繰り返すことによってピッチの時系列データを得ることができますが，プログラムを書くだけになるので割愛します．

注意点を挙げるとすれば，ケプストラムには音の大きさ等の情報が表れてこないため，発声しているかどうかを別途判定しておく必要があるということぐらいでしょうか．

フォルマント推定

同様の考え方でケプストラムを用いてフォルマント推定を行うこともできるのですが，今回はLPC(線形予測法)を使用したフォルマント推定を行います． なぜかといいますと，MATLABのドキュメントにサンプルが乗っかっているからです(すごい)．

こちらのページでは，(たぶんSignal Processing Toolboxの生みの親であるLorenさんの声で)"MATLAB"と発音した音声データをサンプルとして，第一音節の母音 /æ/ のフォルマントを推定しています． 同様の方針でフォルマント推定を行いましょう．

線形予測法は，未来の値を過去の値の重み付け線形和として予測するアルゴリズムです． これに信号を通して重み係数を適応的に更新させることで，ある意味で学習されたフィルタを得ることができます． そのフィルタの極なりなんなりを調べることによって，もとの信号の特性を得ることができるわけです．

詳細なアルゴリズムについては他の記事に譲るとして，MATLABでお手軽にやってみましょう．なんとLPCがコマンド1つで実行できます．

% プリエンファシスフィルタ: 高周波域を強調
preemph = [1 0.97];
yp = filter(1,preemph,y1);

% 窓関数をかける
w = window(@hamming, nwindow);
yw = yp.*w;

% LPC (1コマンドでできちゃう…)
A = lpc(yw, 8);
rts = roots(A); % フィルタの極を求める

% 共役根を除いて極に対応する角度を計算
rts = rts(imag(rts)>=0);
angz = atan2(imag(rts),real(rts));

% 極の角度から周波数を計算して昇順にソート
[frqs,indices] = sort(angz.*(Fsd/(2*pi)));
% バンド幅の計算
bw = -1/2*(Fsd/(2*pi))*log(abs(rts(indices)));

nn = 1;
for kk = 1:length(frqs)
    % フォルマント周波数は90Hzより大きく，バンド幅は400Hz未満．絞り込む
    if (frqs(kk) > 90 && bw(kk) < 400)
        formants(nn) = frqs(kk);
        nn = nn+1;
    end
end
formants

実行結果:

formants =

   1.0e+03 *

    0.5245    1.0946    2.9470

よさそうですね．これも同様に時系列で出すようにする必要がありますが，割愛します． ただし，このまま時系列に対して処理をしようとすると，フォルマントの値がギザギザしてしまったり，突然変な値に飛んだりしてしまいますので，何かしらの処理が必要になります． 無声音の部分でも推定値が乱れますので，もう少し応用的なアルゴリズムを使用するか何かした方がいいかもしれません．

適当な補間処理と無声音・無音に対する処理を加えると，以下のように処理結果が得られます．

おまけ・次回予告

まず，推定したピッチの時系列をもとにノコギリ波で声帯の音を近似してみました．

つぎに，推定したフォルマントでピークが立つようなフィルタを iirpeak コマンドで設計して，この擬似声帯音にかけてみます．

なんとなく「にゃーん」と言っているような気がしますね． では，ピッチを2倍，フォルマントを1.5倍の周波数としてこれを生成するとどうなるでしょうか．

お，それっぽい(美少女には程遠いですが)！でも音が汚くて残念ですね… フィルタリングの過程で高周波も削ぎ落とされてしまっているようです．

この方法に限界を感じて調べていたところ，TD-PSOLA法というシンプルかつ有効なアルゴリズムがあるらしいということを知りましたので，次回はそれを試してみたいと思います． 美少女への道のりは険しい．

To be continued…

Share Comments