MATLABで音楽を生成してみた

明日はクリスマス！

まずはジングルベルをどうぞ！（注：音が出ます）

こちらの曲はMATLABで生成した音声ファイル（m4a）です．この記事ではその過程について説明します．

まず言っておきたいのが，私は音声・音響のプロではありません．アマチュアでもありません．せいぜい知っているのは，音は振動から成るもので正弦波のようなもので定義することによって生成することができる，ということ．

そもそもなぜこれをやることになったかというと，今月あった娘たちの幼稚園での音楽発表会がきっかけです．年長の娘のクラスは，今年の合奏で美女と野獣の「ひとりぼっちの晩餐会」を弾くことになりました．娘は木琴を担当することになったのですが，メインメロディではないためなかなか家での一人練習が物足りなくて，本人も楽しめてなさそうでした．曲は全部で６のパート（楽器）からなっていたので，妻がキーボードで１パートやったとしてもまだ物足りないという状況でした．

そこで，MATLABを使ってすべてのパートを再生してみようと考えたのです．

まずは１音から

まずは，お試しで１音を鳴らしてみます．A4は 440 Hz です．

Code

f_A4 = 440;
fs = 44100; % （Hz）サンプリング

t = 0:1/fs:1;  % 1秒の音
y = sin(2*pi*f_A4*t);

sound(y,fs)

聴いてみる

波形は想像通りのものです．

Code

plot(t,y)
xlim([0.5 0.52])  % 拡大表示

ドレミファソラシド

１音ができたらあとは周波数を変えることによって音楽が作れます！

Code

% C4 D4 E4 F4 G4 A4 B4 C5
freqs = [261.63 293.66 329.63 349.23 392 440 493.88 523.25];
t = 0:1/fs:0.5;
y = [];
for id = 1:length(freqs)
    y = [y, sin(2*pi*freqs(id)*t)];
end

sound(y,fs)

聴いてみる

これで完成！

と，言いたいところですが，ここまできたら少し拘りたくなってきました．

滑らかな音

よーく今の音階を聞いてみると，音程の変化の部分がブツブツ聞こえます．音のプロではありませんが，おそらく周波数が急に変化することによってそう聞こえるのかと思います．

ならば，それぞれの音の始まりと終わりを滑らかにしてみてはどうでしょう．つまり，音の始まりは信号を徐々に増幅し，音の終わりは信号を徐々に減衰させてみます．

これといった科学的根拠があるわけではありませんが，いろいろ試した結果，最終的に対数関数で増幅・減衰させることにしました．

Code

f_A4 = 440;
fs = 44100; % （Hz）サンプリング

t = 0:1/fs:1;  % 1秒の音
y = sin(2*pi*f_A4*t);

len = length(t);
if mod(len,2) == 0  % 偶数サンプル
    m = (1e5-[logspace(5,0,len/2),logspace(0,5,len/2)])/1e5;
else                % 奇数サンプル
    m = (1e5-[logspace(5,0,len/2+1),logspace(0,5,len/2)])/1e5;
end
y = y .* m;

sound(y,fs)

plot(t,y)

聴いてみる

なんか心地よい音ですね．

また音階で聞いてみましょう．

Code

% C4 D4 E4 F4 G4 A4 B4 C5
freqs = [261.63 293.66 329.63 349.23 392 440 493.88 523.25];
t = 0:1/fs:0.5;

len = length(t);
if mod(len,2) == 0  % 偶数サンプル
    m = (1e5-[logspace(5,0,len/2),logspace(0,5,len/2)])/1e5;
else                % 奇数サンプル
    m = (1e5-[logspace(5,0,len/2+1),logspace(0,5,len/2)])/1e5;
end

y = [];
for id = 1:length(freqs)
    y = [y, sin(2*pi*freqs(id)*t).*m];
end

sound(y,fs)

plot(y)

聴いてみる

倍音を追加

結構いいところまで来たかと思いますが，ちょっと機械音っぽいですね．楽器などの音には「倍音」の成分が含まれると聞いたことがあるので，それを追加してみましょう．（再び言いますが）私は音響のプロではないので，ここでは取り敢えず適当に倍音の成分を足し合わせています．

Code

% C4 D4 E4 F4 G4 A4 B4 C5
freqs = [261.63 293.66 329.63 349.23 392 440 493.88 523.25];
t = 0:1/fs:0.5;

len = length(t);
if mod(len,2) == 0  % 偶数サンプル
    m = (1e5-[logspace(5,0,len/2),logspace(0,5,len/2)])/1e5;
else                % 奇数サンプル
    m = (1e5-[logspace(5,0,len/2+1),logspace(0,5,len/2)])/1e5;
end

y = [];
for id = 1:length(freqs)
    x = 0.5*sin(2*pi*freqs(id)*t);
    for ii = 2:8     % 2 から 8 倍までの倍音を追加
        x = x+1/(2^ii)*sin(2*pi*ii*freqs(id)*t);
    end

    y = [y, x.*m];
end

sound(y,fs)

聴いてみる

なんか管楽器のような音ですね．

ジングルベル

さて，これを組み合わせてジングルベルを生成しましょう．楽譜はここのものを使いました．

あらかじめ楽譜から音程とタイミング情報を取り出して，各パートをテーブルに保管し構造体としてまとめたものとして保存しておきました

Code

load jingle_bell music
music

Output

music = 
    Part1: [98x2 table]
    Part2: [82x2 table]
    Part3: [63x2 table]
    Part4: [43x2 table]

パートの一部をみてみると

Code

music.Part1(1:10,:)

	Freq	EndTime
1	"D4"	0.2500
2	"B4"	0.5000
3	"A4"	0.7500
4	"G4"	1
5	"D4"	1.7500
6	"D4"	1.8750
7	"D4"	2
8	"D4"	2.2500
9	"B4"	2.5000
10	"A4"	2.7500

また，音程の周波数も事前に用意しておきます．

Code

load notes_freq notes_freq
notes_freq(1:10,:)

	Note	FrequencyHz
1 rest	"rest"	0
2 C0	"C0"	16.3500
3 C#0	"C#0"	17.3200
4 Db0	"Db0"	17.3200
5 D0	"D0"	18.3500
6 D#0	"D#0"	19.4500
7 Eb0	"Eb0"	19.4500
8 E0	"E0"	20.6000
9 F0	"F0"	21.8300
10 F#0	"F#0"	23.1200

では，各パートごとに音程とそのタイミングにあわせて周波数ベクトルを設定していきます．

Code

y = [];
for iPart = 1:4
    part = music.("Part"+iPart);

    % 0秒から始まるので，0 を追加
    timestamp = [0;part.EndTime];

    % 曲の長さの時間ベクトルを定義
    t = 0:1/fs:timestamp(end);

    % 周波数ベクトルと減衰用のベクトルを定義
    freq = zeros(size(t));
    m = ones(size(t));

    % 音程ごとに周波数と減衰用ベクトルを設定
    for id = 1:height(part)
        % 音程の周波数
        f = notes_freq.FrequencyHz(part.Freq(id));

        % 音程に対応する時間の範囲を抽出（論理配列）
        idid = t>=timestamp(id) & t<timestamp(id+1);

        % 対応する範囲の周波数を設定
        freq(idid) = f;

        % 減衰ベクトルの設定
        len = nnz(idid);
        if mod(len,2) == 0  % 偶数サンプル
            mm = (1e5-[logspace(5,0,len/2),logspace(0,5,len/2)])/1e5;
        else                % 奇数サンプル
            mm = (1e5-[logspace(5,0,len/2+1),logspace(0,5,len/2)])/1e5;
        end
        m(idid) = mm;

    end

    % 倍音成分の追加
    x = 0.5*sin(2*pi*freq.*t);
    for ii = 2:8     % 2 から 8 倍までの倍音を追加
        x = x+1/(2^ii)*sin(2*pi*ii*freq.*t);
    end

    % 各パートを別の行として定義
    y(iPart,:) = x .* m;
end

% 4つのパートを足し合わせる
Y = sum(y);

% -1から1の範囲に正規化
Y = Y / max(abs(Y));

sound(Y,fs)

もう一度ジングルベル！

Merry Christmas!