Miyazawa’s Pukiwiki
MATLAB Note/音声の加工
はすでに存在します。
開始行:
-目次
#contents
**音を作る [#b8b38104]
***純音を作る [#o225da40]
-正弦波(sin波)
--基本式
#mimetex(f(\theta) = b sin \theta )
---b:正弦波の振幅
---θ:角度[゜]
---円運動する点の縦座標を、角度の関数として表現した式であ...
--ここで、
#mimetex(\omega = \theta / t)
---ω:1[s]に回れる角度[゜/s]
---θ:ある時間tに回った角度[゜]
---t:時間[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~...
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){ご指摘ありがと...
#geshi(matlab){{
Fs = 44100; %サンプリング周波数4410...
f = 440; %信号の周波数(音の高さ...
A = 1.0; %信号の振幅(音の大きさ...
time = 2.5; %音の持続時間2.5秒
t = 0 : 1/Fs : time; %時間軸の作成 1/Fs間隔...
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の音を...
data2 = A * sin(2 * pi * f_2 * t); %スピーカ2の音を...
data = ([data1 ; data2])'; %2つの音を縦に結...
sound(data,Fs,16); %信号を再生
}}
--例:高さがじょじょに上がっていく純音を作るには
#geshi(matlab){{
time = 0.5; Fs = 44100; A = 1.0;
t = 0 : 1/44100 : time; %1/44100 間隔で ...
f = linspace(440, 880, length(t)); %周波数値を線形...
data = 1.0 * sin(2 * pi * f .* t); %正弦波信号を作成
sound(data, Fs, 16);
subplot(2,1,1); plot(data); title('時間信号'); xlabel('t...
subplot(2,1,2); spectrogram(data, 200, 'yaxis');
title('スペクトログラム'); xlabel('time');
}}
-参考
--[[特別研究I: MATLAB演習(2005 年度版):http://winnie.ku...
--[[平均律、自然倍音、ピタゴラス音階の比較:http://www.dai...
***母音の合成 [#ja8cf200]
-フリーの Auditory Toolbox を使います。
--[[Auditory Toolbox ダウンロード:http://cobweb.ecn.purdu...
--MakeVowel.m を使います。
**音声変換 [#u1815459]
***無音置換・雑音置換 [#hb36bb05]
-指定した秒数の区間を無音にするには、該当部分の音データの...
#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....
***音声を増幅する [#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'); %スペクトログ...
}}
***再生速度を変える [#ve6ee22c]
-再生速度を変えるには、サンプリング周波数Fsの値を大きくし...
#geshi(matlab){{
load handel; %MATLABにはじめからセットされているハレ...
Fs = Fs * 2; %サンプリング周波数の値を2倍にする
sound(y, Fs); %ロードした音...
}}
--ただしこの方法では、単位時間あたりの波の数も多くなって...
---Fsを一定にして、音声データのサンプルの中抜きや線形補間...
-ピッチを一定に保ちながら再生速度を変えるには、フレームご...
--参考:[[Ackie Sound タイムストレッチ、ピッチシフトのア...
--再生速度を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-frameShiftSamp...
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)) + t...
% フレームの追加
y2 = [y2(1:length(y2)-frameShiftSample) ; ...
y2_shift ; thisData( frameShiftSample+1 : length...
end
startFrame = startFrame + frameShiftSample;
endFrame = startFrame + frameSizeSample - 1;
end
sound(y2, Fs); % 再生
wavwrite(y2, Fs, 'handel_rateUp.wav');
}}
---結果
#ref(handel_rateUp.wav);
---やや減衰が起こっています。これは、フレームの中抜きをし...
--そこでフレーム間の位相を滑らかにするために、[[LSEE-MSTF...
#geshi(matlab){{
%LSEE-MSTFTMアルゴリズム
%【各フレームの振幅スペクトル】の計算
startFrame = 1; endFrame = startFrame + frameSizeSample -...
% 音声の長さが変わったので最大フレーム数も計算し直す
maxFrame = fix((length(y2)-(frameSizeSample-frameShiftSam...
ModifiedSTFTM = zeros(frameSizeSample*2, maxFrame); % 各...
for countFrame = 1 : 1 : maxFrame
thisData = y2(startFrame : endFrame);
window = hamming(frameSizeSample); thisData = thisDat...
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 * ...
y3 = randn(1, wgn_n)'; % 初...
% 50回のループ
for count_i = 1 : 50
y_new = zeros(length(y3), 1);
% 各フレームごとに解析
startFrame = 1; endFrame = startFrame + frameSizeSamp...
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), r...
%【振幅スペクトルとdata_LSEEMSTFTM の位相スペクト...
sound_fft_LSEEMSTFTM = ModifiedSTFTM(:,count_fram...
% 逆フーリエ変換
sound_part_LSEEMSTFTM = real(ifft(sound_fft_LSEEM...
sound_part_LSEEMSTFTM = sound_part_LSEEMSTFTM(1:f...
%音声を足しあわせ
y_new(startFrame:endFrame) = y_new(startFrame:end...
startFrame = startFrame + frameShiftSample;
endFrame = startFrame + frameSizeSample - 1;
end
y3 = y_new;
end
sound(y3, Fs); % 再生(LSEE...
wavwrite(y3, Fs, 'handel_rateUp_LSEE_MSTFTM.wav');
}}
---結果
#ref(handel_rateUp_LSEE_MSTFTM.wav);
***声の高さを変える(ピッチシフト) [#g6af6801]
-[[再生速度を変える:http://speechresearch.fiw-web.net/75....
--参考:[[Ackie Sound タイムストレッチ、ピッチシフトのア...
-以下のプログラムを実行して下さい。
#ref(pitchShift_sampleratedown.m);
#ref(Test.wav);
--実行結果
#ref(Test_pitchDOWN_sampleratedown.wav);
***声質を維持したピッチシフト [#y517f110]
-ここでは声質=声帯成分に由来するスペクトル包絡の形状と考...
--[[LPC解析:http://speechresearch.fiw-web.net/73.html#y6a...
--[[ケプストラム解析:http://speechresearch.fiw-web.net/66...
-ここではケプストラムによる方法をまとめておきます。((ピッ...
--以下を実行して下さい(リフタリングの詳細は[[メル周波数...
#ref(pitchShift.m);
#ref(Test.wav);
--結果(ピッチ上昇)
#ref(Test_pitchUP.png,,80%);
#ref(Test_pitchUP.wav);
--結果(ピッチ下降)
#ref(Test_pitchDOWN.wav);
---サンプルレート変換とフレーム中抜きによる方法&ref(Test_...
***切り出した音声の両端のノイズ(プツッという音)を除去す...
-音声の有声区間を切り出したものを再生すると、プツッという...
-例えば下記は「おばあさん」と発声した音声の/o/の途中から/...
#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://speechresearch.fiw-web...
#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; %純音1の周波数
f_2 = 330; %純音2の周波数
data1 = A * sin(2 * pi * f_1 * t); %純音1の音を生成
data2 = A * sin(2 * pi * f_2 * t); %純音2の音を生成
data = [data1 data2]; %2つの音を結合
sound(data,Fs,16); %信号を再生
}}
***フレーム同士の結合 [#ae15abec]
-[[MATLAB Note/音声の分析/フレーム切り出し:http://speechr...
-フレームシフトによる重複部分を足しあわせながら音声信号を...
#geshi(matlab){{
clear all;
load mtlb;
data_before = mtlb;
data_before = filter([1 -0.97],1,data_before); %プリ...
frameSize = 0.025; % フレーム長:0.025秒(2...
frameShift = 0.010; % フレームシフト長:0.01...
frameSizeSample = fix( Fs * frameSize ); % フレー...
frameShiftSample = fix( Fs * frameShift ); % フレー...
maxFrame = fix((length(data_before)-(frameSizeSample-fram...
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...
subplot(2, 1, 2);
spectrogram(data_after, hamming(64), 32, 256, Fs, 'yaxis'...
%再生
wavplay(data_before, Fs);
wavplay(data_after, Fs);
}}
--結果
#ref(merge1.png,,80%);
---「元の波形」としてプロットしているのはプリエンファシス...
--なお結合する前に、各フレームの窓掛けをすることが重要で...
#ref(merge2.png,,80%);
-各フレームに対してフーリエ変換→フーリエ逆変換を行った後...
#ref(mergeFrameFFT.m);
を参照して下さい。((本プログラムの内容は[[菊池英明研究室:...
---サンプルサイズが不一致でも、サンプリングレートを変換す...
***LSEE-MSTFTMアルゴリズム [#b7c530cc]
-各フレームに対してフーリエ変換を行った後で、(例えばリフ...
--加工後の振幅スペクトルに対して位相スペクトルの形状が一...
--これを解決する方法の一つが LSEE-MSTFTMアルゴリズム([[D...
---日本語の解説:[[東&川又,1999:http://ci.nii.ac.jp/naid/...
---ピッチをシフトさせた振幅スペクトルと、同じ長さのホワイ...
--試しに以下を実行して下さい(リフタリングの処理の詳細は[...
#geshi(matlab){{
clear all;
load mtlb;
data_before = mtlb;
data_before = filter([1 -0.97],1,data_before); %プリ...
frameSize = 0.025; % フレーム長:0.025秒(2...
frameShift = 0.010; % フレームシフト長:0.01...
frameSizeSample = fix( Fs * frameSize ); % フレー...
frameShiftSample = fix( Fs * frameShift ); % フレー...
maxFrame = fix((length(data_before)-(frameSizeSample-fram...
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); % 波形成...
%【各フレームの加工後の振幅スペクトル】を代入
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/matla...
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 + frameSizeSamp...
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), r...
%【加工後振幅スペクトルとdata_LSEEMSTFTM の位相ス...
sound_fft_LSEEMSTFTM = ModifiedSTFTM(:,count_fram...
% 逆フーリエ変換
sound_part_LSEEMSTFTM = real(ifft(sound_fft_LSEEM...
sound_part_LSEEMSTFTM = sound_part_LSEEMSTFTM(1:f...
%音声を足しあわせ
data_new(startFrame:endFrame) = ...
data_new(startFrame:endFrame) + sound_par...
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, 'y...
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://speechresearch.fiw-...
-50回のループを少なくする方法として、スライディングブロッ...
**周波数フィルタをかける [#c9e1d03f]
-参考:[[MATLAB音声信号処理 応用例:ハレルヤに低域通過フ...
//-例:ハレルヤに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; %フィルタリングの値(カッ...
//n = 150; %フィルタ次数
//b = fir1(n,Wn,'low'); %低域通過フィルタを設計(>...
//[h,w] = freqz(b,1,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'); %高域通過フィルタを設計(>>h...
//
//-フィルタリング関数の値を知るには、上に続いて以下を実行
//#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'...
終了行:
-目次
#contents
**音を作る [#b8b38104]
***純音を作る [#o225da40]
-正弦波(sin波)
--基本式
#mimetex(f(\theta) = b sin \theta )
---b:正弦波の振幅
---θ:角度[゜]
---円運動する点の縦座標を、角度の関数として表現した式であ...
--ここで、
#mimetex(\omega = \theta / t)
---ω:1[s]に回れる角度[゜/s]
---θ:ある時間tに回った角度[゜]
---t:時間[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~...
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){ご指摘ありがと...
#geshi(matlab){{
Fs = 44100; %サンプリング周波数4410...
f = 440; %信号の周波数(音の高さ...
A = 1.0; %信号の振幅(音の大きさ...
time = 2.5; %音の持続時間2.5秒
t = 0 : 1/Fs : time; %時間軸の作成 1/Fs間隔...
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の音を...
data2 = A * sin(2 * pi * f_2 * t); %スピーカ2の音を...
data = ([data1 ; data2])'; %2つの音を縦に結...
sound(data,Fs,16); %信号を再生
}}
--例:高さがじょじょに上がっていく純音を作るには
#geshi(matlab){{
time = 0.5; Fs = 44100; A = 1.0;
t = 0 : 1/44100 : time; %1/44100 間隔で ...
f = linspace(440, 880, length(t)); %周波数値を線形...
data = 1.0 * sin(2 * pi * f .* t); %正弦波信号を作成
sound(data, Fs, 16);
subplot(2,1,1); plot(data); title('時間信号'); xlabel('t...
subplot(2,1,2); spectrogram(data, 200, 'yaxis');
title('スペクトログラム'); xlabel('time');
}}
-参考
--[[特別研究I: MATLAB演習(2005 年度版):http://winnie.ku...
--[[平均律、自然倍音、ピタゴラス音階の比較:http://www.dai...
***母音の合成 [#ja8cf200]
-フリーの Auditory Toolbox を使います。
--[[Auditory Toolbox ダウンロード:http://cobweb.ecn.purdu...
--MakeVowel.m を使います。
**音声変換 [#u1815459]
***無音置換・雑音置換 [#hb36bb05]
-指定した秒数の区間を無音にするには、該当部分の音データの...
#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....
***音声を増幅する [#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'); %スペクトログ...
}}
***再生速度を変える [#ve6ee22c]
-再生速度を変えるには、サンプリング周波数Fsの値を大きくし...
#geshi(matlab){{
load handel; %MATLABにはじめからセットされているハレ...
Fs = Fs * 2; %サンプリング周波数の値を2倍にする
sound(y, Fs); %ロードした音...
}}
--ただしこの方法では、単位時間あたりの波の数も多くなって...
---Fsを一定にして、音声データのサンプルの中抜きや線形補間...
-ピッチを一定に保ちながら再生速度を変えるには、フレームご...
--参考:[[Ackie Sound タイムストレッチ、ピッチシフトのア...
--再生速度を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-frameShiftSamp...
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)) + t...
% フレームの追加
y2 = [y2(1:length(y2)-frameShiftSample) ; ...
y2_shift ; thisData( frameShiftSample+1 : length...
end
startFrame = startFrame + frameShiftSample;
endFrame = startFrame + frameSizeSample - 1;
end
sound(y2, Fs); % 再生
wavwrite(y2, Fs, 'handel_rateUp.wav');
}}
---結果
#ref(handel_rateUp.wav);
---やや減衰が起こっています。これは、フレームの中抜きをし...
--そこでフレーム間の位相を滑らかにするために、[[LSEE-MSTF...
#geshi(matlab){{
%LSEE-MSTFTMアルゴリズム
%【各フレームの振幅スペクトル】の計算
startFrame = 1; endFrame = startFrame + frameSizeSample -...
% 音声の長さが変わったので最大フレーム数も計算し直す
maxFrame = fix((length(y2)-(frameSizeSample-frameShiftSam...
ModifiedSTFTM = zeros(frameSizeSample*2, maxFrame); % 各...
for countFrame = 1 : 1 : maxFrame
thisData = y2(startFrame : endFrame);
window = hamming(frameSizeSample); thisData = thisDat...
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 * ...
y3 = randn(1, wgn_n)'; % 初...
% 50回のループ
for count_i = 1 : 50
y_new = zeros(length(y3), 1);
% 各フレームごとに解析
startFrame = 1; endFrame = startFrame + frameSizeSamp...
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), r...
%【振幅スペクトルとdata_LSEEMSTFTM の位相スペクト...
sound_fft_LSEEMSTFTM = ModifiedSTFTM(:,count_fram...
% 逆フーリエ変換
sound_part_LSEEMSTFTM = real(ifft(sound_fft_LSEEM...
sound_part_LSEEMSTFTM = sound_part_LSEEMSTFTM(1:f...
%音声を足しあわせ
y_new(startFrame:endFrame) = y_new(startFrame:end...
startFrame = startFrame + frameShiftSample;
endFrame = startFrame + frameSizeSample - 1;
end
y3 = y_new;
end
sound(y3, Fs); % 再生(LSEE...
wavwrite(y3, Fs, 'handel_rateUp_LSEE_MSTFTM.wav');
}}
---結果
#ref(handel_rateUp_LSEE_MSTFTM.wav);
***声の高さを変える(ピッチシフト) [#g6af6801]
-[[再生速度を変える:http://speechresearch.fiw-web.net/75....
--参考:[[Ackie Sound タイムストレッチ、ピッチシフトのア...
-以下のプログラムを実行して下さい。
#ref(pitchShift_sampleratedown.m);
#ref(Test.wav);
--実行結果
#ref(Test_pitchDOWN_sampleratedown.wav);
***声質を維持したピッチシフト [#y517f110]
-ここでは声質=声帯成分に由来するスペクトル包絡の形状と考...
--[[LPC解析:http://speechresearch.fiw-web.net/73.html#y6a...
--[[ケプストラム解析:http://speechresearch.fiw-web.net/66...
-ここではケプストラムによる方法をまとめておきます。((ピッ...
--以下を実行して下さい(リフタリングの詳細は[[メル周波数...
#ref(pitchShift.m);
#ref(Test.wav);
--結果(ピッチ上昇)
#ref(Test_pitchUP.png,,80%);
#ref(Test_pitchUP.wav);
--結果(ピッチ下降)
#ref(Test_pitchDOWN.wav);
---サンプルレート変換とフレーム中抜きによる方法&ref(Test_...
***切り出した音声の両端のノイズ(プツッという音)を除去す...
-音声の有声区間を切り出したものを再生すると、プツッという...
-例えば下記は「おばあさん」と発声した音声の/o/の途中から/...
#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://speechresearch.fiw-web...
#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; %純音1の周波数
f_2 = 330; %純音2の周波数
data1 = A * sin(2 * pi * f_1 * t); %純音1の音を生成
data2 = A * sin(2 * pi * f_2 * t); %純音2の音を生成
data = [data1 data2]; %2つの音を結合
sound(data,Fs,16); %信号を再生
}}
***フレーム同士の結合 [#ae15abec]
-[[MATLAB Note/音声の分析/フレーム切り出し:http://speechr...
-フレームシフトによる重複部分を足しあわせながら音声信号を...
#geshi(matlab){{
clear all;
load mtlb;
data_before = mtlb;
data_before = filter([1 -0.97],1,data_before); %プリ...
frameSize = 0.025; % フレーム長:0.025秒(2...
frameShift = 0.010; % フレームシフト長:0.01...
frameSizeSample = fix( Fs * frameSize ); % フレー...
frameShiftSample = fix( Fs * frameShift ); % フレー...
maxFrame = fix((length(data_before)-(frameSizeSample-fram...
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...
subplot(2, 1, 2);
spectrogram(data_after, hamming(64), 32, 256, Fs, 'yaxis'...
%再生
wavplay(data_before, Fs);
wavplay(data_after, Fs);
}}
--結果
#ref(merge1.png,,80%);
---「元の波形」としてプロットしているのはプリエンファシス...
--なお結合する前に、各フレームの窓掛けをすることが重要で...
#ref(merge2.png,,80%);
-各フレームに対してフーリエ変換→フーリエ逆変換を行った後...
#ref(mergeFrameFFT.m);
を参照して下さい。((本プログラムの内容は[[菊池英明研究室:...
---サンプルサイズが不一致でも、サンプリングレートを変換す...
***LSEE-MSTFTMアルゴリズム [#b7c530cc]
-各フレームに対してフーリエ変換を行った後で、(例えばリフ...
--加工後の振幅スペクトルに対して位相スペクトルの形状が一...
--これを解決する方法の一つが LSEE-MSTFTMアルゴリズム([[D...
---日本語の解説:[[東&川又,1999:http://ci.nii.ac.jp/naid/...
---ピッチをシフトさせた振幅スペクトルと、同じ長さのホワイ...
--試しに以下を実行して下さい(リフタリングの処理の詳細は[...
#geshi(matlab){{
clear all;
load mtlb;
data_before = mtlb;
data_before = filter([1 -0.97],1,data_before); %プリ...
frameSize = 0.025; % フレーム長:0.025秒(2...
frameShift = 0.010; % フレームシフト長:0.01...
frameSizeSample = fix( Fs * frameSize ); % フレー...
frameShiftSample = fix( Fs * frameShift ); % フレー...
maxFrame = fix((length(data_before)-(frameSizeSample-fram...
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); % 波形成...
%【各フレームの加工後の振幅スペクトル】を代入
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/matla...
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 + frameSizeSamp...
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), r...
%【加工後振幅スペクトルとdata_LSEEMSTFTM の位相ス...
sound_fft_LSEEMSTFTM = ModifiedSTFTM(:,count_fram...
% 逆フーリエ変換
sound_part_LSEEMSTFTM = real(ifft(sound_fft_LSEEM...
sound_part_LSEEMSTFTM = sound_part_LSEEMSTFTM(1:f...
%音声を足しあわせ
data_new(startFrame:endFrame) = ...
data_new(startFrame:endFrame) + sound_par...
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, 'y...
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://speechresearch.fiw-...
-50回のループを少なくする方法として、スライディングブロッ...
**周波数フィルタをかける [#c9e1d03f]
-参考:[[MATLAB音声信号処理 応用例:ハレルヤに低域通過フ...
//-例:ハレルヤに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; %フィルタリングの値(カッ...
//n = 150; %フィルタ次数
//b = fir1(n,Wn,'low'); %低域通過フィルタを設計(>...
//[h,w] = freqz(b,1,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'); %高域通過フィルタを設計(>>h...
//
//-フィルタリング関数の値を知るには、上に続いて以下を実行
//#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'...
ページ名:
既存のページ名で編集する