Miyazawa’s Pukiwiki
MATLAB Note/音声の分析
はすでに存在します。
開始行:
#access
#analog
-目次
#contents
**音声を取り込む・再生する [#x740141a]
-%%WAVREAD Microsoft WAVE (".wav") サウンドファイルの読み...
#geshi(matlab){{
filename = 'a.wav'; %読み込むファイル...
[data,Fs,Bits] = wavread(filename); %dataに音声データ...
sound(data,Fs); %サンプリング周波数 Fs で再生
wavplay(data,Fs); %サンプリング周波数 Fs でWindo...
wavplay(data,Fs,'async'); %音声の再生に平行して処理を続...
disp('再生中...');
}}
-AUDIOREAD Microsoft WAVE (".wav") サウンドファイルの読み...
#geshi(matlab){{
filename = 'a.wav'; %読み込むファイル...
[data,Fs] = audioread(filename); %dataに音声データ、Fs...
sound(data,Fs); %サンプリング周波数 Fs で再生
player = audioplayer(data,Fs);
play(player); %サンプリング周波数 Fs でWindows ...
}}
-デフォルトで、以下のような音声データが用意されています。
#geshi(matlab){{
load laughter; %笑い声
sound(y, Fs);
load handel; %handelのハレルヤ
sound(y, Fs);
load mtlb; %MATLABと発音した英語音声
sound(mtlb,Fs);
}}
**波形をプロットする [#d9cc99ea]
-取り込んだ波形をプロットするには
#geshi(matlab){{
load mtlb; %MATLABと発音した音声をロード（変数mt...
plot(mtlb);
}}
-44Hzの正弦波を作って、波形をプロットするには
#geshi(matlab){{
Fs = 8192; %サンプリング周波数 8192[Hz]
f = 44; %周波数（音の高さ）44[Hz]
A = 1.0; %振幅（音の大きさ）1.0
length = 0.5; %音の再生時間 0.5[秒]
timeline = 0 : 1/Fs : length; %時間軸の作成 1/Fs 間隔で...
data = A * sin(2 * pi * f * timeline); %正弦波信号を生成
subplot(2,1,1); plot(data, '.') %ドットでプロット
subplot(2,1,2); plot(data) %通常のプロット
title('時間信号'), xlabel('time'); %グラフのラベル
}}
**スペクトログラム [#g39862ab]
-スペクトログラムは、パワースペクトル解析の結果を時間軸上...
-「MATLAB」と発音したサンプル英語音声を取り込んでプロット...
#geshi(matlab){{
load mtlb; %MATLABと発音した英語音声
disp(strcat('サンプリング周波数は', int2str(Fs), '[Hz] ...
sound(mtlb, Fs) %再生
wavwrite(mtlb, Fs, 32, 'mtlb.wav'); %量子化ビット数32で...
subplot(2,1,1);
plot(mtlb); xlabel('Time (sample)'); xlim([1 length(mtlb...
subplot(2,1,2);
spectrogram(mtlb, hamming(64), 32, 256, Fs, 'yaxis'); %...
% 64 サンプル（Wavesurferと同じ）の対称なHamming窓で分割...
% 32 サンプルのオーバーラップ
% 256 サンプルのFFT長
}}
--実行結果（音声波形とスペクトログラム）~
#ref(mtlb_mfcc_1.jpg)
---元の音声のサンプリング周波数と量子化ビット数によって、...
-マイクから2.0秒間音を取り込んで、取り込んだ音声のスペク...
#geshi(matlab){{
Fs = 44100; time = 2.0;
data = wavrecord(time*Fs, Fs, 1); %1...
spectrogram(data, hamming(64), 32, 256, Fs, 'yaxis'); %...
% 64 点（Wavesurferと同じ）のHamming窓で解析する、細かい...
% mtlb の各セグメントが 32 点のオーバーラップをする
% 離散フーリエ変換を計算するために使われる周波数点数は 256
% yaxis ： y軸上に周波数をプロットする
}}
-スペクトログラムからある瞬間（1、50、100 フレーム）にお...
#geshi(matlab){{
load mtlb; %MATLABにはじめからセットされている英語音声...
subplot(3,2,[1 2 3 4]);
spectrogram(mtlb, hamming(64), 32, 256, Fs, 'yaxis');
[s,f,t,p] = spectrogram(mtlb, hamming(64), 32, 256, Fs, '...
% f : スペクトログラムにおける周波数ベクトル
% t : 時間ベクトル
% p : 各セグメントのパワースペクトル密度
subplot(3,2,[5 6]);
plot([1:length(p)], p(:,1), [1:length(p)], p(:,50), [1:le...
legend('segment = 1', 'segment = 50', 'segment = 100');
xlabel('周波数'); ylabel('パワー');
}}
**音声波形のフレーム分割および前処理 [#gcaedeef]
-音声解析の前準備として、音声の高域強調・フレーム分割・窓...
***高域強調（プリエンファシス） [#q962d911]
-音声信号は低周波数成分が大きく、周波数が大きくなるにつれ...
-声帯の基本振動、空気中の伝搬の影響を除き、口腔内の形状に...
-プリエンファシスを行うことで、第二・第三フォルマントのよ...
-具体的には、ハイパスフィルタを使い、高周波数成分の振幅を...
--式は以下のとおり。~
#mimetex(y(n) = x(n) - px(n-1))
---x(n)はサンプル数nを変数とするデジタルサウンドデータ（...
---pはプリエンファシス係数。音声認識では0.97を使うことが...
---１サンプル前の値に対して、大きく値が変化した場合（高周...
--MATLAB では以下のようにする。
#geshi(matlab){{
load mtlb;
data = mtlb;
pre_emphasis = 0.97; %プリエンファシス係数
newData = [];
for countData = 2 : 1 : length(data)
%プリエンファシス
thisData = data(countData) - ( pre_emphasis * data(c...
newData = [newData ; thisData ];
end
subplot(2,2,1); plot(data); xlim([0,4500]); title('元の...
subplot(2,2,2); spectrogram(data, hamming(64), 32, 256, ...
subplot(2,2,3); plot(newData); xlim([0,4500]); title('プ...
subplot(2,2,4); spectrogram(newData, hamming(64), 32, 25...
wavwrite(newData, Fs, 'newmtlb.wav');
}}
---以下の関数でも同じ処理ができる。
#geshi(matlab){{
newData = filter([1 -0.97],1,data);
}}
--プリエンファシスの結果は以下のようになる。
#ref(mtlb_mfcc_5.jpg)
---基本周波数成分（スペクトログラム下部のボイスバーの部分...
***フレーム切り出し [#o5c30647]
-連続音声の分析を行う場合は、音声信号を短い時間単位で切り...
--1フレーム内に2周期程度の音声信号が含まれていれば、うま...
--通常フレーム同士は重複するようにとる。重複部分の長さを...
-例えば、お題の音声を、音声認識で一般的に用いられる フレ...
#ref(mtlb_mfcc_2.jpg)
--以下のコードで、フレームの切り出しを行う（&color(red){2...
#geshi(matlab){{
frameSize = 0.025; % フレーム長：0.025秒（...
frameShift = 0.010; % フレームシフト長：0.0...
load mtlb;
data = mtlb;
data = filter([1 -0.97],1,data); %プリエンファシス
frameSizeSample = fix( Fs * frameSize ); % フレー...
frameShiftSample = fix( Fs * frameShift ); % フレー...
disp( strcat('サウンドデータのサンプル数 ',int2str(lengt...
'フレーム長 ',int2str(frameSizeSample),'[sample], ',...
'フレームシフト長 ',int2str(frameShiftSample),'[samp...
%フレームの総数を求める
maxFrame = ...
fix( ( length(data) - (frameSizeSample - frameShiftSamp...
%フレームの個数だけ処理を繰り返す
startThisFrame = 1; % ...
endThisFrame = startThisFrame + frameSizeSample - 1; % ...
for countFrame = 1 : 1 : maxFrame
disp(strcat('フレーム番号 ',int2str(countFrame),' / ...
disp(strcat('このフレームの開始サンプル番号 ',num2st...
disp(strcat('このフレームの終了サンプル番号 ',num2st...
thisData = data(startThisFrame : endThisFrame);
% ここでこのフレームに対する計算を行う
subplot(6, 9, countFrame); plot(thisData); ylim([-3,...
axis off; %目盛りを消す
startThisFrame = startThisFrame + frameShiftSample;
endThisFrame = startThisFrame + frameSizeSample - 1;
end
}}
---結果は以下のようになる。51個のフレームに切り出された。
#ref(mtlb_mfcc_3.png,,80%)
-各フレームを結合して元の音声を復元する方法は [[MATLAB No...
***窓関数による丸め [#kf8e24d7]
-上で取り出したフレームは、単純に音声波形を切り出しただけ...
--そのため、上で切り出した音声に窓関数をかけあわせて、切...
-ここでは、[[ハミング窓:http://ja.wikipedia.org/wiki/%E7%...
--ハミング窓は以下の式で表わされる((Wikipedia より引用。...
#mimetex(w(x)=0.54-0.46cos2\pi x ,\quad if \quad 0 \leq ...
---MATLABでは以下の関数を使って作ることができる。
#geshi(matlab){{
frameSizeSample = 185; % フレーム長：サンプ...
window = hamming(frameSizeSample); % ハミング窓を生成
plot(window);
}}
--前述のスクリプトにハミング窓をかけてみる。
---[[mtlb_mfcc_1.m:http://shower.human.waseda.ac.jp/~m-ko...
---結果は以下のようになる。
#ref(mtlb_mfcc_4.png,,80%);
---各フレームの角がまるまっている。
-各フレームを結合して元の音声を復元する方法は [[MATLAB No...
-音声の両端にのみ窓かけをする方法は [[MATLAB Note/音声の...
**パワースペクトル解析 [#u704d8b6]
-パワースペクトルの詳細は [[メル周波数ケプストラム（MFCC...
-音声のパワースペクトル解析の手順は [[メル周波数ケプスト...
***パワースペクトルによるホワイトノイズの検出 [#a467125b]
-音声のパワースペクトルは低周波数成分のパワーが大きく、白...
#ref(analyzeNoize.m);
**線形予測法（LPC）とフォルマント解析 [#y6a12906]
-FFTによって出力されたスペクトル解析の結果をみると、小さ...
--スペクトル包絡の山は、音声の声道に由来する周波数成分で...
--音声では、スペクトル包絡の山の頂点がフォルマント周波数...
-スペクトル包絡をより明確に求めたい場合、線形予測法（LPC...
--&ref(a.wav); を使って以下のコードを実行
#geshi(matlab){{
[wav, Fs] = wavread('a.wav');
%中央部分を切り出し
center = fix(length(wav) / 2); %中心のサンプル番号
cuttime = 0.04; %切り出す長さ[s]
wavdata = wav(center-fix(cuttime/2*Fs) : center+fix(cutti...
%ハニング窓をかける
han_window = 0.5 - 0.5 * cos(2 * pi * [0 : 1/length(wavda...
wavdata = han_window(1:length(wavdata))' .* wavdata;
lpcsize = 32; %LPC次数
fftsize = 2048; %フーリエ変換の次数, 周波数ポイントの数
[P,f] = pyulear(wavdata, lpcsize, fftsize, Fs); %LPC分析
AP = abs(P) / fftsize; %振幅ス...
%プロット
fscale = linspace(0, Fs, fftsize); %0～Fsまでをfftsize...
%プロット（表示するのはサンプリング周波数の半分まで）
subplot(2,1,1); plot(fscale(1:fftsize/2), AP(1:fftsize/2)...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル（LPC）'); xl...
dft = fft(wavdata, fftsize); %比較用...
Adft = abs(dft) / fftsize; %振幅ス...
subplot(2,1,2); plot(fscale(1:fftsize/2), Adft(1:fftsize/...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル（FFT）'); xl...
}}
---[[LPC次数の決定方法（AICによる）:http://www.cs.takusho...
--結果
#ref(lpc_01.png);
---スペクトル包絡が取り出せました。
-スペクトル包絡のピークピッキングによって、フォルマント周...
--上記のコードに続いて、以下を実行します。
#geshi(matlab){{
peakNum = []; %ピークのフレーム番号
for peakcount = 2 : 1 : (length(AP)-1)
if ((AP(peakcount) >= AP(peakcount-1)) & (AP(peakcount)...
peakNum(peakcount) = 1; % ピークのところだけ1...
end
end
peakposition = find(peakNum == 1); % ピークの位置を決定
fmag = AP(peakposition); % 振幅を求める
frequences = f(peakposition); % フォルマント周波数...
if (length(frequences) >= 1)
F1 = frequences(1); % 第一フォルマント
end
if (length(frequences) >= 2)
F2 = frequences(2); % 第二フォルマント
end
if (length(frequences) >= 3)
F3 = frequences(3); % 第三フォルマント
end
disp(strcat('F1 : ', num2str(F1), ' F2 : ', num2str(F2), ...
}}
--結果
F1 :473.7305 F2 :1098.1934 F3 :3337.6465
---小さな値から第一フォルマント、第二フォルマント、第三フ...
---なお、ここで得られる値は確実にフォルマントを取り出せて...
--導関数（微分法）を使ったピークピッキングの方法について...
-文献
--[[応用音響学: 線形予測分析(1) LPC分析:http://ocw.u-toky...
**基本周波数解析((参考文献：[[MATLABマルチメディア信号処...
-ケプストラム法による方法
--長所 ： 正確 短所 ： ノイズに弱い
--ケプストラム法の詳細は [[メル周波数ケプストラム（MFCC）...
--&ref(a.wav); を使って以下のコードを実行
#geshi(matlab){{
[wav, Fs] = wavread('a.wav');
t = 0 : 1/Fs : length(wav)/Fs;
%中央部分を切り出し
center = fix(length(wav) / 2); %中心のサンプル番号
cuttime = 0.04; %切り出す長さ[s]
wavdata = wav(center-fix(cuttime/2*Fs) : center+fix(cutti...
time = t(center-fix(cuttime/2*Fs) : center+fix(cuttime/2*...
%ハニング窓をかける
han_window = 0.5 - 0.5 * cos(2 * pi * [0 : 1/length(wavda...
wavdata = han_window(1:length(wavdata))' .* wavdata;
fftsize = 2048; % フーリエ変換の次...
dft = fft(wavdata, fftsize); % 離散フーリエ変換
Adft = abs(dft); % 振幅スペクトル
Adft_log = log10(abs(dft)); % 対数振幅スペクトル
cps = real(ifft( Adft_log )); % ケプストラム
quefrency = linspace(0, cuttime*2, cuttime*Fs*2); %ケフレ...
subplot(2,1,1); plot(quefrency(1:length(cps))*1000, cps);...
%基本周波数ピーク抽出
cps(1:length(cps)/20) = 0; % ケプストラムの低...
cps(fftsize/2:fftsize) = 0; % ケプストラムの後...
subplot(2,1,2); plot(quefrency(1:length(cps))*1000, cps);...
xlabel('ケフレンシー[ms]');
}}
---結果
#ref(cps_001.png,,70%);
---上図の低次成分を除去した結果が下図です。基本周波数のピ...
--ピークのケフレンシーの値を求め、周波数[Hz]に変換します。
#geshi(matlab){{
peakQuefrency = quefrency( find(cps == max(cps)) );
f0 = 1 / peakQuefrency;
disp(f0);
}}
---結果
185.2416
---Praatで a.wav のピッチを解析した結果 &ref(a_pitch_praa...
---より正確に求めるには、閾値などを決めてエラー検出をする...
--よりノイズに強いケプストラム法として、'''[[改良ケプスト...
-[[自己相関:http://shower.human.waseda.ac.jp/~m-kouki/puk...
--長所 ： ノイズに強い（ロバスト） 短所 ： ケプストラム...
-'''[[STRAIGHTによる方法（限定公開）:http://shower.human....
-MATLAB以外にも、[[wavesurfer:http://www.f.waseda.jp/kiku...
**メルケプストラム解析 [#gab6732b]
-理論的な詳細は &pgid(,メル周波数ケプストラム（MFCC）); ...
-[[MFCC解析のツール/MATLAB Auditory Toolbox, Praat, HTK ...
-MATLABでメルケプストラム解析を実行するには、Malcolm Slan...
--[[Auditory Toolbox ダウンロード:http://cobweb.ecn.purdu...
--Auditory Toolbox の mfcc.m で解析します。(([[Auditory T...
--ためしに、「MATLAB」と発音したファイルを解析してみまし...
#geshi(matlab){{
%TEST_MFCC Auditory Toolbox によるメルケプストラム解析
load mtlb; %MATLABと発音した英語音声
sound(mtlb, Fs) %再生
%Auditory Toolbox によるMFCC解析
[ceps,freqresp,fb,fbrecon,freqrecon] = mfcc(mtlb, Fs, 10...
mfcc_data = ceps(2:13,:) %MFCC（C0値を除...
figure(1),subplot(2,2,1),imagesc(flipud(freqresp)),title...
figure(1),subplot(2,2,2),imagesc(flipud(fb)),title('fb')
figure(1),subplot(2,2,3),imagesc(flipud(ceps)),title('ce...
figure(1),subplot(2,2,4),imagesc(flipud(ceps(2:13,:))),t...
%ファイル書き出し
csvwrite('mtlb.txt', ceps);
}}
--結果は以下のようになります。~
&ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/mfc...
-【補足1】mfcc.m でのフレーム長・フレームシフト長について
--フレーム長は定数値を与えており、mfcc.m の40行目「window...
Fsは7418（1秒間に7418サンプル）なので、フレーム長は256/74...
--フレームシフト長は mfcc.m の引数である %%samplingRate,%...
---&color(red){frameRate は「1秒間をいくつのフレームに分...
---mfcc.m の111行目「windowStep = samplingRate/frameRate;...
---上の例では、frameRateは100を指定していますから、「1秒...
--mfcc.m の112行目「cols = fix((length(input)-windowSize)...
size(mfcc_data) = [12 50] になり、12次元のMFCCが50フレー...
-【補足2】フレーム長 0.025秒、フレームシフト長 0.010秒 で...
--WindowSize の値を Fs * 0.025 に書き換えます。
--フレームシフト長は1/100秒（1秒を100のフレームに分割）な...
**フィルタバンク解析 [#tb4d530b]
-'''[[デジタルフィルタ:http://shower.human.waseda.ac.jp/~...
--'''[[低域通過（ローパス）フィルタ:http://shower.human.w...
--'''[[高域通過（ハイパス）フィルタ:http://shower.human.w...
--'''[[帯域通過（バンドパス）フィルタ:http://shower.human...
---'''[[任意の形状のゲイン特性を指定したバンドパスフィル...
--'''[[帯域阻止（バンドストップ）フィルタ:http://shower.h...
-フィルタバンク解析とは
--フィルタバンク：バンドパスフィルタの集合体。入力信号を...
--重要でない帯域のサンプリング周波数を減らし、再び結合す...
--各帯域のスペクトルの積分を求めることで、スペクトルの次...
-以下のフィルタバンク解析において、バンドパスフィルタの設...
***2分割フィルタバンク((バンドパスフィルタの設計手法は [[...
-以下のようにします。
#geshi(matlab){{
load handel;
n = 150; % 150次のローパ...
bandpassWidth = [1 2000 ; 2001 4000]; % 1～2000Hzと、2...
%順番に処理する
subplotcount = 1;
for count = 1 : 1 : length(bandpassWidth)
% インパルス応答を求める
Wn = [ bandpassWidth(count,1)/(Fs/2) bandpassWidth(c...
b = fir1(n,Wn,'bandpass');
% ゲイン特性を求める
[h,w] = freqz(b,1,512);
subplot(length(bandpassWidth),3,subplotcount);
plot((w/pi)*(Fs/2),abs(h)); xlabel('Hz'); legend('ゲ...
subplotcount = subplotcount + 1;
% フィルタをかける
y2 = filter(b,1,y); % フィルタをかけ...
subplot(length(bandpassWidth),3,[subplotcount subplo...
specgram(y2,512,Fs,kaiser(500,5),475);
subplotcount = subplotcount + 2;
end
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
***オクターブフィルタバンク（オクターブバンド）((バンドパ...
-[[オクターブスケールとは？:http://shower.human.waseda.ac...
-オクターブフィルタバンクとは？
--[[JIS規格による、1/1、1/3オクターブバンドの上下遮断周波...
---1/1オクターブバンド：1/1オクターブの通過帯域をもつフィ...
---1/3オクターブバンド：1/3オクターブの通過帯域をもつフィ...
-まず、オクターブバンドの各フィルタのインパルス応答とゲイ...
--あらかじめ、指定したベクトルの中で、指定した値に一番近...
--以下を実行します。
#geshi(matlab){{
Fs = 8192; % サンプリング周波数（load handel の Fs ...
n = 150; % 150次のローパスデジタルフィルタ
octWidth = 1/1; % フィルタバンクの幅（オクターブ）
%周波数スケール軸と、オクターブスケール軸を定義する
freqScale = (0:10:Fs/2); % 1Hzから10Hzおきに...
octaveScale = log2(freqScale + 1); % 周波数が二倍になっ...
%バンドパスフィルタの幅を決める
bandpassFreq = [];
for count = min(octaveScale) : octWidth : max(octaveScale)
startOct = count; % このフィルタの開...
endOct = count + octWidth; % このフィルタの終...
% ベクトルoctaveScaleの中で、startOctに最も近い値の...
[startNum startData] = getNearNum(octaveScale, start...
% ベクトルoctaveScaleの中で、endOctに最も近い値の番...
[endNum endData] = getNearNum(octaveScale, endOct);
% 周波数スケールに変換
% 開始周波数が終了周波数より小さいとき（周波数が小...
if (freqScale(startNum) < freqScale(endNum))
bandpassFreq = [bandpassFreq ; freqScale(startNu...
end
end
%各フィルタバンクのインパルス応答とゲイン特性を求める
bBank = [];
hBank = [];
wBank = [];
for count = 1 : 1 : length(bandpassFreq)
% インパルス応答を求める
Wn = [ bandpassFreq(count,1)/(Fs/2) bandpassFreq(cou...
b = fir1(n,Wn,'bandpass');
bBank = [bBank ; b];
% ゲイン特性を求める
[h,w] = freqz(b,1,512);
hBank = [hBank h];
wBank = [wBank w];
end
% ゲイン特性をプロットする
filtersize = length(bandpassFreq); ...
set(gcf,'DefaultAxesColorOrder', [linspace(0,1,filtersiz...
linspace(0.5,0,filtersize)' linspace(1,0,filtersize)'])...
subplot(2,1,1); plot((wBank/pi)*(Fs/2),abs(hBank));
xlim([min(freqScale) max(freqScale)]);
xlabel('周波数[Hz]');
title( strcat(num2str(octWidth), 'オクターブフィルタバン...
subplot(2,1,2); plot( log2((wBank/pi)*(Fs/2)), abs(hBank...
xlim([min(octaveScale) max(octaveScale)]);
xlabel('オクターブ');
}}
--実行結果は以下のようになります（1/1オクターブフィルタバ...
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
---周波数軸上で対数間隔、オクターブスケール上で等間隔にフ...
---低い周波数では、ゲイン特性が予想とは異なる形になってし...
--1/3オクターブフィルタバンクを作りたいときは、変数octWid...
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
-フィルタバンクに音声データを入力します。
--以下を実行します。
#geshi(matlab){{
load handel;
% 順番に処理する
yAll = [];
for count = 1 : 1 : filtersize
% フィルタをかける
y2 = filter(bBank(count,:),1,y);
yAll = [yAll y2];
% プロット
subplot((filtersize/2),2,count);
specgram(y2,512,Fs,kaiser(500,5),475);
end
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
---フィルタの1番目（ bBank(1,:) ）は、全てNaNのため出力も...
---なお、7番目～10番目のフィルタのゲイン応答と、出力波形...
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
-分割した出力波形を足しあわせて、もとの波形を再現すること...
#geshi(matlab){{
yAllSize = size(yAll);
wav = sum(yAll(:, 2:yAllSize(2)),2); % 1番目のフィル...
wavplay(wav,Fs); % 再生
specgram(wav,512,Fs,kaiser(500,5),475); % スペクトログラム
}}
***メルフィルタバンク((バンドパスフィルタの設計手法は [[M...
-[[メルスケールとは？:http://shower.human.waseda.ac.jp/~m...
-メルバンドの各フィルタのインパルス応答とゲイン特性を求め...
--はじめに、メルスケールの変換関数 [[mellog.m:http://show...
--以下を実行します。((&color(red){このやり方だとサンプリ...
#geshi(matlab){{
Fs = 8192; % サンプリング周波数（load handel ...
n = 150; % 150次のローパスデジタルフィルタ
melFilterNum = 20; % フィルタバンクの分割数
%周波数スケール軸と、メルスケール軸を定義する
freqScale = (1:10:Fs/2); % 1Hzから10Hzおきに...
melScale = mellog(freqScale);
%フィルタバンクの幅を求める
% メルスケール上で等間隔、ただし半分は隣と重複しているの...
% 最後に足す0.5は右端の三角窓の右半分の長さ
melWidth = max(melScale) / ( melFilterNum / 2 + 0.5 );
%バンドパスフィルタの幅を決める
bandpassFreq = [];
countFilterNum = 0;
for count = min(melScale) : melWidth/2 : max(melScale)
countFilterNum = countFilterNum + 1;
if countFilterNum <= melFilterNum
startMel = count; % このフィルタ...
endMel = count + melWidth; % このフィルタ...
% ベクトルmelScaleの中で、startMelに最も近い値の...
[startNum startData] = getNearNum(melScale, star...
% ベクトルmelScaleの中で、endMelに最も近い値の番...
[endNum endData] = getNearNum(melScale, endMel);
% 周波数スケールに変換
bandpassFreq = [bandpassFreq ; freqScale(startNu...
end
end
%各フィルタバンクのインパルス応答とゲイン特性を求める
mBank = [];
bBank = [];
hBank = [];
wBank = [];
for count = 1 : 1 : length(bandpassFreq)
% 三角窓のゲイン特性を作る
startFreq = bandpassFreq(count,1);
endFreq = bandpassFreq(count,2);
% 10Hzおきに、長さ(endFreq - startFreq)Hzの三角窓を...
triWin = triang((endFreq - startFreq)/10)';
% ゲインの値を初期化
m = zeros(1,length(freqScale));
% startFreqHzからendFreqHzの区間のゲインを triWin に...
m(ceil(startFreq/10):ceil(endFreq/10-1)) = triWin;
mBank = [mBank ; m];
% 周波数軸を0～1の範囲に修正
f = freqScale - min(freqScale); % 最小を0に
f = f / max(f); % 最大を1に
% インパルス応答を求める
b = fir2(n,f,m);
bBank = [bBank ; b];
% ゲイン特性を求める
[h,w] = freqz(b,1,512);
hBank = [hBank h];
wBank = [wBank w];
end
% ゲイン特性をプロットする
filtersize = length(bandpassFreq); ...
set(gcf,'DefaultAxesColorOrder', [linspace(0,1,filtersiz...
linspace(0.5,0,filtersize)' linspace(1,0,filtersize)'])...
subplot(3,1,1); plot(freqScale, mBank);
title( strcat(num2str(melWidth), 'メルフィルタバンク') );
xlim([min(freqScale) max(freqScale)]);
xlabel('周波数[Hz](理想応答)');
subplot(3,1,2); plot((wBank/pi)*(Fs/2),abs(hBank));
xlim([min(freqScale) max(freqScale)]);
xlabel('周波数[Hz](fir2による設計)');
subplot(3,1,3); plot( mellog((wBank/pi)*(Fs/2)), abs(hBa...
xlim([min(melScale) max(melScale)]);
xlabel('メル[mel]');
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
---メル周波数スケール上で等間隔な、20個のバンドパスフィル...
---周波数スケール上では、低周波数成分ほど幅の狭い（＝周波...
---隣の三角窓と半分重複しているのは、音声認識の特徴量の慣...
-フィルタバンクに音声データを入力します。
--以下を実行します。
#geshi(matlab){{
load handel;
% 順番に処理する
yAll = [];
for count = 1 : 1 : filtersize
% フィルタをかける
y2 = filter(bBank(count,:),1,y);
yAll = [yAll y2];
% プロット
subplot((filtersize/4),4,count);
specgram(y2,512,Fs,kaiser(500,5),475);
end
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
-分割した出力波形を足しあわせて、もとの波形を再現すること...
***メルフィルタバンク解析（音声波形を分割するやり方） [#q...
-上で分割した出力波形から[[MFCC:http://shower.human.wased...
-まず、それぞれの帯域の出力波形のスペクトル解析を行い、振...
--スペクトル解析のスクリプト（サウンドデータにハニング窓...
--上に続けて、以下を実行します（このスクリプトでは、hande...
#geshi(matlab){{
% 切り出した各帯域の音声波形の中央部分を切り出して、振幅...
cuttime = 0.04; %切り出す長さ[s]
timeScale = 0 : 1/Fs : length(y)/Fs;
pre_emphasis = 0.97; %プリエンファシス係数
fftsize = 2048; %フーリエ変換の次数,...
fscale = linspace(0, Fs, fftsize); %周波数スケール（0～...
AdftSum = [];
for count = 1 : 1 : filtersize
wav = yAll(:,count);
% 中央部分を切り出し
center = fix(length(thisWav) / 2); %中心のサンプル...
wavdata = wav(center-fix(cuttime/2*Fs) : center+fix(...
time = timeScale(center-fix(cuttime/2*Fs) : center+f...
figure(1); subplot((filtersize/4),4,count);
plot(time*1000, wavdata); xlabel('時間[ms]');
% 振幅スペクトルを求め、和をとる
Adft = calcSpectrum( wavdata, fftsize, pre_emphasis );
figure(2); subplot((filtersize/4),4,count);
plot(fscale(1:fftsize/2), Adft(1:fftsize/2)); xlim([...
AdftSum = [AdftSum sum(Adft)];
end
}}
--結果は以下のようになります。
---帯域ごとに切り出した音声波形
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
---各帯域の振幅スペクトル
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
--続いて、20次元に圧縮された振幅スペクトルをプロットしま...
#geshi(matlab){{
%フィルタバンクによって melFilterNum 次元に圧縮された対...
bandpassMedianFreq = median(bandpassFreq,2); % バンド...
AdftSum_log = log10(AdftSum);
% 比較用、もとの音声波形の振幅スペクトルと対数振幅スペ...
center = fix(length(y) / 2);
wavdata = y(center-fix(cuttime/2*Fs) : center+fix(cuttim...
time = timeScale(center-fix(cuttime/2*Fs) : center+fix(c...
Adft = calcSpectrum( wavdata, fftsize, pre_emphasis ); A...
figure(3);
subplot(2,2,1); plot(fscale(1:fftsize/2), Adft(1:fftsize...
xlim([0,5000]); title(strcat('もとの音声波形の振幅スペク...
subplot(2,2,2); plot(fscale(1:fftsize/2), Adft_log(1:fft...
xlim([0,5000]); title(strcat('もとの音声波形の対数振幅ス...
subplot(2,2,3); plot(bandpassMedianFreq, AdftSum, '.', b...
xlim([0,5000]); title(strcat(int2str(melFilterNum),' 次...
subplot(2,2,4); plot(bandpassMedianFreq, AdftSum_log, '....
xlim([0,5000]); title(strcat(int2str(melFilterNum),' 次...
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
-メルスケールで等間隔な20次元に圧縮された対数振幅スペクト...
--以下を実行
#geshi(matlab){{
%コサイン変換
cpst = 12; % ケプストラム係数（低次成分何次元を取り...
AdftCpst = dct(AdftSum_log);
AdftCpst_low = AdftCpst(1:cpst);
figure(4); plot(AdftCpst, '.'); hold on; plot(AdftCpst_l...
xlabel('ケフレンシー'); ylabel('メル周波数ケプストラム係...
}}
--結果は以下のようになります（赤の部分が低次成分12次元）...
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
***メルフィルタバンク解析（振幅スペクトルを分割するやり方...
-音声波形に対してデジタルフィルタ（fir 関数）を適用するの...
--この場合、「''振幅スペクトルを求める''」→「''理想応答の...
--メルスケールの変換関数 [[mellog.m:http://shower.human.w...
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
**自己相関関数 [#oaa5df31]
-'''[[準備中:http://shower.human.waseda.ac.jp/~m-kouki/pu...
終了行:
#access
#analog
-目次
#contents
**音声を取り込む・再生する [#x740141a]
-%%WAVREAD Microsoft WAVE (".wav") サウンドファイルの読み...
#geshi(matlab){{
filename = 'a.wav'; %読み込むファイル...
[data,Fs,Bits] = wavread(filename); %dataに音声データ...
sound(data,Fs); %サンプリング周波数 Fs で再生
wavplay(data,Fs); %サンプリング周波数 Fs でWindo...
wavplay(data,Fs,'async'); %音声の再生に平行して処理を続...
disp('再生中...');
}}
-AUDIOREAD Microsoft WAVE (".wav") サウンドファイルの読み...
#geshi(matlab){{
filename = 'a.wav'; %読み込むファイル...
[data,Fs] = audioread(filename); %dataに音声データ、Fs...
sound(data,Fs); %サンプリング周波数 Fs で再生
player = audioplayer(data,Fs);
play(player); %サンプリング周波数 Fs でWindows ...
}}
-デフォルトで、以下のような音声データが用意されています。
#geshi(matlab){{
load laughter; %笑い声
sound(y, Fs);
load handel; %handelのハレルヤ
sound(y, Fs);
load mtlb; %MATLABと発音した英語音声
sound(mtlb,Fs);
}}
**波形をプロットする [#d9cc99ea]
-取り込んだ波形をプロットするには
#geshi(matlab){{
load mtlb; %MATLABと発音した音声をロード（変数mt...
plot(mtlb);
}}
-44Hzの正弦波を作って、波形をプロットするには
#geshi(matlab){{
Fs = 8192; %サンプリング周波数 8192[Hz]
f = 44; %周波数（音の高さ）44[Hz]
A = 1.0; %振幅（音の大きさ）1.0
length = 0.5; %音の再生時間 0.5[秒]
timeline = 0 : 1/Fs : length; %時間軸の作成 1/Fs 間隔で...
data = A * sin(2 * pi * f * timeline); %正弦波信号を生成
subplot(2,1,1); plot(data, '.') %ドットでプロット
subplot(2,1,2); plot(data) %通常のプロット
title('時間信号'), xlabel('time'); %グラフのラベル
}}
**スペクトログラム [#g39862ab]
-スペクトログラムは、パワースペクトル解析の結果を時間軸上...
-「MATLAB」と発音したサンプル英語音声を取り込んでプロット...
#geshi(matlab){{
load mtlb; %MATLABと発音した英語音声
disp(strcat('サンプリング周波数は', int2str(Fs), '[Hz] ...
sound(mtlb, Fs) %再生
wavwrite(mtlb, Fs, 32, 'mtlb.wav'); %量子化ビット数32で...
subplot(2,1,1);
plot(mtlb); xlabel('Time (sample)'); xlim([1 length(mtlb...
subplot(2,1,2);
spectrogram(mtlb, hamming(64), 32, 256, Fs, 'yaxis'); %...
% 64 サンプル（Wavesurferと同じ）の対称なHamming窓で分割...
% 32 サンプルのオーバーラップ
% 256 サンプルのFFT長
}}
--実行結果（音声波形とスペクトログラム）~
#ref(mtlb_mfcc_1.jpg)
---元の音声のサンプリング周波数と量子化ビット数によって、...
-マイクから2.0秒間音を取り込んで、取り込んだ音声のスペク...
#geshi(matlab){{
Fs = 44100; time = 2.0;
data = wavrecord(time*Fs, Fs, 1); %1...
spectrogram(data, hamming(64), 32, 256, Fs, 'yaxis'); %...
% 64 点（Wavesurferと同じ）のHamming窓で解析する、細かい...
% mtlb の各セグメントが 32 点のオーバーラップをする
% 離散フーリエ変換を計算するために使われる周波数点数は 256
% yaxis ： y軸上に周波数をプロットする
}}
-スペクトログラムからある瞬間（1、50、100 フレーム）にお...
#geshi(matlab){{
load mtlb; %MATLABにはじめからセットされている英語音声...
subplot(3,2,[1 2 3 4]);
spectrogram(mtlb, hamming(64), 32, 256, Fs, 'yaxis');
[s,f,t,p] = spectrogram(mtlb, hamming(64), 32, 256, Fs, '...
% f : スペクトログラムにおける周波数ベクトル
% t : 時間ベクトル
% p : 各セグメントのパワースペクトル密度
subplot(3,2,[5 6]);
plot([1:length(p)], p(:,1), [1:length(p)], p(:,50), [1:le...
legend('segment = 1', 'segment = 50', 'segment = 100');
xlabel('周波数'); ylabel('パワー');
}}
**音声波形のフレーム分割および前処理 [#gcaedeef]
-音声解析の前準備として、音声の高域強調・フレーム分割・窓...
***高域強調（プリエンファシス） [#q962d911]
-音声信号は低周波数成分が大きく、周波数が大きくなるにつれ...
-声帯の基本振動、空気中の伝搬の影響を除き、口腔内の形状に...
-プリエンファシスを行うことで、第二・第三フォルマントのよ...
-具体的には、ハイパスフィルタを使い、高周波数成分の振幅を...
--式は以下のとおり。~
#mimetex(y(n) = x(n) - px(n-1))
---x(n)はサンプル数nを変数とするデジタルサウンドデータ（...
---pはプリエンファシス係数。音声認識では0.97を使うことが...
---１サンプル前の値に対して、大きく値が変化した場合（高周...
--MATLAB では以下のようにする。
#geshi(matlab){{
load mtlb;
data = mtlb;
pre_emphasis = 0.97; %プリエンファシス係数
newData = [];
for countData = 2 : 1 : length(data)
%プリエンファシス
thisData = data(countData) - ( pre_emphasis * data(c...
newData = [newData ; thisData ];
end
subplot(2,2,1); plot(data); xlim([0,4500]); title('元の...
subplot(2,2,2); spectrogram(data, hamming(64), 32, 256, ...
subplot(2,2,3); plot(newData); xlim([0,4500]); title('プ...
subplot(2,2,4); spectrogram(newData, hamming(64), 32, 25...
wavwrite(newData, Fs, 'newmtlb.wav');
}}
---以下の関数でも同じ処理ができる。
#geshi(matlab){{
newData = filter([1 -0.97],1,data);
}}
--プリエンファシスの結果は以下のようになる。
#ref(mtlb_mfcc_5.jpg)
---基本周波数成分（スペクトログラム下部のボイスバーの部分...
***フレーム切り出し [#o5c30647]
-連続音声の分析を行う場合は、音声信号を短い時間単位で切り...
--1フレーム内に2周期程度の音声信号が含まれていれば、うま...
--通常フレーム同士は重複するようにとる。重複部分の長さを...
-例えば、お題の音声を、音声認識で一般的に用いられる フレ...
#ref(mtlb_mfcc_2.jpg)
--以下のコードで、フレームの切り出しを行う（&color(red){2...
#geshi(matlab){{
frameSize = 0.025; % フレーム長：0.025秒（...
frameShift = 0.010; % フレームシフト長：0.0...
load mtlb;
data = mtlb;
data = filter([1 -0.97],1,data); %プリエンファシス
frameSizeSample = fix( Fs * frameSize ); % フレー...
frameShiftSample = fix( Fs * frameShift ); % フレー...
disp( strcat('サウンドデータのサンプル数 ',int2str(lengt...
'フレーム長 ',int2str(frameSizeSample),'[sample], ',...
'フレームシフト長 ',int2str(frameShiftSample),'[samp...
%フレームの総数を求める
maxFrame = ...
fix( ( length(data) - (frameSizeSample - frameShiftSamp...
%フレームの個数だけ処理を繰り返す
startThisFrame = 1; % ...
endThisFrame = startThisFrame + frameSizeSample - 1; % ...
for countFrame = 1 : 1 : maxFrame
disp(strcat('フレーム番号 ',int2str(countFrame),' / ...
disp(strcat('このフレームの開始サンプル番号 ',num2st...
disp(strcat('このフレームの終了サンプル番号 ',num2st...
thisData = data(startThisFrame : endThisFrame);
% ここでこのフレームに対する計算を行う
subplot(6, 9, countFrame); plot(thisData); ylim([-3,...
axis off; %目盛りを消す
startThisFrame = startThisFrame + frameShiftSample;
endThisFrame = startThisFrame + frameSizeSample - 1;
end
}}
---結果は以下のようになる。51個のフレームに切り出された。
#ref(mtlb_mfcc_3.png,,80%)
-各フレームを結合して元の音声を復元する方法は [[MATLAB No...
***窓関数による丸め [#kf8e24d7]
-上で取り出したフレームは、単純に音声波形を切り出しただけ...
--そのため、上で切り出した音声に窓関数をかけあわせて、切...
-ここでは、[[ハミング窓:http://ja.wikipedia.org/wiki/%E7%...
--ハミング窓は以下の式で表わされる((Wikipedia より引用。...
#mimetex(w(x)=0.54-0.46cos2\pi x ,\quad if \quad 0 \leq ...
---MATLABでは以下の関数を使って作ることができる。
#geshi(matlab){{
frameSizeSample = 185; % フレーム長：サンプ...
window = hamming(frameSizeSample); % ハミング窓を生成
plot(window);
}}
--前述のスクリプトにハミング窓をかけてみる。
---[[mtlb_mfcc_1.m:http://shower.human.waseda.ac.jp/~m-ko...
---結果は以下のようになる。
#ref(mtlb_mfcc_4.png,,80%);
---各フレームの角がまるまっている。
-各フレームを結合して元の音声を復元する方法は [[MATLAB No...
-音声の両端にのみ窓かけをする方法は [[MATLAB Note/音声の...
**パワースペクトル解析 [#u704d8b6]
-パワースペクトルの詳細は [[メル周波数ケプストラム（MFCC...
-音声のパワースペクトル解析の手順は [[メル周波数ケプスト...
***パワースペクトルによるホワイトノイズの検出 [#a467125b]
-音声のパワースペクトルは低周波数成分のパワーが大きく、白...
#ref(analyzeNoize.m);
**線形予測法（LPC）とフォルマント解析 [#y6a12906]
-FFTによって出力されたスペクトル解析の結果をみると、小さ...
--スペクトル包絡の山は、音声の声道に由来する周波数成分で...
--音声では、スペクトル包絡の山の頂点がフォルマント周波数...
-スペクトル包絡をより明確に求めたい場合、線形予測法（LPC...
--&ref(a.wav); を使って以下のコードを実行
#geshi(matlab){{
[wav, Fs] = wavread('a.wav');
%中央部分を切り出し
center = fix(length(wav) / 2); %中心のサンプル番号
cuttime = 0.04; %切り出す長さ[s]
wavdata = wav(center-fix(cuttime/2*Fs) : center+fix(cutti...
%ハニング窓をかける
han_window = 0.5 - 0.5 * cos(2 * pi * [0 : 1/length(wavda...
wavdata = han_window(1:length(wavdata))' .* wavdata;
lpcsize = 32; %LPC次数
fftsize = 2048; %フーリエ変換の次数, 周波数ポイントの数
[P,f] = pyulear(wavdata, lpcsize, fftsize, Fs); %LPC分析
AP = abs(P) / fftsize; %振幅ス...
%プロット
fscale = linspace(0, Fs, fftsize); %0～Fsまでをfftsize...
%プロット（表示するのはサンプリング周波数の半分まで）
subplot(2,1,1); plot(fscale(1:fftsize/2), AP(1:fftsize/2)...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル（LPC）'); xl...
dft = fft(wavdata, fftsize); %比較用...
Adft = abs(dft) / fftsize; %振幅ス...
subplot(2,1,2); plot(fscale(1:fftsize/2), Adft(1:fftsize/...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル（FFT）'); xl...
}}
---[[LPC次数の決定方法（AICによる）:http://www.cs.takusho...
--結果
#ref(lpc_01.png);
---スペクトル包絡が取り出せました。
-スペクトル包絡のピークピッキングによって、フォルマント周...
--上記のコードに続いて、以下を実行します。
#geshi(matlab){{
peakNum = []; %ピークのフレーム番号
for peakcount = 2 : 1 : (length(AP)-1)
if ((AP(peakcount) >= AP(peakcount-1)) & (AP(peakcount)...
peakNum(peakcount) = 1; % ピークのところだけ1...
end
end
peakposition = find(peakNum == 1); % ピークの位置を決定
fmag = AP(peakposition); % 振幅を求める
frequences = f(peakposition); % フォルマント周波数...
if (length(frequences) >= 1)
F1 = frequences(1); % 第一フォルマント
end
if (length(frequences) >= 2)
F2 = frequences(2); % 第二フォルマント
end
if (length(frequences) >= 3)
F3 = frequences(3); % 第三フォルマント
end
disp(strcat('F1 : ', num2str(F1), ' F2 : ', num2str(F2), ...
}}
--結果
F1 :473.7305 F2 :1098.1934 F3 :3337.6465
---小さな値から第一フォルマント、第二フォルマント、第三フ...
---なお、ここで得られる値は確実にフォルマントを取り出せて...
--導関数（微分法）を使ったピークピッキングの方法について...
-文献
--[[応用音響学: 線形予測分析(1) LPC分析:http://ocw.u-toky...
**基本周波数解析((参考文献：[[MATLABマルチメディア信号処...
-ケプストラム法による方法
--長所 ： 正確 短所 ： ノイズに弱い
--ケプストラム法の詳細は [[メル周波数ケプストラム（MFCC）...
--&ref(a.wav); を使って以下のコードを実行
#geshi(matlab){{
[wav, Fs] = wavread('a.wav');
t = 0 : 1/Fs : length(wav)/Fs;
%中央部分を切り出し
center = fix(length(wav) / 2); %中心のサンプル番号
cuttime = 0.04; %切り出す長さ[s]
wavdata = wav(center-fix(cuttime/2*Fs) : center+fix(cutti...
time = t(center-fix(cuttime/2*Fs) : center+fix(cuttime/2*...
%ハニング窓をかける
han_window = 0.5 - 0.5 * cos(2 * pi * [0 : 1/length(wavda...
wavdata = han_window(1:length(wavdata))' .* wavdata;
fftsize = 2048; % フーリエ変換の次...
dft = fft(wavdata, fftsize); % 離散フーリエ変換
Adft = abs(dft); % 振幅スペクトル
Adft_log = log10(abs(dft)); % 対数振幅スペクトル
cps = real(ifft( Adft_log )); % ケプストラム
quefrency = linspace(0, cuttime*2, cuttime*Fs*2); %ケフレ...
subplot(2,1,1); plot(quefrency(1:length(cps))*1000, cps);...
%基本周波数ピーク抽出
cps(1:length(cps)/20) = 0; % ケプストラムの低...
cps(fftsize/2:fftsize) = 0; % ケプストラムの後...
subplot(2,1,2); plot(quefrency(1:length(cps))*1000, cps);...
xlabel('ケフレンシー[ms]');
}}
---結果
#ref(cps_001.png,,70%);
---上図の低次成分を除去した結果が下図です。基本周波数のピ...
--ピークのケフレンシーの値を求め、周波数[Hz]に変換します。
#geshi(matlab){{
peakQuefrency = quefrency( find(cps == max(cps)) );
f0 = 1 / peakQuefrency;
disp(f0);
}}
---結果
185.2416
---Praatで a.wav のピッチを解析した結果 &ref(a_pitch_praa...
---より正確に求めるには、閾値などを決めてエラー検出をする...
--よりノイズに強いケプストラム法として、'''[[改良ケプスト...
-[[自己相関:http://shower.human.waseda.ac.jp/~m-kouki/puk...
--長所 ： ノイズに強い（ロバスト） 短所 ： ケプストラム...
-'''[[STRAIGHTによる方法（限定公開）:http://shower.human....
-MATLAB以外にも、[[wavesurfer:http://www.f.waseda.jp/kiku...
**メルケプストラム解析 [#gab6732b]
-理論的な詳細は &pgid(,メル周波数ケプストラム（MFCC）); ...
-[[MFCC解析のツール/MATLAB Auditory Toolbox, Praat, HTK ...
-MATLABでメルケプストラム解析を実行するには、Malcolm Slan...
--[[Auditory Toolbox ダウンロード:http://cobweb.ecn.purdu...
--Auditory Toolbox の mfcc.m で解析します。(([[Auditory T...
--ためしに、「MATLAB」と発音したファイルを解析してみまし...
#geshi(matlab){{
%TEST_MFCC Auditory Toolbox によるメルケプストラム解析
load mtlb; %MATLABと発音した英語音声
sound(mtlb, Fs) %再生
%Auditory Toolbox によるMFCC解析
[ceps,freqresp,fb,fbrecon,freqrecon] = mfcc(mtlb, Fs, 10...
mfcc_data = ceps(2:13,:) %MFCC（C0値を除...
figure(1),subplot(2,2,1),imagesc(flipud(freqresp)),title...
figure(1),subplot(2,2,2),imagesc(flipud(fb)),title('fb')
figure(1),subplot(2,2,3),imagesc(flipud(ceps)),title('ce...
figure(1),subplot(2,2,4),imagesc(flipud(ceps(2:13,:))),t...
%ファイル書き出し
csvwrite('mtlb.txt', ceps);
}}
--結果は以下のようになります。~
&ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/mfc...
-【補足1】mfcc.m でのフレーム長・フレームシフト長について
--フレーム長は定数値を与えており、mfcc.m の40行目「window...
Fsは7418（1秒間に7418サンプル）なので、フレーム長は256/74...
--フレームシフト長は mfcc.m の引数である %%samplingRate,%...
---&color(red){frameRate は「1秒間をいくつのフレームに分...
---mfcc.m の111行目「windowStep = samplingRate/frameRate;...
---上の例では、frameRateは100を指定していますから、「1秒...
--mfcc.m の112行目「cols = fix((length(input)-windowSize)...
size(mfcc_data) = [12 50] になり、12次元のMFCCが50フレー...
-【補足2】フレーム長 0.025秒、フレームシフト長 0.010秒 で...
--WindowSize の値を Fs * 0.025 に書き換えます。
--フレームシフト長は1/100秒（1秒を100のフレームに分割）な...
**フィルタバンク解析 [#tb4d530b]
-'''[[デジタルフィルタ:http://shower.human.waseda.ac.jp/~...
--'''[[低域通過（ローパス）フィルタ:http://shower.human.w...
--'''[[高域通過（ハイパス）フィルタ:http://shower.human.w...
--'''[[帯域通過（バンドパス）フィルタ:http://shower.human...
---'''[[任意の形状のゲイン特性を指定したバンドパスフィル...
--'''[[帯域阻止（バンドストップ）フィルタ:http://shower.h...
-フィルタバンク解析とは
--フィルタバンク：バンドパスフィルタの集合体。入力信号を...
--重要でない帯域のサンプリング周波数を減らし、再び結合す...
--各帯域のスペクトルの積分を求めることで、スペクトルの次...
-以下のフィルタバンク解析において、バンドパスフィルタの設...
***2分割フィルタバンク((バンドパスフィルタの設計手法は [[...
-以下のようにします。
#geshi(matlab){{
load handel;
n = 150; % 150次のローパ...
bandpassWidth = [1 2000 ; 2001 4000]; % 1～2000Hzと、2...
%順番に処理する
subplotcount = 1;
for count = 1 : 1 : length(bandpassWidth)
% インパルス応答を求める
Wn = [ bandpassWidth(count,1)/(Fs/2) bandpassWidth(c...
b = fir1(n,Wn,'bandpass');
% ゲイン特性を求める
[h,w] = freqz(b,1,512);
subplot(length(bandpassWidth),3,subplotcount);
plot((w/pi)*(Fs/2),abs(h)); xlabel('Hz'); legend('ゲ...
subplotcount = subplotcount + 1;
% フィルタをかける
y2 = filter(b,1,y); % フィルタをかけ...
subplot(length(bandpassWidth),3,[subplotcount subplo...
specgram(y2,512,Fs,kaiser(500,5),475);
subplotcount = subplotcount + 2;
end
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
***オクターブフィルタバンク（オクターブバンド）((バンドパ...
-[[オクターブスケールとは？:http://shower.human.waseda.ac...
-オクターブフィルタバンクとは？
--[[JIS規格による、1/1、1/3オクターブバンドの上下遮断周波...
---1/1オクターブバンド：1/1オクターブの通過帯域をもつフィ...
---1/3オクターブバンド：1/3オクターブの通過帯域をもつフィ...
-まず、オクターブバンドの各フィルタのインパルス応答とゲイ...
--あらかじめ、指定したベクトルの中で、指定した値に一番近...
--以下を実行します。
#geshi(matlab){{
Fs = 8192; % サンプリング周波数（load handel の Fs ...
n = 150; % 150次のローパスデジタルフィルタ
octWidth = 1/1; % フィルタバンクの幅（オクターブ）
%周波数スケール軸と、オクターブスケール軸を定義する
freqScale = (0:10:Fs/2); % 1Hzから10Hzおきに...
octaveScale = log2(freqScale + 1); % 周波数が二倍になっ...
%バンドパスフィルタの幅を決める
bandpassFreq = [];
for count = min(octaveScale) : octWidth : max(octaveScale)
startOct = count; % このフィルタの開...
endOct = count + octWidth; % このフィルタの終...
% ベクトルoctaveScaleの中で、startOctに最も近い値の...
[startNum startData] = getNearNum(octaveScale, start...
% ベクトルoctaveScaleの中で、endOctに最も近い値の番...
[endNum endData] = getNearNum(octaveScale, endOct);
% 周波数スケールに変換
% 開始周波数が終了周波数より小さいとき（周波数が小...
if (freqScale(startNum) < freqScale(endNum))
bandpassFreq = [bandpassFreq ; freqScale(startNu...
end
end
%各フィルタバンクのインパルス応答とゲイン特性を求める
bBank = [];
hBank = [];
wBank = [];
for count = 1 : 1 : length(bandpassFreq)
% インパルス応答を求める
Wn = [ bandpassFreq(count,1)/(Fs/2) bandpassFreq(cou...
b = fir1(n,Wn,'bandpass');
bBank = [bBank ; b];
% ゲイン特性を求める
[h,w] = freqz(b,1,512);
hBank = [hBank h];
wBank = [wBank w];
end
% ゲイン特性をプロットする
filtersize = length(bandpassFreq); ...
set(gcf,'DefaultAxesColorOrder', [linspace(0,1,filtersiz...
linspace(0.5,0,filtersize)' linspace(1,0,filtersize)'])...
subplot(2,1,1); plot((wBank/pi)*(Fs/2),abs(hBank));
xlim([min(freqScale) max(freqScale)]);
xlabel('周波数[Hz]');
title( strcat(num2str(octWidth), 'オクターブフィルタバン...
subplot(2,1,2); plot( log2((wBank/pi)*(Fs/2)), abs(hBank...
xlim([min(octaveScale) max(octaveScale)]);
xlabel('オクターブ');
}}
--実行結果は以下のようになります（1/1オクターブフィルタバ...
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
---周波数軸上で対数間隔、オクターブスケール上で等間隔にフ...
---低い周波数では、ゲイン特性が予想とは異なる形になってし...
--1/3オクターブフィルタバンクを作りたいときは、変数octWid...
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
-フィルタバンクに音声データを入力します。
--以下を実行します。
#geshi(matlab){{
load handel;
% 順番に処理する
yAll = [];
for count = 1 : 1 : filtersize
% フィルタをかける
y2 = filter(bBank(count,:),1,y);
yAll = [yAll y2];
% プロット
subplot((filtersize/2),2,count);
specgram(y2,512,Fs,kaiser(500,5),475);
end
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
---フィルタの1番目（ bBank(1,:) ）は、全てNaNのため出力も...
---なお、7番目～10番目のフィルタのゲイン応答と、出力波形...
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
-分割した出力波形を足しあわせて、もとの波形を再現すること...
#geshi(matlab){{
yAllSize = size(yAll);
wav = sum(yAll(:, 2:yAllSize(2)),2); % 1番目のフィル...
wavplay(wav,Fs); % 再生
specgram(wav,512,Fs,kaiser(500,5),475); % スペクトログラム
}}
***メルフィルタバンク((バンドパスフィルタの設計手法は [[M...
-[[メルスケールとは？:http://shower.human.waseda.ac.jp/~m...
-メルバンドの各フィルタのインパルス応答とゲイン特性を求め...
--はじめに、メルスケールの変換関数 [[mellog.m:http://show...
--以下を実行します。((&color(red){このやり方だとサンプリ...
#geshi(matlab){{
Fs = 8192; % サンプリング周波数（load handel ...
n = 150; % 150次のローパスデジタルフィルタ
melFilterNum = 20; % フィルタバンクの分割数
%周波数スケール軸と、メルスケール軸を定義する
freqScale = (1:10:Fs/2); % 1Hzから10Hzおきに...
melScale = mellog(freqScale);
%フィルタバンクの幅を求める
% メルスケール上で等間隔、ただし半分は隣と重複しているの...
% 最後に足す0.5は右端の三角窓の右半分の長さ
melWidth = max(melScale) / ( melFilterNum / 2 + 0.5 );
%バンドパスフィルタの幅を決める
bandpassFreq = [];
countFilterNum = 0;
for count = min(melScale) : melWidth/2 : max(melScale)
countFilterNum = countFilterNum + 1;
if countFilterNum <= melFilterNum
startMel = count; % このフィルタ...
endMel = count + melWidth; % このフィルタ...
% ベクトルmelScaleの中で、startMelに最も近い値の...
[startNum startData] = getNearNum(melScale, star...
% ベクトルmelScaleの中で、endMelに最も近い値の番...
[endNum endData] = getNearNum(melScale, endMel);
% 周波数スケールに変換
bandpassFreq = [bandpassFreq ; freqScale(startNu...
end
end
%各フィルタバンクのインパルス応答とゲイン特性を求める
mBank = [];
bBank = [];
hBank = [];
wBank = [];
for count = 1 : 1 : length(bandpassFreq)
% 三角窓のゲイン特性を作る
startFreq = bandpassFreq(count,1);
endFreq = bandpassFreq(count,2);
% 10Hzおきに、長さ(endFreq - startFreq)Hzの三角窓を...
triWin = triang((endFreq - startFreq)/10)';
% ゲインの値を初期化
m = zeros(1,length(freqScale));
% startFreqHzからendFreqHzの区間のゲインを triWin に...
m(ceil(startFreq/10):ceil(endFreq/10-1)) = triWin;
mBank = [mBank ; m];
% 周波数軸を0～1の範囲に修正
f = freqScale - min(freqScale); % 最小を0に
f = f / max(f); % 最大を1に
% インパルス応答を求める
b = fir2(n,f,m);
bBank = [bBank ; b];
% ゲイン特性を求める
[h,w] = freqz(b,1,512);
hBank = [hBank h];
wBank = [wBank w];
end
% ゲイン特性をプロットする
filtersize = length(bandpassFreq); ...
set(gcf,'DefaultAxesColorOrder', [linspace(0,1,filtersiz...
linspace(0.5,0,filtersize)' linspace(1,0,filtersize)'])...
subplot(3,1,1); plot(freqScale, mBank);
title( strcat(num2str(melWidth), 'メルフィルタバンク') );
xlim([min(freqScale) max(freqScale)]);
xlabel('周波数[Hz](理想応答)');
subplot(3,1,2); plot((wBank/pi)*(Fs/2),abs(hBank));
xlim([min(freqScale) max(freqScale)]);
xlabel('周波数[Hz](fir2による設計)');
subplot(3,1,3); plot( mellog((wBank/pi)*(Fs/2)), abs(hBa...
xlim([min(melScale) max(melScale)]);
xlabel('メル[mel]');
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
---メル周波数スケール上で等間隔な、20個のバンドパスフィル...
---周波数スケール上では、低周波数成分ほど幅の狭い（＝周波...
---隣の三角窓と半分重複しているのは、音声認識の特徴量の慣...
-フィルタバンクに音声データを入力します。
--以下を実行します。
#geshi(matlab){{
load handel;
% 順番に処理する
yAll = [];
for count = 1 : 1 : filtersize
% フィルタをかける
y2 = filter(bBank(count,:),1,y);
yAll = [yAll y2];
% プロット
subplot((filtersize/4),4,count);
specgram(y2,512,Fs,kaiser(500,5),475);
end
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
-分割した出力波形を足しあわせて、もとの波形を再現すること...
***メルフィルタバンク解析（音声波形を分割するやり方） [#q...
-上で分割した出力波形から[[MFCC:http://shower.human.wased...
-まず、それぞれの帯域の出力波形のスペクトル解析を行い、振...
--スペクトル解析のスクリプト（サウンドデータにハニング窓...
--上に続けて、以下を実行します（このスクリプトでは、hande...
#geshi(matlab){{
% 切り出した各帯域の音声波形の中央部分を切り出して、振幅...
cuttime = 0.04; %切り出す長さ[s]
timeScale = 0 : 1/Fs : length(y)/Fs;
pre_emphasis = 0.97; %プリエンファシス係数
fftsize = 2048; %フーリエ変換の次数,...
fscale = linspace(0, Fs, fftsize); %周波数スケール（0～...
AdftSum = [];
for count = 1 : 1 : filtersize
wav = yAll(:,count);
% 中央部分を切り出し
center = fix(length(thisWav) / 2); %中心のサンプル...
wavdata = wav(center-fix(cuttime/2*Fs) : center+fix(...
time = timeScale(center-fix(cuttime/2*Fs) : center+f...
figure(1); subplot((filtersize/4),4,count);
plot(time*1000, wavdata); xlabel('時間[ms]');
% 振幅スペクトルを求め、和をとる
Adft = calcSpectrum( wavdata, fftsize, pre_emphasis );
figure(2); subplot((filtersize/4),4,count);
plot(fscale(1:fftsize/2), Adft(1:fftsize/2)); xlim([...
AdftSum = [AdftSum sum(Adft)];
end
}}
--結果は以下のようになります。
---帯域ごとに切り出した音声波形
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
---各帯域の振幅スペクトル
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
--続いて、20次元に圧縮された振幅スペクトルをプロットしま...
#geshi(matlab){{
%フィルタバンクによって melFilterNum 次元に圧縮された対...
bandpassMedianFreq = median(bandpassFreq,2); % バンド...
AdftSum_log = log10(AdftSum);
% 比較用、もとの音声波形の振幅スペクトルと対数振幅スペ...
center = fix(length(y) / 2);
wavdata = y(center-fix(cuttime/2*Fs) : center+fix(cuttim...
time = timeScale(center-fix(cuttime/2*Fs) : center+fix(c...
Adft = calcSpectrum( wavdata, fftsize, pre_emphasis ); A...
figure(3);
subplot(2,2,1); plot(fscale(1:fftsize/2), Adft(1:fftsize...
xlim([0,5000]); title(strcat('もとの音声波形の振幅スペク...
subplot(2,2,2); plot(fscale(1:fftsize/2), Adft_log(1:fft...
xlim([0,5000]); title(strcat('もとの音声波形の対数振幅ス...
subplot(2,2,3); plot(bandpassMedianFreq, AdftSum, '.', b...
xlim([0,5000]); title(strcat(int2str(melFilterNum),' 次...
subplot(2,2,4); plot(bandpassMedianFreq, AdftSum_log, '....
xlim([0,5000]); title(strcat(int2str(melFilterNum),' 次...
}}
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
-メルスケールで等間隔な20次元に圧縮された対数振幅スペクト...
--以下を実行
#geshi(matlab){{
%コサイン変換
cpst = 12; % ケプストラム係数（低次成分何次元を取り...
AdftCpst = dct(AdftSum_log);
AdftCpst_low = AdftCpst(1:cpst);
figure(4); plot(AdftCpst, '.'); hold on; plot(AdftCpst_l...
xlabel('ケフレンシー'); ylabel('メル周波数ケプストラム係...
}}
--結果は以下のようになります（赤の部分が低次成分12次元）...
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
***メルフィルタバンク解析（振幅スペクトルを分割するやり方...
-音声波形に対してデジタルフィルタ（fir 関数）を適用するの...
--この場合、「''振幅スペクトルを求める''」→「''理想応答の...
--メルスケールの変換関数 [[mellog.m:http://shower.human.w...
--結果は以下のようになります。
#ref(http://shower.human.waseda.ac.jp/~m-kouki/matlab/fil...
**自己相関関数 [#oaa5df31]
-'''[[準備中:http://shower.human.waseda.ac.jp/~m-kouki/pu...
ページ名:
既存のページ名で編集する
