サウンド処理のプログラミング

科目名：　メディアネットワーク実験IIA（2007年〜2013年）
対象：　メディアネットワークコース3年目
実験室：　M棟1階計算機室
レポート提出締切：　次週の水曜日
レポート提出先：　情報エレクトロニクス棟6階6-08
連絡先：　青木 直史（Tel: 706-6532）（E-mail: aoki@ime.ist.hokudai.ac.jp）

目的

　音はマルチメディアコンテンツを構成する重要な要素である．本実験は，MATLABによるプログラミングを通して，サウンド処理に対する理解を深めることを目的としている．

1．はじめに

　MATLABは，Mファイルと呼ばれるテキストファイルにプログラムを記述し，これを一括して実行させることができるようになっている．

　この実験では，Z:\MyDoc_xp\MyDocにMファイルを保存するためのiia_4フォルダを作成し，このフォルダを作業用のディレクトリに指定して，以降の実験を行うものとする．

2．実験1：人工的な音データ

1．ex01.m〜ex05.mを実行し，音データを生成しなさい．これらの音データを聞き比べ，音色に違いがあることを確認せよ．

% ex01.m
clear;
fs=8000;
bits=16;
length_of_s=8000;
s=zeros(1,length_of_s);
time=zeros(1,length_of_s);
for n=1:length_of_s,
    s(n)=0.5*sin(2*pi*250*(n-1)/fs);
    time(n)=n-1;
end
plot(time,s);
axis([0 128 -1 1]);
sound(s,fs);
wavwrite(s,fs,bits,'ex01.wav');

% ex02.m
clear;
fs=8000;
bits=16;
length_of_s=8000;
s=zeros(1,length_of_s);
time=zeros(1,length_of_s);
k=0;
for n=1:length_of_s,
    s(n)=0.5*(1-2*k/31);
    k=k+1;
    if k==32
        k=0;
    end
    time(n)=n-1;
end
plot(time,s);
axis([0 128 -1 1]);
sound(s,fs);
wavwrite(s,fs,bits,'ex02.wav');

% ex03.m
clear;
fs=8000;
bits=16;
length_of_s=8000;
s=zeros(1,length_of_s);
time=zeros(1,length_of_s);
k=0;
for n=1:length_of_s,
    if k<16
        s(n)=0.5;
    else
        s(n)=-0.5;
    end
    k=k+1;
    if k==32
        k=0;
    end
    time(n)=n-1;
end
plot(time,s);
axis([0 128 -1 1]);
sound(s,fs);
wavwrite(s,fs,bits,'ex03.wav');

% ex04.m
clear;
fs=8000;
bits=16;
length_of_s=8000;
s=zeros(1,length_of_s);
time=zeros(1,length_of_s);
k=8;
for n=1:length_of_s,
    if k<16
        s(n)=0.5*(-1+2*k/16);
    else
        s(n)=0.5*(1-2*(k-16)/16);
    end
    k=k+1;
    if k==32
        k=0;
    end
    time(n)=n-1;
end
plot(time,s);
axis([0 128 -1 1]);
sound(s,fs);
wavwrite(s,fs,bits,'ex04.wav');

% ex05.m
clear;
fs=8000;
bits=16;
length_of_s=8000;
s=zeros(1,length_of_s);
time=zeros(1,length_of_s);
for n=1:length_of_s,
    s(n)=rand-0.5;
    time(n)=n-1;
end
plot(time,s);
axis([0 128 -1 1]);
sound(s,fs);
wavwrite(s,fs,bits,'ex05.wav');

2．ex06.mのように，これらの音データに対してそれぞれFFT（Fast Fourier Transform）を実行し，周波数特性を表示しなさい．周波数という視点から，音色や音の高さといった音データの特徴を調べよ．なお，音の高さは波形のどのような特徴に対応するか考察せよ．

% ex06.m
clear;
fs=8000;
dft_size=1024;
s=wavread('ex01.wav');
x=zeros(1,dft_size);
for n=1:dft_size,
    x(n)=s(n);
end
X=fft(x,dft_size);
A=zeros(1,dft_size/2+1);
frequency=zeros(1,dft_size/2+1);
for k=1:dft_size/2+1,
    A(k)=abs(X(k));
    frequency(k)=(k-1)*fs/dft_size;
end
plot(frequency,A);

　図1は，周期が4msの音データの波形と，その周波数特性である．この図に示すように，この音データは250Hzの整数倍の周波数成分から構成されていることがわかる．

　こうした周波数特性を「倍音構造」と呼ぶ．倍音構造を構成する周波数成分は，小さいものから順番に，「基本音」，「2倍音」，「3倍音」と呼ばれる．基本音は音の高さに対応し，倍音の分布は音色に対応する．

図1．(a) 周期性を示す音データの波形，(b) 周波数特性

　図2は，周期性を示さない音データの波形と周波数特性である．この音データには倍音構造が見られない．そのため，音の高さを定義することができない．

図2．(a) 周期性を示さない音データの波形，(b) 周波数特性

　1980年代の家庭用ゲーム機は，こうした人工的な音データを音源として，BGMを演奏していた．ex07.mとex08.mは，その一例である．

% ex07.m
clear;
fs=44100;
bits=16;
length_of_s=35000;
s=zeros(1,length_of_s);
x0=SquareWave_(5000,45,0.1,0.25); % B5
x1=TriangleWave_(5000,150,0.2); % D4
offset=0;
for n=1:5000,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(5000,38,0.1,0.25); % D6
x1=TriangleWave_(5000,113,0.2); % G4
offset=5000;
for n=1:5000,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(5000,28,0.1,0.25); % G6
x1=TriangleWave_(5000,89,0.2); % B4
offset=10000;
for n=1:5000,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(5000,22,0.1,0.25); % B6
x1=TriangleWave_(5000,75,0.2); % D5
offset=15000;
for n=1:5000,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(15000,28,0.1,0.25); %G6
x1=TriangleWave_(15000,89,0.2); % B4
offset=20000;
for n=1:15000,
    s(offset+n)=x0(n)+x1(n);
end
sound(s,fs);
wavwrite(s,fs,bits,'ex07.wav');

% ex08.m
clear;
fs=44100;
bits=16;
length_of_s=112000;
s=zeros(1,length_of_s);
x0=SquareWave_(6200,67,0.1,0.5); % E5
x1=TriangleWave_(6200,300,0.5); % D3
offset=0;
for n=1:6200,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(6200,67,0.1,0.5); % E5
x1=TriangleWave_(6200,300,0.5); % D3
offset=7000;
for n=1:6200,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(6200,67,0.1,0.5); % E5
x1=TriangleWave_(6200,300,0.5); % D3
offset=21000;
for n=1:6200,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(6200,84,0.1,0.5); % C5
x1=TriangleWave_(6200,300,0.5); % D3
offset=35000;
for n=1:6200,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(6200,67,0.1,0.5); % E5
x1=TriangleWave_(6200,300,0.5); % D3
offset=42000;
for n=1:6200,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(6200,56,0.1,0.5); % G5
x1=TriangleWave_(6200,225,0.5); % G3
offset=56000;
for n=1:6200,
    s(offset+n)=x0(n)+x1(n);
end
x0=SquareWave_(6200,112,0.1,0.5); % G4
x1=TriangleWave_(6200,225,0.5); % G3
offset=84000;
for n=1:6200,
    s(offset+n)=x0(n)+x1(n);
end
x2=NoiseWave_(4000,0.3);
offset=0;
for n=1:4000,
    s(offset+n)=s(offset+n)+x2(n);
end
x2=NoiseWave_(1000,0.3);
offset=14000;
for n=1:1000,
    s(offset+n)=s(offset+n)+x2(n);
end
x2=NoiseWave_(4000,0.3);
offset=21000;
for n=1:4000,
    s(offset+n)=s(offset+n)+x2(n);
end
x2=NoiseWave_(1000,0.3);
offset=35000;
for n=1:1000,
    s(offset+n)=s(offset+n)+x2(n);
end
x2=NoiseWave_(4000,0.3);
offset=42000;
for n=1:4000,
    s(offset+n)=s(offset+n)+x2(n);
end
x2=NoiseWave_(4000,0.3);
offset=56000;
for n=1:4000,
    s(offset+n)=s(offset+n)+x2(n);
end
x2=NoiseWave_(4000,0.3);
offset=77000;
for n=1:4000,
    s(offset+n)=s(offset+n)+x2(n);
end
x2=NoiseWave_(1000,0.3);
offset=91000;
for n=1:1000,
    s(offset+n)=s(offset+n)+x2(n);
end
x2=NoiseWave_(1000,0.3);
offset=98000;
for n=1:1000,
    s(offset+n)=s(offset+n)+x2(n);
end
x2=NoiseWave_(1000,0.3);
offset=105000;
for n=1:1000,
    s(offset+n)=s(offset+n)+x2(n);
end
sound(s,fs);
wavwrite(s,fs,bits,'ex08.wav');

　なお，これらのプログラムは，次に示すSquareWave_関数，TriangleWave_関数，NoiseWave_関数を利用している．

% SquareWave_.m
function x=SquareWave_(duration,period,a,duty)
x=zeros(1,duration);
k=0;
for n=1:duration,
    if k<period*duty
        x(n)=a;
    else
        x(n)=-a;
    end
    k=k+1;
    if k>=period
        k=k-period;
    end
end

% TriangleWave_.m
function x=TriangleWave_(duration,period,a)
x=zeros(1,duration);
k=period/4;
for n=1:duration,
    if k<period/2
        x(n)=a*(-1+2*k/(period/2));
    else
        x(n)=a*(1-2*(k-(period/2))/(period/2));
    end
    k=k+1;
    if k>=period
        k=k-period;
    end
end

% NoiseWave_.m
function x=NoiseWave_(duration,a)
x=zeros(1,duration);
for n=1:duration,
    x(n)=a*(rand-0.5);
end

3．実験2：音データの周波数分析

　実験1では，人工的な音データの特徴について調べてみた．ここでは，実際の楽器から録音した音データの特徴について調べてみよう．

1．ex09.mのように，guitar_a4_1024.wav〜guitar_a5_1024.wav，recorder_a4_1024.wav〜recorder_a5_1024.wavに対してFFTを実行し，周波数特性を表示しなさい．周波数という視点から，音色や音の高さといった音データの特徴を調べよ．

% ex09.m
clear;
fs=8000;
dft_size=1024;
s=wavread('guitar_a4_1024.wav');
x=zeros(1,dft_size);
w=HanningWindow_(dft_size);
for n=1:dft_size,
    x(n)=s(n)*w(n);
end
X=fft(x,dft_size);
A=zeros(1,dft_size/2+1);
frequency=zeros(1,dft_size/2+1);
for k=1:dft_size/2+1,
    A(k)=abs(X(k));
    frequency(k)=(k-1)*fs/dft_size;
end
plot(frequency,A);

2．ex09.mは，FFTの前処理として，ハニング窓と呼ばれる窓関数を音データにかけている．なぜ，こうした窓処理が必要なのか考察せよ．

% HanningWindow_.m
function w=HanningWindow_(window_size)
w=zeros(1,window_size);
if mod(window_size,2)==0
    for n=1:window_size,
        w(n)=0.5-0.5*cos(2*pi*(n-1)/window_size);
    end 
else
    for n=1:window_size,
        w(n)=0.5-0.5*cos(2*pi*(n-0.5)/window_size);
    end 
end

　一般に，楽器は「12平均律音階」によって調律される．12平均律音階では，同じアルファベットの音階は，音程が1オクターブ高くなると音の高さが2倍になる．たとえば，図3に示すように，A4より音程が1オクターブ高いA5の音の高さはA4の2倍となっている．

図3．12平均律音階

4．実験3：サウンドエフェクト処理

　音楽の製作現場では，音データの音色に独特の変化をつけるサウンドエフェクト処理が利用されている．ここでは，エレキギターの代表的なサウンドエフェクト処理であるディストーションについて調べてみよう．

1．ex10.mを実行し，guitar01.wavに対してディストーションをかけ，処理前後の音を聞き比べなさい．

% ex10.m
clear;
fs=8000;
bits=16;
s0=wavread('guitar01.wav');
length_of_s=length(s0);
s1=zeros(1,length_of_s);
for n=1:length_of_s,
    s1(n)=s0(n)*100;
    if s1(n)>=0.2
        s1(n)=0.2;
    elseif s1(n)<=-0.2
        s1(n)=-0.2;
    end
end
sound(s0,fs);
pause(9);
sound(s1,fs);
wavwrite(s1,fs,bits,'ex10.wav');

2．ex11.mを実行し，ex01.wavに対してディストーションをかけ，処理前後の音を聞き比べなさい．処理前後の波形と周波数特性を調べ，ディストーションがどのような処理であるか考察せよ．

% ex11.m
clear;
fs=8000;
bits=16;
s0=wavread('ex01.wav');
length_of_s=length(s0);
s1=zeros(1,length_of_s);
for n=1:length_of_s,
    s1(n)=s0(n)*100;
    if s1(n)>=0.2
        s1(n)=0.2;
    elseif s1(n)<=-0.2
        s1(n)=-0.2;
    end
end
sound(s0,fs);
pause(2);
sound(s1,fs);
wavwrite(s1,fs,bits,'ex11.wav');

　ディストーションは非線形処理の代表例である．こうした非線形処理を行うと，本来の音データにはなかった周波数成分が現れることになる．これが，線形処理とは異なる大きな特徴である．

5．実験4：ディジタルフィルタ

　ディジタルフィルタを音データに適用すると，本来の周波数特性が変化し，結果として音色が変化する．ディジタルフィルタはFIRフィルタとIIRフィルタに大別されるが，ここではそれぞれのフィルタの処理について調べてみよう．

1．図4にFIRフィルタのブロック図を示す．ex12.mを実行し，ex02.wavに対してFIRフィルタをかけ，処理前後の音を聞き比べなさい．なお，フィルタ係数はつぎのように設定するものとする．

b=[0.0000 0.0007 0.0030 0.0045 -0.0000 -0.0131 -0.0282 -0.0295 0.0000 0.0649 0.1493 0.2215 0.2500 0.2215 0.1493 0.0649 0.0000 -0.0295 -0.0282 -0.0131 -0.0000 0.0045 0.0030 0.0007 0.0000];

処理前後の波形と周波数特性を調べ，このフィルタの特性について考察せよ．

% ex12.m
clear;
fs=8000;
bits=16;
J=24;
b=[0.0000 0.0007 0.0030 0.0045 -0.0000 -0.0131 -0.0282 -0.0295 0.0000 0.0649 0.1493 0.2215 0.2500 0.2215 0.1493 0.0649 0.0000 -0.0295 -0.0282 -0.0131 -0.0000 0.0045 0.0030 0.0007 0.0000];
x=wavread('ex02.wav');
length_of_x=length(x);
y=zeros(1,length_of_x);
for n=1:length_of_x,
    for m=1:J+1,
        if n-m+1>0
            y(n)=y(n)+b(m)*x(n-m+1);
        end
    end
end
sound(x,fs);
pause(2);
sound(y,fs);
wavwrite(y,fs,bits,'ex12.wav');

図4．FIRフィルタのブロック図

2．図5にIIRフィルタのブロック図を示す．ex13.mを実行し，ex02.wavに対してIIRフィルタをかけ，処理前後の音を聞き比べなさい．なお，フィルタ係数はつぎのように設定するものとする．

a=[1.0000 -0.9428 0.3333];
b=[0.0976 0.1953 0.0976];

処理前後の波形と周波数特性を調べ，このフィルタの特性について考察せよ．

% ex13.m
clear;
fs=8000;
bits=16;
I=2;
J=2;
a=[1.0000 -0.9428 0.3333];
b=[0.0976 0.1953 0.0976];
x=wavread('ex02.wav');
length_of_x=length(x);
y=zeros(1,length_of_x);
for n=1:length_of_x,
    for m=1:J+1,
        if n-m+1>0
            y(n)=y(n)+b(m)*x(n-m+1);
        end
    end
    for m=2:I+1,
        if n-m+1>0
            y(n)=y(n)-a(m)*y(n-m+1);
        end
    end
end
sound(x,fs);
pause(2);
sound(y,fs);
wavwrite(y,fs,bits,'ex13.wav');

図5．IIRフィルタのブロック図

6．実験5：音声に関連したサウンド処理

1．CDに収録された音楽データは，左右2チャンネルで録音されたステレオの音データとなっている．こうした音データに対して，ex14.mを実行すると，音楽データからボーカルだけを取り除くことができる．こうしたボーカルキャンセラの原理について調べよ．

% ex14.m
clear;
fs=44100;
bits=16;
s0=wavread('stereo.wav');
s0_left=s0(:,1);
s0_right=s0(:,2);
s1=s0_left-s0_right;
sound(s0,fs);
pause(35);
sound(s1,fs);
wavwrite(s1,fs,bits,'ex14.wav');

2．テレビ番組ではプライバシー保護のために，インタビューの場面の声質を加工する場合がある．ex15.mは，その一例である．こうしたボイスチェンジャの原理について調べよ．

% ex15.m
clear;
fs=44100;
bits=16;
s0=wavread('voice01.wav');
length_of_s=length(s0);
s1=zeros(1,length_of_s);
for n=1:length_of_s,
    s1(n)=s0(n)*sin(2*pi*300*(n-1)/fs);
end
sound(s0,fs);
pause(4);
sound(s1,fs);
wavwrite(s1,fs,bits,'ex15.wav');

3．ex16.mのように，本来とは異なる標本化周波数で音データを再生すると，声質が変化する．なぜ，標本化周波数を変化させると，こうした声質の変化がおきるのか考察せよ．

% ex16.m
clear;
fs=35280; % 44100 * 0.8
bits=16;
s=wavread('monoral.wav');
sound(s,fs);
wavwrite(s,fs,bits,'ex16.wav');

7．実験6：再生速度の制御

　ビデオ録画装置のなかには，音の高さを変化させずに早送り再生できるものがある．一方，音楽プレーヤーのなかには，外国語の学習や楽器の練習のために，音の高さを変化させずにスロー再生できるものがある．ここでは，こうしたサウンド処理のしくみについて調べてみよう．

1．ex17.mを実行して，音データを早送り再生しなさい．このプログラムは，音データの周期性を利用して，音データを短縮するというアルゴリズムを採用している．このアルゴリズムがどのような処理を行っているか考察せよ．

% ex17.m
clear;
fs=44100;
bits=16;

s0=wavread('voice01.wav');
length_of_s0=length(s0);

rate=1.5; % rate > 1.0

length_of_s1=ceil(length_of_s0/rate);
s1=zeros(1,length_of_s1);

template_size = floor(fs * 0.01); % 10 ms
x=zeros(1,template_size);
y=zeros(1,template_size);

offset0=0;
offset1=0;

pmin = floor(fs * 0.005); % 5ms
pmax = floor(fs * 0.02); % 20ms
r=zeros(1,pmax);

while offset0+pmax*2<=length_of_s0,
    for n=1:template_size,
        x(n)=s0(offset0+n);
    end
    max_of_r=0;
    p=pmin;
    for m=pmin:pmax,
        for n=1:template_size,
            y(n)=s0(offset0+m+n);
        end
        r(m)=0;
        for n=1:template_size,
            r(m)=r(m)+x(n)*y(n);
        end
        if r(m)>max_of_r
            max_of_r=r(m);
            p=m;
        end
    end
    for n=1:p,
        s1(offset1+n)=s0(offset0+n)*(p-n)/p;
        s1(offset1+n)=s1(offset1+n)+s0(offset0+p+n)*n/p;
    end
    q=round(p/(rate-1));
    for n=p+1:q,
        if offset0+p+n>length_of_s0
            break;
        end
        s1(offset1+n)=s0(offset0+p+n);
    end
    offset0=offset0+p+q;
    offset1=offset1+q;
end

sound(s0,fs);
pause(4);
sound(s1,fs);
wavwrite(s1,fs,bits,'ex17.wav');

　図6に早送り再生の原理を示す．このアルゴリズムは，はじめに相互相関関数を計算することで，テンプレート波形に似た部分を探し出し，波形の周期を求める．次に，オーバーラップアドによって波形をブレンドしながら，波形を短縮する．つづいて，処理を開始する時刻を更新し，以上の処理を繰り返し実行するという手順になっている．

図6．早送り再生の原理：(a) 本来の音データ，(b) 3/2倍に早送り再生する場合の音データ

2．ex18.mを実行して，音データをスロー再生しなさい．このプログラムは，音データの周期性を利用して，音データを伸長するというアルゴリズムを採用している．このアルゴリズムがどのような処理を行っているか考察せよ．

% ex18.m
clear;
fs=44100;
bits=16;

s0=wavread('voice01.wav');
length_of_s0=length(s0);

rate=0.66; % 0.5 <= rate < 1.0

length_of_s1=ceil(length_of_s0/rate);
s1=zeros(1,length_of_s1);

template_size = floor(fs * 0.01); % 10 ms
x=zeros(1,template_size);
y=zeros(1,template_size);

offset0=0;
offset1=0;

pmin = floor(fs * 0.005); % 5ms
pmax = floor(fs * 0.02); % 20ms
r=zeros(1,pmax);

while offset0+pmax*2<=length_of_s0,
    for n=1:template_size,
        x(n)=s0(offset0+n);
    end
    max_of_r=0;
    p=pmin;
    for m=pmin:pmax,
        for n=1:template_size,
            y(n)=s0(offset0+m+n);
        end
        r(m)=0;
        for n=1:template_size,
            r(m)=r(m)+x(n)*y(n);
        end
        if r(m)>max_of_r
            max_of_r=r(m);
            p=m;
        end
    end
    for n=1:p,
        s1(offset1+n)=s0(offset0+n);
    end
    for n=1:p,
        s1(offset1+p+n)=s0(offset0+p+n)*(p-n)/p;
        s1(offset1+p+n)=s1(offset1+p+n)+s0(offset0+n)*n/p;
    end
    q=round(p*rate/(1-rate));
    for n=p+1:q,
        if offset0+n>length_of_s0
            break;
        end
        s1(offset1+p+n)=s0(offset0+n);
    end
    offset0=offset0+q;
    offset1=offset1+p+q;
end

sound(s0,fs);
pause(4);
sound(s1,fs);
wavwrite(s1,fs,bits,'ex18.wav');

　図7にスロー再生の原理を示す．このアルゴリズムは，はじめに相互相関関数を計算することで，テンプレート波形に似た部分を探し出し，波形の周期を求める．次に，オーバーラップアドによって波形をブレンドしながら，波形を伸長する．つづいて，処理を開始する時刻を更新し，以上の処理を繰り返し実行するという手順になっている．

図7．スロー再生の原理：(a) 本来の音データ，(b) 2/3倍でスロー再生する場合の音データ

8．レポートについて

下記3点についてレポートを書きなさい．

1．波形と周波数特性の関係について説明しなさい．（5点）

2．FIRフィルタとIIRフィルタの特徴について説明しなさい．（5点）

3．サウンド処理に関して，自ら課題をひとつ設定し，調査しなさい．音に関連していればどのような課題でもよい．課題の例としては以下を挙げておく．具体的な考察があるものを高く評価します．（10点）

(a) 家庭用ゲーム機やシンセサイザーにおける音源について
(b) 音楽制作におけるサウンド処理について
(c) 携帯電話におけるサウンド処理について
(d) お気に入りの楽曲とその理由について

9．創成課題

1．音声または楽器音を生成するMATLABのプログラムを作成しなさい．現実の音をお手本として，できる限りリアルな音を生成することをコンテストの課題とします．

2．発表に向けて，プレゼンテーションのためのスライドを作成しなさい．どのような音を生成したか，かならずデモをしてください．優秀な作品は表彰します．

参考文献

青木 直史, ``サウンドプログラミング入門 - 音響合成の基本とC言語による実装 - ,'' 技術評論社, 2013.

青木 直史, ``C言語ではじめる音のプログラミング - サウンドエフェクトの信号処理 - ,'' オーム社, 2008.

青木 直史, ``ディジタル・サウンド処理入門 - 音のプログラミングとMATLAB（Octave・Scilab）における実際 - ,'' CQ出版社, 2006.

Last Modified: May 2 12:00 JST 2014 by Naofumi Aoki
E-mail: aoki@nis-ei.eng.hokudai.ac.jp