サウンドプログラミング1

科目名：　メディアネットワーク実験IA（2014年～）
対象：　メディアネットワークコース3年目
日時：　7月4日（火） - 7月5日（水）13:00～18:00
場所：　M棟1階計算機室
レポート提出締切：　7月11日（火）13:00
レポート提出先：　情報エレクトロニクス棟6階6-08
連絡先：　青木 直史（Tel: 706-6532）（E-mail: aoki@ime.ist.hokudai.ac.jp）

目的

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

1．はじめに

　MATLABは，Mファイルと呼ばれるテキストファイルにプログラムを記述し，これを一括して実行させることができるようになっている．本実験は，下記のMファイルを実行しながら進めるものとする．なお，音を確認する場合はヘッドフォンを使うこと．

2．実験1：音の特徴

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

　なお，最新版のMATLABの場合は，wavwrite関数のかわりにaudiowrite関数を使って，音データをWAVEファイルに保存すること．たとえば，音データsをex01.wavに保存する場合は，該当箇所をつぎのように書き換える必要がある．

audiowrite('ex01.wav', s, fs, 'BitsPerSample', bits);

% ex01.m
clear;
fs=8000;
bits=16;
length_of_s=8000;
s=zeros(1,length_of_s);
time=zeros(1,length_of_s);
f0=250;
gain=0.5;
for n=1:length_of_s,
    s(n)=sin(2*pi*f0*(n-1)/fs)*gain;
    time(n)=(n-1)/fs;
end
plot(time,s);
axis([0 0.016 -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);
f0=250;
gain=0.5;
t0=fs/f0;
m=0;
for n=1:length_of_s,
    s(n)=(1-2*m/t0)*gain;
    m=m+1;
    if m>=t0
        m=0;
    end
    time(n)=(n-1)/fs;
end
plot(time,s);
axis([0 0.016 -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);
f0=250;
gain=0.5;
t0=fs/f0;
m=0;
for n=1:length_of_s,
    if m<t0/2
        s(n)=1*gain;
    else
        s(n)=-1*gain;
    end
    m=m+1;
    if m>=t0
        m=0;
    end
    time(n)=(n-1)/fs;
end
plot(time,s);
axis([0 0.016 -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);
f0=250;
gain=0.5;
t0=fs/f0;
m=0;
for n=1:length_of_s,
    if m<t0/2
        s(n)=(-1+4*m/t0)*gain;
    else
        s(n)=(3-4*m/t0)*gain;
    end
    m=m+1;
    if m>=t0
        m=0;
    end
    time(n)=(n-1)/fs;
end
plot(time,s);
axis([0 0.016 -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);
gain=0.5;
for n=1:length_of_s,
    s(n)=(rand*2-1)*gain;
    time(n)=(n-1)/fs;
end
plot(time,s);
axis([0 0.016 -1 1]);
sound(s,fs);
wavwrite(s,fs,bits,'ex05.wav');

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

　なお，最新版のMATLABの場合は，wavread関数のかわりにaudioread関数を使って，WAVEファイルから音データを読み込むこと．たとえば，ex01.wavの音データを読み込み，sに格納する場合は，該当箇所をつぎのように書き換える必要がある．

[s,fs]=audioread('ex01.wav');

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

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

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

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

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

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

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

% ex07.m
clear;
fs=44100;
bits=16;
length_of_s=7000*16;
s=zeros(1,length_of_s);

% melody
x=square_wave(fs,659.26,0.1,6125); % E5
for n=1:6125,
    s(n+7000*0)=s(n+7000*0)+x(n);
end
x=square_wave(fs,659.26,0.1,6125); % E5
for n=1:6125,
    s(n+7000*1)=s(n+7000*1)+x(n);
end
x=square_wave(fs,659.26,0.1,6125); % E5
for n=1:6125,
    s(n+7000*3)=s(n+7000*3)+x(n);
end
x=square_wave(fs,523.25,0.1,6125); % C5
for n=1:6125,
    s(n+7000*5)=s(n+7000*5)+x(n);
end
x=square_wave(fs,659.26,0.1,6125); % E5
for n=1:6125,
    s(n+7000*6)=s(n+7000*6)+x(n);
end
x=square_wave(fs,783.99,0.1,6125); % G5
for n=1:6125,
    s(n+7000*8)=s(n+7000*8)+x(n);
end
x=square_wave(fs,392.00,0.1,6125); % G4
for n=1:6125,
    s(n+7000*12)=s(n+7000*12)+x(n);
end

% bass
x=triangle_wave(fs,146.83,0.2,6125); % D3
for n=1:6125,
    s(n+7000*0)=s(n+7000*0)+x(n);
end
x=triangle_wave(fs,146.83,0.2,6125); % D3
for n=1:6125,
    s(n+7000*1)=s(n+7000*1)+x(n);
end
x=triangle_wave(fs,146.83,0.2,6125); % D3
for n=1:6125,
    s(n+7000*3)=s(n+7000*3)+x(n);
end
x=triangle_wave(fs,146.83,0.2,6125); % D3
for n=1:6125,
    s(n+7000*5)=s(n+7000*5)+x(n);
end
x=triangle_wave(fs,146.83,0.2,6125); % D3
for n=1:6125,
    s(n+7000*6)=s(n+7000*6)+x(n);
end
x=triangle_wave(fs,196.00,0.2,6125); % G3
for n=1:6125,
    s(n+7000*8)=s(n+7000*8)+x(n);
end
x=triangle_wave(fs,196.00,0.2,6125); % G3
for n=1:6125,
    s(n+7000*12)=s(n+7000*12)+x(n);
end

% percussion
x=white_noise(0.1,4000);
for n=1:4000,
    s(n+7000*0)=s(n+7000*0)+x(n);
end
x=white_noise(0.1,1000);
for n=1:1000,
    s(n+7000*2)=s(n+7000*2)+x(n);
end
x=white_noise(0.1,4000);
for n=1:4000,
    s(n+7000*3)=s(n+7000*3)+x(n);
end
x=white_noise(0.1,1000);
for n=1:1000,
    s(n+7000*5)=s(n+7000*5)+x(n);
end
x=white_noise(0.1,4000);
for n=1:4000,
    s(n+7000*6)=s(n+7000*6)+x(n);
end
x=white_noise(0.1,4000);
for n=1:4000,
    s(n+7000*8)=s(n+7000*8)+x(n);
end
x=white_noise(0.1,4000);
for n=1:4000,
    s(n+7000*11)=s(n+7000*11)+x(n);
end
x=white_noise(0.1,1000);
for n=1:1000,
    s(n+7000*13)=s(n+7000*13)+x(n);
end
x=white_noise(0.1,1000);
for n=1:1000,
    s(n+7000*14)=s(n+7000*14)+x(n);
end
x=white_noise(0.1,1000);
for n=1:1000,
    s(n+7000*15)=s(n+7000*15)+x(n);
end

sound(s,fs);
wavwrite(s,fs,bits,'ex07.wav');

　なお，このプログラムは，次に示すsquare_wave関数，triangle_wave関数，white_noise関数を利用している．

% square_wave.m
function x=square_wave(fs,f0,gain,duration)
x=zeros(1,duration);
t0=fs/f0;
m=0;
for n=1:duration,
    if m<t0/2
        x(n)=1*gain;
    else
        x(n)=-1*gain;
    end
    m=m+1;
    if m>=t0
        m=0;
    end
end

% triangle_wave.m
function x=triangle_wave(fs,f0,gain,duration)
x=zeros(1,duration);
t0=fs/f0;
m=0;
for n=1:duration,
    if m<t0/2
        x(n)=(-1+4*m/t0)*gain;
    else
        x(n)=(3-4*m/t0)*gain;
    end
    m=m+1;
    if m>=t0
        m=0;
    end
end

% white_noise.m
function x=white_noise(gain,duration)
x=zeros(1,duration);
for n=1:duration,
    x(n)=(rand*2-1)*gain;
end

3．実験2：楽器音の特徴

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

1．ex08.mのように，ギターのA4からA5までの8音，リコーダーのA4からA5までの8音に対してFFTを実行し，周波数特性を表示しなさい．周波数という視点から，音色や音の高さといった音データの特徴を調べよ．

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

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

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

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

図3．12平均律音階

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

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

1．図4にFIRフィルタのブロック図を示す．ex09.mを実行し，ex05.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];

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

% ex09.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('ex05.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,'ex09.wav');

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

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

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

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

% ex10.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('ex05.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,'ex10.wav');

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

5．実験4：ギターのサウンドエフェクト処理

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

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

% ex11.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,'ex11.wav');

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

　ディストーションは非線形処理の代表例である．本来の周波数成分の配合比率が変化するだけで，新しい周波数成分が発生しないのが線形処理の特徴であるが，一方，新しい周波数成分が発生するのが非線形処理の特徴になっている．

6．実験5：音声のサウンドエフェクト処理

1．音楽CDの音データに対してex12.mを実行すると，ボーカルの歌声を取り除くことができる．こうしたボーカルキャンセラの原理について調べよ．

% ex12.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,'ex12.wav');

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

% ex13.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,'ex13.wav');

3．ex14.mのように，本来とは異なる標本化周波数で音データを再生すると声質が変化する．こうした声質の変化が生じる理由について考察せよ．

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

7．実験6：早送り再生とスロー再生

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

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

% ex15.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
    rmax=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)>rmax
            rmax=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,'ex15.wav');

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

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

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

% ex16.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
    rmax=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)>rmax
            rmax=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,'ex16.wav');

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

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

8．レポートについて

下記4点についてレポートを作成しなさい．締め切りは7月11日（火）13:00とする．

1．波形と周波数特性の関係について具体的な例をあげて説明しなさい．（5点）

2．FIRフィルタとIIRフィルタの特徴について具体的な例をあげて説明しなさい．（5点）

3．サウンドエフェクト処理の具体的な方法について調査し，そのしくみについて説明しなさい．（5点）

4．MATLABを使ってサウンドエフェクト処理のプログラムを作成しなさい．また，どのように音を加工したか，プログラムの内容について説明しなさい．処理前と処理後の音データおよびプログラムについてはメールで提出すること．（10点）

Facebookページ

https://www.facebook.com/almostanysound

参考文献

青木 直史, ``ゼロからはじめる音響学,'' 講談社, 2014.

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

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

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

Last Modified: July 4 12:00 JST 2017 by Naofumi Aoki
E-mail: aoki@ime.ist.hokudai.ac.jp