- 追加された行はこの色です。
- 削除された行はこの色です。
-目次
#contents
**音を作る [#b8b38104]
***純音を作る [#o225da40]
-正弦波（sin波）
--基本式
#mimetex(f(\theta) = b sin \theta )
---b：正弦波の振幅
---θ：角度[゜]
---円運動する点の縦座標を、角度の関数として表現した式である。
--ここで、
#mimetex(\omega = \theta / t)
---ω：１[s]に回れる角度[゜/s]
---θ：ある時間ｔに回った角度[゜]
---ｔ：時間[s]
--したがって、
#mimetex(f(t) = b sin \omega t )
---円運動する点の縦座標を、時間の関数として表現した式である。
--度(degree)を孤度(radian)に変換（360゜= 2π）すると、角速度と周波数の間には、以下の関係が成り立つ。
#mimetex(\omega = 2 \pi f)
---f：1[s]間に上下する波の回数（= 円運動の回数）
--したがって、正弦波の式は、以下のようになる。
#mimetex(f(t) = b sin 2 \pi f t )
---MATLABで正弦波を作る（デジタル波）
#geshi(matlab){{
b = 1.0; %正弦波の振幅
f = 12; %周波数（1[s]の間に何回同じ波形が出現するか）
Fs = 8820; %サンプリング周波数
t = 0 : 1/Fs : 1.0; %[時間軸の作成] 1/Fs 間隔で 0.0～1.0[s] まで
data = b * sin(2 * pi * f * t); %pi=円周率π
plot(data);
title('時間信号'), xlabel('time[sample]');
}}
-余弦波（cos波）
--余弦波の式は、以下のようになる。
#mimetex(f(t) = a cos 2 \pi f t )
---a：正弦波の振幅
---円運動する点の横座標を、時間の関数として表現した式である。
---MATLABで余弦波を作る（デジタル波）
#geshi(matlab){{
a = 1.0; f = 12; Fs = 8820;
t = 0 : 1/Fs : 1.0;
data = a * cos(2 * pi * f * t);
plot(data); title('時間信号'), xlabel('time[sample]');
}}
-実例
--例：周波数440Hz（「ラ」の音((&color(red){ご指摘ありがとうございました。};))）、2.5秒間のモノラル音を作るには
#geshi(matlab){{
Fs = 44100; %サンプリング周波数44100Hz
f = 440; %信号の周波数（音の高さ）440Hz
A = 1.0; %信号の振幅（音の大きさ）1.0
time = 2.5; %音の持続時間2.5秒
t = 0 : 1/Fs : time; %時間軸の作成 1/Fs間隔でtimeまで
data = A * sin(2 * pi * f * t); %正弦波信号を生成 pi=円周率π
sound(data,Fs,16) %信号を再生 16=量子化ビット数
}}
--例：ステレオ音を作るには
#geshi(matlab){{
time = 5.0; %再生時間
Fs = 44100; %サンプリング周波数
A = 1.0; %振幅（音の大きさ）
f_1 = 440; %純音A（スピーカ1）の周波数（音の高さ）
f_2 = 330; %純音B（スピーカ2）の周波数（音の高さ）
t = 0 : 1/Fs : time; %時間軸の作成
data1 = A * sin(2 * pi * f_1 * t); %スピーカ1の音を生成（1行t列）
data2 = A * sin(2 * pi * f_2 * t); %スピーカ2の音を生成（1行t列）
data = ([data1 ; data2])'; %2つの音を縦に結合してひっくり返す（t行2列）
sound(data,Fs,16); %信号を再生
}}
--例：高さがじょじょに上がっていく純音を作るには
#geshi(matlab){{
time = 0.5; Fs = 44100; A = 1.0;
t = 0 : 1/44100 : time; %1/44100 間隔で 0～time[秒]まで時間軸を作成
f = linspace(440, 880, length(t)); %周波数値を線形に変化させる（440Hz→880Hz）
data = 1.0 * sin(2 * pi * f .* t); %正弦波信号を作成
sound(data, Fs, 16);
subplot(2,1,1); plot(data); title('時間信号'); xlabel('time');
subplot(2,1,2); spectrogram(data, 200, 'yaxis');
title('スペクトログラム'); xlabel('time');
}}
-参考
--[[特別研究I: MATLAB演習（2005 年度版）:http://winnie.kuis.kyoto-u.ac.jp/~kitahara/memos/matlab-2005.pdf]]（北原鉄朗先生）
--[[平均律、自然倍音、ピタゴラス音階の比較:http://www.dai3gen.net/onkai11.htm]]
***母音の合成 [#ja8cf200]
-フリーの Auditory Toolbox を使います。
--[[Auditory Toolbox ダウンロード:http://cobweb.ecn.purdue.edu/~malcolm/interval/1998-010/]]
--MakeVowel.m を使います。
**音声変換 [#u1815459]
***無音置換・雑音置換 [#hb36bb05]
-指定した秒数の区間を無音にするには、該当部分の音データの振幅を0に置き換えます。
#geshi(matlab){{
Fs = 44100; A = 1.0; time = 1.0; f = 440;
t = 0 : 1/Fs : time; %時間軸の作成
data = A * sin(2 * pi * f * t); %正弦波信号を生成
spaceTime = 0.1; %無音部の長さ[秒]
spaceStart = 0.5; %無音部の開始位置[秒]
spaceStartFs = (spaceStart * Fs); %無音部の開始位置を計算
spaceEnd = (spaceStartFs + (spaceTime * Fs)); %無音部の終了位置を計算
data(spaceStart:spaceEnd) = 0; %指定箇所を無音に置き換え
sound(data,Fs,16); %信号を再生
plot(data, '.') %プロット
}}
-指定した秒数の区間にノイズを加えるには
--参考：[[Matlab虎の巻：音声の提示:http://www.h6.dion.ne.jp/~fff/old/technique/matlab_intro/sound.html]]
***音声を増幅する [#a415f4f2]
-音の大きさを変えるには、波形データの振幅の値を大きくしたり、小さくしたりします。
#geshi(matlab){{
load handel; %MATLABにはじめからセットされているハレルヤをロード
y = y * 2; %振幅の値を2倍にする
sound(y, Fs); %ロードした音を再生
subplot(2,1,1); plot(y); ylim([-1,1]); %波形（縦軸のサイズを指定）
subplot(2,1,2);
spectrogram(y, 570, 'yaxis'); %スペクトログラム（570セグメントに分割）
}}
***再生速度を変える [#ve6ee22c]
-再生速度を変えるには、サンプリング周波数Fsの値を大きくしたり、小さくしたりします。
#geshi(matlab){{
load handel; %MATLABにはじめからセットされているハレルヤをロード
Fs = Fs * 2; %サンプリング周波数の値を2倍にする
sound(y, Fs); %ロードした音を再生
}}
--ただしこの方法では、単位時間あたりの波の数も多くなってしまう（周波数が大きくなってしまう）ため、声の高さ（ピッチ）も高くなってしまいます。
---Fsを一定にして、音声データのサンプルの中抜きや線形補間を行っても同様の問題が起きます。
-ピッチを一定に保ちながら再生速度を変えるには、フレームごとに音声波形を増やしたり、減らしたりします。
--参考：[[Ackie Sound タイムストレッチ、ピッチシフトのアルゴリズム:http://ackiesound.ifdef.jp/tech/timestretch.html]]
--再生速度を4倍にするプログラム（4フレームにつき1回の割合でフレーム再結合）
#geshi(matlab){{
clear all;
load handel; % ハレルヤをロード
frameSize = 0.025; % フレーム長：0.025秒（25ms）
frameShift = 0.010; % フレームシフト長：0.010秒（10ms）
frameSizeSample = fix( Fs * frameSize ); % フレーム長：サンプル換算
frameShiftSample = fix( Fs * frameShift ); % フレームシフト長：サンプル換算
maxFrame = fix((length(y)-(frameSizeSample-frameShiftSample))/frameShiftSample)-1;
startFrame = 1; % フレームの開始サンプル番号
endFrame = startFrame + frameSizeSample - 1; % フレームの終了サンプル番号
% 再生速度を何倍にするか（正の整数で指定して下さい）
speechRate = 4;
% 修正後の音声データの格納用
y2 = zeros(frameSizeSample, 1);
for countFrame = 1 : 1 : maxFrame
thisData = y(startFrame : endFrame);
window = hamming(frameSizeSample); % ハミング窓を生成
thisData = thisData .* window; % 切り出したフレームにハミング窓をかける
% countFrame が speechRate で割り切れる時だけ、切り出したフレームを【足し合わせる】
if mod(countFrame, speechRate) == 0
% フレームシフト部の足しあわせ
y2_shift = ...
y2(length(y2)-frameShiftSample+1:length(y2)) + thisData(1:frameShiftSample);
% フレームの追加
y2 = [y2(1:length(y2)-frameShiftSample) ; ...
y2_shift ; thisData( frameShiftSample+1 : length(thisData) )];
end
startFrame = startFrame + frameShiftSample;
endFrame = startFrame + frameSizeSample - 1;
end
sound(y2, Fs); % 再生
wavwrite(y2, Fs, 'handel_rateUp.wav');
}}
---結果
#ref(handel_rateUp.wav);
---やや減衰が起こっています。これは、フレームの中抜きをしたため、フレームの境界がガタガタになってしまった（位相の不一致が起きた）ためと考えられます。
--そこでフレーム間の位相を滑らかにするために、[[LSEE-MSTFTMアルゴリズム:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#b7c530cc]] を使います。上記に続いて以下を実行します。
#geshi(matlab){{
%LSEE-MSTFTMアルゴリズム
%【各フレームの振幅スペクトル】の計算
startFrame = 1; endFrame = startFrame + frameSizeSample - 1;
% 音声の長さが変わったので最大フレーム数も計算し直す
maxFrame = fix((length(y2)-(frameSizeSample-frameShiftSample))/frameShiftSample)-1;
ModifiedSTFTM = zeros(frameSizeSample*2, maxFrame); % 各フレームの振幅スペクトルの格納用
for countFrame = 1 : 1 : maxFrame
thisData = y2(startFrame : endFrame);
window = hamming(frameSizeSample); thisData = thisData .* window;
fftsize = frameSizeSample * 2; % FFT次数
dft = fft(thisData, fftsize); % フーリエ変換
Adft = abs(dft); % 振幅スペクトル
ModifiedSTFTM(:, countFrame) = Adft; % 各フレームの振幅スペクトルを記録
startFrame = startFrame + frameShiftSample;
endFrame = startFrame + frameSizeSample - 1;
end
%【各フレームの位相スペクトル】の計算
wgn_sf = Fs; wgn_n = length(y2); rand('state', sum(100 * clock));
y3 = randn(1, wgn_n)'; % 初期信号（ホワイトノイズ）の生成
% 50回のループ
for count_i = 1 : 50
y_new = zeros(length(y3), 1);
% 各フレームごとに解析
startFrame = 1; endFrame = startFrame + frameSizeSample - 1;
for count_frame = 1 : maxFrame
thisData = y3(startFrame : endFrame);
window = hamming(frameSizeSample);
thisData = thisData .* window;
dft_LSEEMSTFTM = fft(thisData, fftsize);
%【data_LSEEMSTFTM の位相スペクトル】
Angdft_LSEEMSTFTM = atan2(imag(dft_LSEEMSTFTM), real(dft_LSEEMSTFTM));
%【振幅スペクトルとdata_LSEEMSTFTM の位相スペクトルから音声波形を作る】
sound_fft_LSEEMSTFTM = ModifiedSTFTM(:,count_frame) .* exp(Angdft_LSEEMSTFTM*j);
% 逆フーリエ変換
sound_part_LSEEMSTFTM = real(ifft(sound_fft_LSEEMSTFTM, fftsize));
sound_part_LSEEMSTFTM = sound_part_LSEEMSTFTM(1:fftsize/2);
%音声を足しあわせ
y_new(startFrame:endFrame) = y_new(startFrame:endFrame) + sound_part_LSEEMSTFTM;
startFrame = startFrame + frameShiftSample;
endFrame = startFrame + frameSizeSample - 1;
end
y3 = y_new;
end
sound(y3, Fs); % 再生（LSEE_MSTFTMアルゴリズム適用後）
wavwrite(y3, Fs, 'handel_rateUp_LSEE_MSTFTM.wav');
}}
---結果
#ref(handel_rateUp_LSEE_MSTFTM.wav);
***声の高さを変える（ピッチシフト） [#g6af6801]
-[[再生速度を変える:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#ve6ee22c]] との組み合わせで実現可能です。サンプリングレートを変えてピッチと再生速度を変える→再生速度だけ元の速度に戻す、の手続きを行います。
--参考：[[Ackie Sound タイムストレッチ、ピッチシフトのアルゴリズム:http://ackiesound.ifdef.jp/tech/timestretch.html]]
-以下のプログラムを実行して下さい。
--[[pitchShift_sampleratedown.m:http://shower.human.waseda.ac.jp/~m-kouki/matlab/pitchShift_sampleratedown.m.txt]]
#ref(Test.wav);
--実行結果
#ref(Test_pitchDOWN_sampleratedown.wav);
***声質を維持したピッチシフト [#y517f110]
-ここでは声質＝声帯成分に由来するスペクトル包絡の形状と考えて、スペクトル包絡を保存したピッチシフトについて検討します。
--[[LPC解析:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#y6a12906]]＋[[LSEE-MSTFTMアルゴリズム:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#b7c530cc]] による方法（LPCボコーダ） ： [[水野,小野,西本&蘇我山,2009:https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CDUQFjAA&url=http%3A%2F%2Fhil.t.u-tokyo.ac.jp%2Fpublications%2Fdownload.php%3Fbib%3DMizuno2009H10.pdf&ei=rj3cUMSEPIS0kAXbrYGwAg&usg=AFQjCNGEMU1zPh4wX2cU7CVumS2V5jcXHw&sig2=ThsZ_VsgH73HkNAXwbQ7eA&bvm=bv.1355325884,d.dGI]] など
--[[ケプストラム解析:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/66.html#qee1e455]]＋[[LSEE-MSTFTMアルゴリズム:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#b7c530cc]] による方法 ： [[澤,1996:http://library.naist.jp/dspace/handle/10061/2182]] など
-ここではケプストラムによる方法をまとめておきます。((ピッチ成分の加工方法は[[菊池英明研究室:http://www.f.waseda.jp/kikuchi/]]の江成君が検討し、実装してくれました。まことにありがとうございます。))
--以下を実行して下さい（リフタリングの詳細は[[メル周波数ケプストラム（MFCC）/音声波形の特性:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/66.html#yc248b2b]] ～ [[音声再合成:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/66.html#f411ed6a]] を参照。音声結合部の詳細はこのページの [[フレーム同士の結合:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#ae15abec]] および [[LSEE-MSTFTMアルゴリズム:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#b7c530cc]] を参照）。
---[[pitchShift.m:http://shower.human.waseda.ac.jp/~m-kouki/matlab/pitchShift.m.txt]]
#ref(Test.wav);
--結果（ピッチ上昇）
#ref(Test_pitchUP.png,,80%);
#ref(Test_pitchUP.wav);
--結果（ピッチ下降）
#ref(Test_pitchDOWN.wav);
---サンプルレート変換とフレーム中抜きによる方法&ref(Test_pitchDOWN_sampleratedown.wav); （[[声の高さを変える（ピッチシフト）:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#g6af6801]] 参照）と比べてみると、今回の方法の方が話者性が維持されています。
***切り出した音声の両端のノイズ（プツッという音）を除去する [#dcf35a5e]
-音声の有声区間を切り出したものを再生すると、プツッというノイズが入ります。
-例えば下記は「おばあさん」と発声した音声の/o/の途中から/sa/の途中までを切り出したもの。
#ref(obaHsan.wav);
-
#geshi(matlab){{
[obaHsan, fs, nbits] = wavread('obaHsan.wav');
soundsc(obaHsan,fs);
spectrogram(obaHsan, hamming(64), 32, 256, fs, 'yaxis');
}}
#ref(obaHsan1.png,,55%);
--両端が（ぶつ切りの）有声部から始まっているためノイズになります。
-そこで、[[窓関数による丸め:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#kf8e24d7]] の要領で両端に窓をかけます。以下の関数を使用します。
#ref(edgeWindow.m);
-以下のように使います。
#geshi(matlab){{
[obaHsan, fs, nbits] = wavread('obaHsan.wav');
obaHsan = edgeWindow( obaHsan, 'hamming', 600 );
soundsc(obaHsan,fs);
spectrogram(obaHsan, hamming(64), 32, 256, fs, 'yaxis');
wavwrite(obaHsan,fs,'obaHsan2.wav');
}}
#ref(obaHsan2.png,,60%);
-
#ref(obaHsan2.wav);
--角が丸まってノイズが消えました。
**音声を結合する [#l34c4f56]
***単純な結合 [#s826dba3]
-基本的には、配列の操作と同じやり方です。
#geshi(matlab){{
Fs = 44100; A = 1.0; time = 1.0;
t = 0 : 1/Fs : time; %時間軸の作成
f_1 = 440; %純音１の周波数
f_2 = 330; %純音２の周波数
data1 = A * sin(2 * pi * f_1 * t); %純音１の音を生成
data2 = A * sin(2 * pi * f_2 * t); %純音２の音を生成
data = [data1 data2]; %2つの音を結合
sound(data,Fs,16); %信号を再生
}}
***フレーム同士の結合 [#ae15abec]
-[[MATLAB Note/音声の分析/フレーム切り出し:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#o5c30647]] によって切り出した各音声フレームを再結合して元の音声に戻すには
-フレームシフトによる重複部分を足しあわせながら音声信号を復元していきます。
#geshi(matlab){{
clear all;
load mtlb;
data_before = mtlb;
data_before = filter([1 -0.97],1,data_before); %プリエンファシス
frameSize = 0.025; % フレーム長：0.025秒（25ms）
frameShift = 0.010; % フレームシフト長：0.010秒（10ms）
frameSizeSample = fix( Fs * frameSize ); % フレーム長：サンプル換算
frameShiftSample = fix( Fs * frameShift ); % フレームシフト長：サンプル換算
maxFrame = fix((length(data_before)-(frameSizeSample-frameShiftSample))/frameShiftSample)-1;
startThisFrame = 1; % フレームの開始サンプル番号
endThisFrame = startThisFrame + frameSizeSample - 1; % フレームの終了サンプル番号
%結合後の音声データの格納用
data_after = zeros(length(data_before), 1);
for countFrame = 1 : 1 : maxFrame
thisData = data_before(startThisFrame : endThisFrame);
% ここでこのフレームに対する計算を行う
%切り出したフレームの端は窓掛けで丸めておかないと、足し合わせたとき不自然な音になります。
window = hamming(frameSizeSample); % ハミング窓を生成
thisData = thisData .* window; % 切り出したフレームにハミング窓をかける
%切り出したフレームを【足し合わせる】
data_after(startThisFrame : endThisFrame) = ...
data_after(startThisFrame : endThisFrame) + thisData;
startThisFrame = startThisFrame + frameShiftSample;
endThisFrame = startThisFrame + frameSizeSample - 1;
end
%プロット
subplot(2, 1, 1);
spectrogram(data_before, hamming(64), 32, 256, Fs, 'yaxis'); title('元の波形');
subplot(2, 1, 2);
spectrogram(data_after, hamming(64), 32, 256, Fs, 'yaxis'); title('再結合した波形');
%再生
wavplay(data_before, Fs);
wavplay(data_after, Fs);
}}
--結果
#ref(merge1.png,,80%);
---「元の波形」としてプロットしているのはプリエンファシス後です。プリエンファシスを行う前の波形とは形や音圧が異なっています。
--なお結合する前に、各フレームの窓掛けをすることが重要です。窓掛けをしないとノイズ混じりの音声になります。
#ref(merge2.png,,80%);
-各フレームに対してフーリエ変換→フーリエ逆変換を行った後で、フレームを結合して元の音声波形を復元したい場合は、処理前と処理後の音声波形のサンプルサイズが一致させておきます。
--[[mergeFrameFFT.m:http://shower.human.waseda.ac.jp/~m-kouki/matlab/mergeFrameFFT.m.txt]] を参照して下さい。((本プログラムの内容は[[菊池英明研究室:http://www.f.waseda.jp/kikuchi/]]の江成君が調べてまとめてくれました。ありがとうございます。))
---サンプルサイズが不一致でも、サンプリングレートを変換するなどして対処できると思われます（未確認）
***LSEE-MSTFTMアルゴリズム [#b7c530cc]
-各フレームに対してフーリエ変換を行った後で、（例えばリフタリング処理などで）振幅スペクトルの形状を加工した場合、フレームを結合して音声波形を復元した際にノイズが入ります。
--加工後の振幅スペクトルに対して位相スペクトルの形状が一致しないことが原因です。
--これを解決する方法の一つが LSEE-MSTFTMアルゴリズム（[[DANIEL W. GRIFFIN, JAE S. LIM, "Signal estimation from modified short-time Fourier transform", ICASSP 1984:http://libra.msra.cn/Publication/1319417/signal-estimation-from-modified-short-time-fourier-transform]]） です。((法政大学 伊藤克亘先生、伊藤研 小泉さんのご指摘によります。まことにありがとうございました。))
---日本語の解説：[[東&川又,1999:http://ci.nii.ac.jp/naid/110003298786]]、[[澤,1996:http://library.naist.jp/dspace/handle/10061/2182]]、[[土井,1997:http://library.naist.jp/dspace/handle/10061/2239]]、[[水野,小野&嵯峨,2009:https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&ved=0CDUQFjAA&url=http%3A%2F%2Fhil.t.u-tokyo.ac.jp%2Fpublications%2Fdownload.php%3Fbib%3DMizuno2009ASJ09.pdf&ei=jizcUPiTOMbOlAXn2ICQAg&usg=AFQjCNGya8G_jcLXAmUmQEAyL-qPCpgT4A&sig2=rxYcVMiP4Ygcl-ozf8eK6g&bvm=bv.1355325884,d.dGI]]
---ピッチをシフトさせた振幅スペクトルと、同じ長さのホワイトノイズから作った位相スペクトルを結合して合成音声を作り、その音声の各フレームの位相スペクトルを求めて前と同じ振幅スペクトルと結合して合成音声を作り、再度その音声の各フレームの位相スペクトルを求めて前と同じ振幅スペクトルと結合して・・・を50回繰り返す。最初はフレーム間の位相スペクトルがガタガタなので合成音声にもノイズが入るが、音声全体に対するフーリエ変換を何度も繰り返すことで、位相がならされて滑らかになる。
--試しに以下を実行して下さい（リフタリングの処理の詳細は[[メル周波数ケプストラム（MFCC）/音声波形の特性:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/66.html#yc248b2b]] ～ [[音声再合成:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/66.html#f411ed6a]] を参照。音声結合部の詳細はこのページの [[フレーム同士の結合:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#ae15abec]] を参照）。
#geshi(matlab){{
clear all;
load mtlb;
data_before = mtlb;
data_before = filter([1 -0.97],1,data_before); %プリエンファシス
frameSize = 0.025; % フレーム長：0.025秒（25ms）
frameShift = 0.010; % フレームシフト長：0.010秒（10ms）
frameSizeSample = fix( Fs * frameSize ); % フレーム長：サンプル換算
frameShiftSample = fix( Fs * frameShift ); % フレームシフト長：サンプル換算
maxFrame = fix((length(data_before)-(frameSizeSample-frameShiftSample))/frameShiftSample)-1;
startThisFrame = 1; % フレームの開始サンプル番号
endThisFrame = startThisFrame + frameSizeSample - 1; % フレームの終了サンプル番号
%結合後の音声データの格納用
data_after = zeros(length(data_before), 1);
%【各フレームの加工後の振幅スペクトル】を格納
ModifiedSTFTM = zeros(frameSizeSample * 2, maxFrame);
for countFrame = 1 : 1 : maxFrame
thisData = data_before(startThisFrame : endThisFrame);
window = hamming(frameSizeSample); % ハミング窓を生成
thisData = thisData .* window;
%高速フーリエ変換
fftsize = frameSizeSample * 2; % FFT次数
dft = fft(thisData, fftsize); % フーリエ変換
%振幅スペクトル
Adft = abs(dft);
%対数振幅スペクトル
Adft_log = log10(abs(dft));
%位相スペクトル
Angdft = atan2(imag(dft), real(dft));
%ケプストラム
cps = real(ifft(Adft_log));
%リフタリング
cps_lif_P = cps;
cps_lif_P(11:length(cps)-9) = 0; %【声帯成分除去】
%plot(cps_lif_P);
%pause();
% sin波成分とcos波成分に戻す
dftSpc_P = fft(cps_lif_P, fftsize);
AdftSpc_P = abs(10 .^ dftSpc_P);
% sin波成分とcos波成分に戻す
sound_fft = AdftSpc_P .* exp(Angdft*j);
%逆フーリエ変換をして実部を取り出す
sound_part = real(ifft(sound_fft, fftsize));
sound_part_cut = sound_part(1:fftsize/2); % 波形成分は1番目～半分まで
%【各フレームの加工後の振幅スペクトル】を代入
ModifiedSTFTM(:, countFrame) = AdftSpc_P;
% 音声を結合
data_after(startThisFrame:endThisFrame) = ...
data_after(startThisFrame:endThisFrame) + sound_part_cut;
startThisFrame = startThisFrame + frameShiftSample;
endThisFrame = startThisFrame + frameSizeSample - 1;
end
%LSEE-MSTFTMアルゴリズム
% 初期信号生成（ホワイトノイズ）
% 参考：http://www.h6.dion.ne.jp/~fff/old/technique/matlab_intro/sound.html
wgn_sf = Fs;
wgn_n = length(data_after);
rand('state', sum(100 * clock));
data_LSEEMSTFTM = randn(1, wgn_n)'; % ホワイトノイズの生成
% 50回のループ
for count_i = 1 : 50
% data_LSEEMSTFTM の位相スペクトルを解析して、
% 加工後の振幅スペクトルと結合して新しい音声data_newを作る
data_new = zeros(length(data_LSEEMSTFTM), 1);
% 各フレームごとに解析
% フレームの個数だけ処理を繰り返す
startFrame = 1; endFrame = startFrame + frameSizeSample - 1;
for count_frame = 1 : maxFrame
thisData = data_LSEEMSTFTM(startFrame : endFrame);
window = hamming(frameSizeSample);
thisData = thisData .* window;
dft_LSEEMSTFTM = fft(thisData, fftsize);
%【data_LSEEMSTFTM の位相スペクトル】
Angdft_LSEEMSTFTM = atan2(imag(dft_LSEEMSTFTM), real(dft_LSEEMSTFTM));
%【加工後振幅スペクトルとdata_LSEEMSTFTM の位相スペクトルから音声波形を作る】
sound_fft_LSEEMSTFTM = ModifiedSTFTM(:,count_frame) .* exp(Angdft_LSEEMSTFTM*j);
% 逆フーリエ変換
sound_part_LSEEMSTFTM = real(ifft(sound_fft_LSEEMSTFTM, fftsize));
sound_part_LSEEMSTFTM = sound_part_LSEEMSTFTM(1:fftsize/2);
%音声を足しあわせ
data_new(startFrame:endFrame) = ...
data_new(startFrame:endFrame) + sound_part_LSEEMSTFTM;
startFrame = startFrame + frameShiftSample;
endFrame = startFrame + frameSizeSample - 1;
end
data_LSEEMSTFTM = data_new;
end
%プロット
subplot(3, 1, 1);
spectrogram(data_before, hamming(64), 32, 256, Fs, 'yaxis');
title('元の波形');
subplot(3, 1, 2);
spectrogram(data_after, hamming(64), 32, 256, Fs, 'yaxis');
title('声帯成分を除去して再結合した波形');
subplot(3, 1, 3);
spectrogram(data_LSEEMSTFTM, hamming(64), 32, 256, Fs, 'yaxis');
title('LSEE-MSTFTMアルゴリズムにかけて再結合した波形');
%再生
wavplay(data_before, Fs);
wavplay(data_after, Fs);
wavplay(data_LSEEMSTFTM, Fs);
wavwrite(data_LSEEMSTFTM, Fs, 'mtlb_F0delete.wav');
}}
--結果
#ref(LSEE_MSTFTM.png,,80%);
#ref(mtlb_F0delete.wav);
--[[声質を維持したピッチシフト:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#y517f110]] も実行してみてください。こちらの方がLSEE-MSTFTMの効果が分かりやすいです。
-50回のループを少なくする方法として、スライディングブロック法が提案されています。LSEE_MSTFTMの初期値をホワイトノイズではなく、より適切な波形を指定するというものです。詳細は [[東&川又,1999:http://ci.nii.ac.jp/naid/110003298786]] を参照してください。
**周波数フィルタをかける [#c9e1d03f]
-参考：[[MATLAB音声信号処理 応用例：ハレルヤに低域通過フィルタをかけてみる（北海道大学 豊村先生）:http://lis2.huie.hokudai.ac.jp/~toyo/MATLAB/#5-11]]
//-例：ハレルヤに800Hzの低域通過フィルタをかける
//#geshi(matlab)
//filterRate = 800;
//load handel;
//subplot(4,1,1); plot(y); ylim([-1,1]);
//subplot(4,1,2); specgram(y);
//Fn = Fs / 2; %ナイキスト周波数 = Fsの 1/2
//Wn = filterRate / Fn; %フィルタリングの値（カットオフ周波数、0＜Wn＜1）
//n = 150; %フィルタ次数
//b = fir1(n,Wn,'low'); %低域通過フィルタを設計（>>help fir1を参照）
//[h,w] = freqz(b,1,512); %周波数応答を計算（ビットレート = 512）
//y2 = filter(b,1,y); %フィルタをかける
//subplot(4,1,3); plot(y2); ylim([-1,1]); %波形をプロット
//subplot(4,1,4); spectrogram(y2, 570, 'yaxis'); %スペクトログラムをプロット
//wavwrite(y2, Fs, 'handel_filter.wav'); sound(y2, Fs); %ファイル出力と再生
//
//-高域通過フィルタを設計するには、上の8行目を以下に変更
//#geshi(matlab)
//b = fir1(n,Wn,'high'); %高域通過フィルタを設計（>>help fir1を参照）
//
//-フィルタリング関数の値を知るには、上に続いて以下を実行
//#geshi(matlab)
//figure(2);
//subplot(2,2,1); plot(b); title('インパルス応答');
//subplot(2,2,2); plot(h);
//subplot(2,2,3); plot(w);
//subplot(2,2,4); plot((w/pi)*(Fs/2),abs(h)); xlabel('Hz'); title('フィルタ特性');
