トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索   ヘルプ   最終更新のRSS

MATLAB Note/音声の加工

Last-modified: 2019-11-07 (木) 19:43:23
Top / MATLAB Note / 音声の加工

音を作る †

↑

純音を作る †

  • 正弦波（sin波）
    • 基本式
      a math image
      • b：正弦波の振幅
      • θ：角度[゜]
      • 円運動する点の縦座標を、角度の関数として表現した式である。
    • ここで、
      a math image
      • ω：１[s]に回れる角度[゜/s]
      • θ：ある時間ｔに回った角度[゜]
      • ｔ：時間[s]
    • したがって、
      a math image
      • 円運動する点の縦座標を、時間の関数として表現した式である。
    • 度(degree)を孤度(radian)に変換（360゜= 2π）すると、角速度と周波数の間には、以下の関係が成り立つ。
      a math image
      • f：1[s]間に上下する波の回数（= 円運動の回数）
    • したがって、正弦波の式は、以下のようになる。
      a math image
      • 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波）
    • 余弦波の式は、以下のようになる。
      a math image
      • a：正弦波の振幅
      • 円運動する点の横座標を、時間の関数として表現した式である。
      • 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（「ラ」の音*1）、2.5秒間のモノラル音を作るには 
       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=量子化ビット数
    • 例：ステレオ音を作るには 
       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);                      %信号を再生
    • 例：高さがじょじょに上がっていく純音を作るには 
       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');
↑

母音の合成 †

↑

音声変換 †

↑

無音置換・雑音置換 †

  • 指定した秒数の区間を無音にするには、該当部分の音データの振幅を0に置き換えます。 
    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, '.')					%プロット
↑

音声を増幅する †

  • 音の大きさを変えるには、波形データの振幅の値を大きくしたり、小さくしたりします。 
     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セグメントに分割）
↑

再生速度を変える †

  • 再生速度を変えるには、サンプリング周波数Fsの値を大きくしたり、小さくしたりします。 
     load handel;    %MATLABにはじめからセットされているハレルヤをロード
 Fs = Fs * 2;    %サンプリング周波数の値を2倍にする
 sound(y, Fs);                             %ロードした音を再生
    • ただしこの方法では、単位時間あたりの波の数も多くなってしまう（周波数が大きくなってしまう）ため、声の高さ（ピッチ）も高くなってしまいます。
      • Fsを一定にして、音声データのサンプルの中抜きや線形補間を行っても同様の問題が起きます。
  • ピッチを一定に保ちながら再生速度を変えるには、フレームごとに音声波形を増やしたり、減らしたりします。
    • 参考：Ackie Sound タイムストレッチ、ピッチシフトのアルゴリズム
    • 再生速度を4倍にするプログラム（4フレームにつき1回の割合でフレーム再結合） 
      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');
      • 結果
        filehandel_rateUp.wav
      • やや減衰が起こっています。これは、フレームの中抜きをしたため、フレームの境界がガタガタになってしまった（位相の不一致が起きた）ためと考えられます。
    • そこでフレーム間の位相を滑らかにするために、LSEE-MSTFTMアルゴリズム を使います。上記に続いて以下を実行します。 
      %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');
↑

声の高さを変える（ピッチシフト） †

↑

声質を維持したピッチシフト †

↑

切り出した音声の両端のノイズ（プツッという音）を除去する †

  • 音声の有声区間を切り出したものを再生すると、プツッというノイズが入ります。
  • 例えば下記は「おばあさん」と発声した音声の/o/の途中から/sa/の途中までを切り出したもの。
    fileobaHsan.wav
    [obaHsan, fs, nbits] = wavread('obaHsan.wav');
soundsc(obaHsan,fs);
spectrogram(obaHsan, hamming(64), 32, 256, fs, 'yaxis');
obaHsan1.png
  • 両端が（ぶつ切りの）有声部から始まっているためノイズになります。
  • 以下のように使います。 
    [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');
obaHsan2.png
  • 角が丸まってノイズが消えました。
↑

音声を結合する †

↑

単純な結合 †

  • 基本的には、配列の操作と同じやり方です。 
    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);			%信号を再生
↑

フレーム同士の結合 †

  • MATLAB Note/音声の分析/フレーム切り出し によって切り出した各音声フレームを再結合して元の音声に戻すには
  • フレームシフトによる重複部分を足しあわせながら音声信号を復元していきます。 
    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);
    • 結果
      merge1.png
      • 「元の波形」としてプロットしているのはプリエンファシス後です。プリエンファシスを行う前の波形とは形や音圧が異なっています。
  • なお結合する前に、各フレームの窓掛けをすることが重要です。窓掛けをしないとノイズ混じりの音声になります。
    merge2.png
  • 各フレームに対してフーリエ変換→フーリエ逆変換を行った後で、フレームを結合して元の音声波形を復元したい場合は、処理前と処理後の音声波形のサンプルサイズが一致させておきます。
    filemergeFrameFFT.m
    を参照して下さい。*3
    • サンプルサイズが不一致でも、サンプリングレートを変換するなどして対処できると思われます（未確認）
↑

LSEE-MSTFTMアルゴリズム †

  • 各フレームに対してフーリエ変換を行った後で、（例えばリフタリング処理などで）振幅スペクトルの形状を加工した場合、フレームを結合して音声波形を復元した際にノイズが入ります。
    • 加工後の振幅スペクトルに対して位相スペクトルの形状が一致しないことが原因です。
    • これを解決する方法の一つが LSEE-MSTFTMアルゴリズム（DANIEL W. GRIFFIN, JAE S. LIM, "Signal estimation from modified short-time Fourier transform", ICASSP 1984） です。*4
      • 日本語の解説：東&川又,1999、澤,1996、土井,1997、水野,小野&嵯峨,2009
      • ピッチをシフトさせた振幅スペクトルと、同じ長さのホワイトノイズから作った位相スペクトルを結合して合成音声を作り、その音声の各フレームの位相スペクトルを求めて前と同じ振幅スペクトルと結合して合成音声を作り、再度その音声の各フレームの位相スペクトルを求めて前と同じ振幅スペクトルと結合して・・・を50回繰り返す。最初はフレーム間の位相スペクトルがガタガタなので合成音声にもノイズが入るが、音声全体に対するフーリエ変換を何度も繰り返すことで、位相がならされて滑らかになる。
  • 試しに以下を実行して下さい（リフタリングの処理の詳細はメル周波数ケプストラム（MFCC）/音声波形の特性 ～ 音声再合成 を参照。音声結合部の詳細はこのページの フレーム同士の結合 を参照）。 
    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');
  • 結果
    LSEE_MSTFTM.png
    filemtlb_F0delete.wav
  • 50回のループを少なくする方法として、スライディングブロック法が提案されています。LSEE_MSTFTMの初期値をホワイトノイズではなく、より適切な波形を指定するというものです。詳細は 東&川又,1999 を参照してください。
↑

周波数フィルタをかける †

*1 ご指摘ありがとうございました。
*2 ピッチ成分の加工方法は菊池英明研究室の江成君が検討し、実装してくれました。まことにありがとうございます。
*3 本プログラムの内容は菊池英明研究室の江成君が調べてまとめてくれました。ありがとうございます。
*4 法政大学 伊藤克亘先生、伊藤研 小泉さんのご指摘によります。まことにありがとうございました。