音声1

覚えること

ピッチ周波数
フォルマント周波数
ケプストラム、リフタ、ケフレンシー

音声生成過程
線形予測分析
ケプストラム
課題
- 補足 1) MATLABサブルーチンの作り方
参考文献

6-1 音声生成過程

音のスペクトル解析を行った所では、振幅スペクトルと位相スペクトルについて
紹介して来ました。横軸が周波数で、縦軸が振幅や位相の成分でした。音のデータ
を特に音声データに限ると特徴的なスペクトルが現れます。特徴的に現れるスペク
トルの「音声のパラメータ」について今回は学びます。

まず音声生成の仕組みを見てみましょう。 音声成過程において、肺(lungs)、声門(glottis)は音源、声道(vocaltracts)は共振、 口唇(lips)は放射としての役割を持ちます。声門と声道(声門から口唇まで)のイメージ図は以下の通りです。声道は鼻につながる 経路がありますが図では省略してあります。

音源
音源は声門は緊張して閉じた状態になり呼吸流により断続的な振動を起こします。
男性の場合1秒間に100から150回、女声の場合250から300回ほど振動します。
この振動数をピッチ周波数、あるいは単にピッチと呼び、声の高さを表します。
共振
共振(共鳴)は物理学的には外部エネルギーによって固有振動が起こる状態を言います。
音叉（おんさ）をたたいて発する音は小さいものですが、これを箱状の共鳴器に接続
すると、驚くほど音が大きくなります。これも共振（共鳴）現象を利用したもので、
音叉の振動数に合わせたサイズの共鳴器によって音圧が著しく高まるのです。

声道を通るときに「ア」「イ」などの響きがつけられます。これは、口が形を変えて
共振の特性を変えているから「ア」や「イ」に聴えるのです。この振幅スペクトルの
ピークである共振周波数をフォルマント周波数と呼び、低い方から順に第1フォルマン
ト(F1)、第2フォルマント(F2)、第3フォルマント(F3)、、、と名付けられています
(下図参照:ただし図は振幅スペクトルの包絡線のみを描いている)。

「ア」などの母音はフォルマント周波数で特徴つけることが出来ます。次の図は
F1とF2平面に日本語の5母音を配置したものです。フォルマント周波数は性別、年齢、
話者によってかなり変動します。図よりF1-F2平面だけでは5母音を区別出来ない
ことが分かります。単独で発話される母音はF1-F2-F3空間でほぼ完全に区別出来ます。

実際のスペクトルは以下の図の様になります。これは/o/を発音しているときのスペクトルです。
青線がFFTによる振幅スペクトル。緑線がケプストラム分析による振幅スペクトルの包絡線です。
F1が400Hz付近に、F2が1100Hz付近にピークとして現れている様子が分かります。
周波数スペクトルに見る音源と共振

上で少し紹介した様に、振幅スペクトル包絡のピーク(共振部分)は母音を区別するための
特徴でした。子音の場合は共振に加えて、鼻へ通じる声道の分岐によって生じる反共振部分
(スペクトル包絡の谷)も特徴量となります。この様に振幅スペクトルにおいて、スペクトル
包絡、つまりスペクトルの空間周波数成分の低周波の部分が、声道の共振を表すと考えられ
ます。

一方、音源の特徴はスペクトルの空間周波数成分の高周波の部分に現れます。下の音声ファイルのサウンド
スペクトログラムを見て下さい。赤い色はスペクトルの振幅値が大きいことを示しています。
母音でピッチがはっきりと出るところは、周波数上で等間隔に赤い線が並んでいます。
上の図でしたら、ほぼ等間隔に出現している青線の小刻みな高周波のピークがスペクトロ
グラムでの赤い色になるわけです。このことから、音源を周波数軸上に等間隔に配置した
インパルスでモデル化する方法が一般的に使われています。
ではスペクトルの空間周波数成分の低周波(共振情報)や高周波(音源情報)をパラメータで表せないでしょうか。
音声のパラメータのうち代表的なものとして、LPC(線形予測分析)とケプストラム
があります。今回はケプストラムを中心に音声分析の具体例を見て行きます。

6-2LPC(線形予測分析)

LPCについては火曜日のゼミの本で詳しくやりますので、今回は説明を省略します。

6-3 ケプストラム

音声波x(n)は音源波g(n)と声道のインパルス応答v(n)との畳み込みとして表現出来る(放射の影響は
ないものとする)と仮定する。すなわち

x(n)=g(n)*v(n)

この信号x(n)に離散フーリエ変換の絶対値の対数演算操作Dを掛けると、対数は掛け算を和の形に
するので以下の様になる。

D{x(n)}=D{g(n)*v(n)}=D{g(n)}+D{v(n)}

(注:)積が和に置き換えられるのは対数演算のためである。

このときD{x(n)}の逆離散フーリエ変換をケプストラムc(n)と呼ぶ。対数演算によりスペクトルの
高周波と低周波が分離され、逆離散フーリエ変換を施すと時間軸の低次のところにスペクトル低周波
の共振情報が、時間軸の高次のところに音源情報が現れる。(逆離散フーリエ変換を行うので、ケプス
トラム変数は時間の次元を持つ。)

D{v(n)}はケプストラムに窓を掛けて短時間成分(つまり周波数領域で低周波成分)だけを取りだす。
これに離散フーリエ変換を行うとスペクトル包絡(共振情報)が求まる。
D{g(n)}のケプストラムのピッチ成分は、現れたり現れなかったりする。出現するときは、ピッチは
周波数軸上等間隔(調波)なので、ある高次のケプストラム次数で大きなピークとして現れる。
ケプストラムは時間領域なので、次元は時間である。が、一般に横軸を
時間(time)と呼ばず、周波数(frequency)をもじってケフレンシー(quefrency)と呼ぶ。

課題

次のMATLABプログラムを実行しなさい。またこのプログラムについて以下の質問に答えなさい。

============================================================================

フレーム表示のMOVIEを作るプログラム

波形表示 (WavFrame.m)
- 質問1 以下に挙げるフレームの条件設定を答えなさい。
  - フレーム長 (? msec)
  - フレームシフト(? msec)
  - 窓関数( 窓かけなしか、ブラックマン窓か、ハミング窓か、ハニング窓か )
FFT振幅スペクトル表示 (SpecFrame.m)
- 質問2 振幅スペクトルは線形表示とdB表示のどちらか?
ケプストラム表示 (CepFrame.m)
- 質問3 ピッチはケプストラム次数がいくらの所で現れているか? フレームにより
  変動するのでおおよその値を答えなさい。
- 質問4 ピッチがきれいに現れるのは以下のどの区間か。
  1. 無音
  2. 母音
  3. 子音
低次ケプストラムによる振幅スペクトル表示 (CepSpecFrame.m)
- 注意:参考のためFFT振幅スペクトルも表示してあります。
- 質問5 リフタの次数(低次ケプストラム係数の次数)はいくらか。
- 質問6 リフタの次数は60でもよいか? 答と共に理由ものべなさい。


============================================================================

課題が出来たら、答えを山本(eli@sys.wakayama-u.ac.jp)までメイルで送って下さい。

締め切りは12/7(木)17:00までです。

プログラム1(WavFrame.m)

コード
%===================================================================

close all;  % 画面クリア
clear all;  % データクリア

[x,fs]=aiffread('kousyouAD16.aiff'); %データ読み込み

t=[0:1/fs:length(x)/fs-1/fs]; %時間軸の設定
wlen=fs*0.032; %窓長
wshift=fs*0.008; % 窓シフト

wnum=length(x)/wshift-3; %フレーム数の計算
wnum=floor(wnum); % 整数に直す

f=figure; % 図のID取得
set(f,'doublebuffer','on');%アニメーション表示設定
M=moviein(wnum,f); %ムービーバッファ初期化

hw=hanning(wlen);


for n=0:wnum-1

  S=n*wshift; %フレーム開始点
  xtmp=x(S+1:S+wlen); % フレームデータ取得
  xtmp2=hw.*xtmp; %窓掛け
  Z(1:wlen)=  xtmp2'; % データの転置
   
  subplot(2,1,1);plot(t(1:S+wlen),x(1:S+wlen), ...
              t(S+wlen+1:end),x(S+wlen+1:end));
  xlabel('time');ylabel('Amplitude');
    
  subplot(2,1,2);plot(Z);axis([1 wlen -10^4 10^4]);
  xlabel('Data Points');ylabel('Amplitude');
        
  M(:,n+1)=getframe(f);%ムービー画像の作成
  disp(['n=' int2str(n)]); %フレーム番号表示

end

movie(f,M);%ムービーの再生

%===================================================================

画面
注意事項
- 上の波形で色が変化している所がフレーム開始点です。
- 下段図の横軸はデータポイント数ですが、上と同じく次元は時間です。

プログラム2(SpecFrame.m)

コード
%===================================================================

close all;  % 画面クリア
clear all;  % データクリア

[x,fs]=aiffread('kousyouAD16.aiff'); %データ読み込み

t=[0:1/fs:length(x)/fs-1/fs]; %時間軸の設定
wlen=fs*0.032; %窓長
wshift=fs*0.008; % 窓シフト

wnum=length(x)/wshift-3; %フレーム数の計算
wnum=floor(wnum); % 整数に直す

f=figure; % 図のID取得
set(f,'doublebuffer','on');%アニメーション表示設定
M=moviein(wnum,f); %ムービーバッファ初期化

hw=hanning(wlen);
fq=(1:wlen)*fs/wlen;  fq=fq(1:wlen/2);%周波数軸の設定


for n=0:wnum-1

  S=n*wshift; %フレーム開始点
  xtmp=x(S+1:S+wlen); % フレームデータ取得
  xtmp2=hw.*xtmp;
  
  y=fft(xtmp2');
  Z= 20*log10(abs(y(1:wlen/2))) ;
     
  subplot(2,1,1);plot(t(1:S+wlen),x(1:S+wlen), ...
              t(S+wlen+1:end),x(S+wlen+1:end));
  xlabel('time');ylabel('Amplitude');
   
  subplot(2,1,2);plot(fq,Z);axis([1 fs/2 0 120]);
  xlabel('Frequency');ylabel('Magnitude'); 
          
  M(:,n+1)=getframe(f);%ムービー画像の作成
  disp(['n=' int2str(n)]); %フレーム番号表示

end

movie(f,M);%ムービーの再生

%===================================================================

画面
注意事項

プログラム3(CepFrame.m)

コード
%===================================================================

close all;  % 画面クリア
clear all;  % データクリア

[x,fs]=aiffread('kousyouAD16.aiff'); %データ読み込み

t=[0:1/fs:length(x)/fs-1/fs]; %時間軸の設定
wlen=fs*0.032; %窓長
wshift=fs*0.008; % 窓シフト

wnum=length(x)/wshift-3; %フレーム数の計算
wnum=floor(wnum); % 整数に直す

f=figure; % 図のID取得
set(f,'doublebuffer','on');%アニメーション表示設定
M=moviein(wnum,f); %ムービーバッファ初期化

hw=hanning(wlen);


for n=0:wnum-1

  S=n*wshift; %フレーム開始点
  xtmp=x(S+1:S+wlen); % フレームデータ取得
  xtmp2=hw.*xtmp;
  
  y=fft(xtmp2');
  Z=real(ifft(log(abs(y))));% ケプストラム
   
  subplot(2,1,1);plot(t(1:S+wlen),x(1:S+wlen), ...
              t(S+wlen+1:end),x(S+wlen+1:end));
  xlabel('time');ylabel('Amplitude');
    
  subplot(2,1,2);plot(Z);axis([1 wlen -1 1]);
  xlabel('Quefrency');ylabel('Amplitude');
        
  M(:,n+1)=getframe(f);%ムービー画像の作成
  disp(['n=' int2str(n)]); %フレーム番号表示

end

movie(f,M);%ムービーの再生

%===================================================================

画面
注意事項
- 左右対照

プログラム4(CepSpec.m)

コード
%===================================================================

close all;  % 画面クリア
clear all;  % データクリア

[x,fs]=aiffread('kousyouAD16.aiff'); %データ読み込み

t=[0:1/fs:length(x)/fs-1/fs]; %時間軸の設定
wlen=fs*0.032; %窓長
wshift=fs*0.008; % 窓シフト

wnum=length(x)/wshift-3; %フレーム数の計算
wnum=floor(wnum); % 整数に直す

f=figure; % 図のID取得
set(f,'doublebuffer','on');%アニメーション表示設定
M=moviein(wnum,f); %ムービーバッファ初期化

hw=hanning(wlen);
fq=(1:wlen)*fs/wlen;  fq=fq(1:wlen/2);%周波数軸の設定

for n=0:wnum-1

  S=n*wshift; %フレーム開始点
  xtmp=x(S+1:S+wlen); % フレームデータ取得
  xtmp2=hw.*xtmp;
  
  y=fft(xtmp2');
  N=30;
  cep=real(ifft(log(abs(y))));
  lifter=[ones(1,N) zeros(1,length(cep)-2*N) ones(1,N)];%リフタの作成
  lowcep=cep.*lifter;%リフタの重みつけによる低次ケプストラムの抽出
   
  yy=fft(lowcep);

  Z1= 20*log10(abs(y(1:wlen/2))) ;
  Z2= 20*log10(exp(abs(yy(1:wlen/2))) );

   
  subplot(2,1,1);plot(t(1:S+wlen),x(1:S+wlen), ...
              t(S+wlen+1:end),x(S+wlen+1:end));
  xlabel('time');ylabel('Amplitude');
    
  subplot(2,1,2);plot(fq,Z1,fq,Z2);axis([1 fs/2 0 120]);
  xlabel('Frequency');ylabel('Magnitude');
        
  M(:,n+1)=getframe(f);%ムービー画像の作成
  disp(['n=' int2str(n)]); %フレーム番号表示

end

movie(f,M);%ムービーの再生

%===================================================================

画面
注意事項
- 青線=FFT振幅スペクトル
- 緑線=低次ケプストラムによる振幅スペクトル

参考文献

MATLAB Getting Started with MATLAB
ユーザーズディジタル信号処理

音声1

覚えること

6-1音声生成過程

6-3 ケプストラム

音声波x(n)は音源波g(n)と声道のインパルス応答v(n)との畳み込みとして表現出来る(放射の影響は ないものとする)と仮定する。すなわち

課題

プログラム1(WavFrame.m)

プログラム2(SpecFrame.m)

プログラム3(CepFrame.m)

プログラム4(CepSpec.m)

参考文献

6-1 音声生成過程

音声波x(n)は音源波g(n)と声道のインパルス応答v(n)との畳み込みとして表現出来る(放射の影響は
ないものとする)と仮定する。すなわち