MATLABを使って位相変調で音を作ったお話【KawazAdventCalendar12/10】

KawazAdventCalendar

Previous Entry: 三国殺――中国カードゲーム輸入せよ!


Kawaz Advent Calendar 10日目は@Toshがお送りします。

突然ですが、以前僕がHSPでwavを直接書き出すことで音を作ってから早1年半

EXS24でファミコン風音源を作って作曲したお話(前編)
EXS24でファミコン風音源を作って作曲したお話(後編)

音の作り方もわかってきて、音そのものに関連したプログラミングをあまりしなくなった矢先、 それを学校の実験でやることになりました

こいつで

MATLAB

MATLAB(マトラブ)は、アメリカ合衆国のMathWorks社が開発している数値解析ソフトウェアであり、その中で使うプログラミング言語の名称でもある。MATLABは、行列計算、関数とデータの可視化、アルゴリズム開発、グラフィカルインターフェイスや、他言語(C/C++/Java/Python)とのインターフェイスの機能を有している。 (from Wikipedia)

行列計算の得意なソフトウェアで、以下のような書き方をします

    % コメント
    a = 0 % 代入
    b = 0; % 「;」はecho (実行内容の表示)を防止する
    c = [1 0; 0 1]; % 2x2 単位行列
    d = [5 0; 2 3];
    e = [-1 2; -4 7];
    f = d * e; % 行列の積 output: [-5 10; -14 25]

「配列同士」の掛け算ができるように見えますが、MATLABでは 配列は「行列」とみなされており、行列の掛け算が行われます。

音も出せる

MATLABは行列の中のデータを利用して音を出すことが可能です。そのデータは-1から1の実数で示された振幅となっています。 サンプリング周波数も指定できるのが良い。

ただ ちと高い

MATLABはシェアウェアです。学生版でMATLABのみの価格が¥4,900(12/2調べ)。こういう 遊び 実験には ちと高い感じ。大学のPCには全OS1に普通に入っていることが多いんですけどね。そこで登場するのがこの言語

GNU Octave

GNU Octave は、主に数値解析を目的とした高レベルプログラミング言語である。Octaveは線形ならびに非線形問題を数値的に解くためのコマンドライン·インタフェースを提供する。また、 MATLABとほぼ互換性のある、数値実験を行うためのプログラミング言語として使用することができる。 (from Wikipedia)

なんと行列計算ができる言語は1つではありませんでした。こいつは 無料 でMacならパッケージマネージャから手に入ります。さらに重要なのが

MATLABとほぼ互換性のある

この部分で、実験の時に使ったコードがほぼそのまま使える感じです。この記事の作成にはこっちを使うこととしましょう。

データ仕様

今回作るプログラム、というよりはMATLAB / Octaveでは、出力するために音を作る際は (チャンネル数) x (サンプル数) 行列を扱います。各行列のデータは-1から1の実数で振幅を表します。今回は以下の仕様で作ってみましょう。

チャンネル: 1 (mono)
サンプリング周波数: 44100Hz (一般によく使われるやつ)

Sine Waveのデータ

だいたい想像通りです。

% createSine.m

function [y,z] = createSine(f0, length)

    % 周波数0以下のサイン波はエラーにする
    if f0 <= 0
        error('[OOPS] InvalidArgumentException thrown from createSine: Can''t create sine wave with negative frequency');
    end

    % 出力先。yはサイン波、zは対応するサンプルのタイムコード
    y=zeros(1,length);
    z=zeros(1,length);

    % サイン波...とは言ったがこの後の都合でcosで出す
    for n=1:length,
        y(n)=cos(2*pi*f0*n/44100)*1.0;
        z(n)=(n-1)/length;
    end

end

ちなみにOctaveではグラフの描画も可能です。

[sine t] = createSine(440, 44100):
figure();
plot(t, sine);
axis([-1 1 0 1]):

Kicking up a notch

しかしこれでは前と同じですしあんまりゲームと関係ありません。 今回はもうちょっとややこしい、そしてゲームの音屋っぽいことをしてみましょう。それは FM音源の再現 です。メガドラとかに入ってるアレです。 それには以下の式で表される 位相変調 を使用します。

6406df5b3a95d4c924752c40152bfadb.png

(定義上はm(t)はsine waveですが、今回はこの関数ででた波形を再度モジュレータ波形として使用することがあるのでこうしています)

本来この位相変調は携帯の通信に使われるものなのですがそれを音作りに悪用利用しようというものです。最初にこれをやったのはJohn Chowningという人だとか。

どうやって位相変調するか?

まずcos関数の中身を計算します(theta(t))mod_waveにモジュレータの波形を与え、sensitivityで波形の影響度を決めます。その結果をまたmod_waveに突っ込むと、さらなる変調が可能です。carrier_freqは音の基本の周波数となります2

% pmApply.m

function s = pmApply(length, carrier_freq, mod_wave, sensitivity)

if ~isvector(mod_wave)
   error('[OOPS] InvalidArgumentException thrown from pmApply: Parameter must be a vector');
end

theta = zeros(1, length);
ttime = zeros(1, length);
for t = 1 : length,
    % caution: 実時間との違いに注意
    rate = (t-1) / 44100;
    theta(t) = 2 * pi * carrier_freq * rate + sensitivity * mod_wave(t);
    ttime(t) = (t-1)/44100;
end

s = zeros(1, length);
for t = 1 : length,
    s(t) = cos(theta(t));
end

end

音の整形

FM音源は大抵、金属的な音を生みます。そしてその原因になるのは音に混じった高い音(高周波成分)。せっかくいい音なのに金属音のように聞こえる時は、高周波成分を取ってしまいましょう。というわけで ローパスフィルタ (LPF)も実装しちゃいます

function [b,a]=IIR_LPF(fc,Q)
fc=tan(pi*fc)/(2*pi);
a(1)=1+2*pi*fc/Q+4*pi*pi*fc*fc;
a(2)=(8*pi*pi*fc*fc-2)/a(1);
a(3)=(1-2*pi*fc/Q+4*pi*pi*fc*fc)/a(1);
b(1)=4*pi*pi*fc*fc/a(1);
b(2)=8*pi*pi*fc*fc/a(1);
b(3)=4*pi*pi*fc*fc/a(1);
a(1)=1;

ここで実装したのは離散時間 IIR (無限インパルス応答)ローパスフィルタです。簡単に言うと、何サンプルか前の入力と出力を元に現在の出力が決まるタイプのフィルタです。ここでは2遅延まで入力しています。

ADSR

そして音の音量は一定ではありません。電子楽器でこれを表現するのに、ADSRという概念が用いられます。 これは信号をA(立ち上がり) D(減衰) S(持続音量) R(余韻) の4つで制御する機能です。これも実装しちゃいましょう。

function e=ADSR(A,D,S,R,gate,duration)
e=zeros(1,duration);
if A~=0
    for n=1:A,
        e(n)=1-exp(-5*(n-1)/A);
    end
end
if D~=0
    for n=A+1:gate,
        e(n)=S+(1-S)*exp(-5*(n-1-A)/D);
    end
else
    for n=A+1:gate,
        e(n)=S;
    end  
end
if R~=0
    for n=gate+1:duration,
        e(n)=e(gate)*exp(-5*(n-gate+2)/R);
    end
end

Sは0.0から1.0、それ以外はサンプル数で指定します。durationは全体の長さ、gateはRが始まるまでの長さを指定します。
この関数については、返された行列を音の各サンプルに掛ける形で処理します。

音を作ってみた

ここまでくれば音を作りに行けます。ここまでで作ったコードで、音を組み立てていきます。

Electric Bass

2つの波形を合成してみました。片方に変調を適用しています。

function y=Bass(freq, length)
    length = round(length);

    [s1, ~] = createSine(freq / 8, length);

    s1 = pmApply(length, freq, s1.*ADSR(0, length, 0.5, 0, length, length), 12);

    s0 = s1.*ADSR(441,30870,0.2,round(length*0.05),round(length*0.98),round(length));

    [s2, ~] = createSine(freq / 8, length);
    s2 = s2.*ADSR(882,30870,0.2,round(length*0.05),round(length*0.98),round(length));

    y = s0*0.6+ s2*0.4;

    %figure(1);
    %plot(s0);
    %axis([0 size(s0,2) -1 1]);
    %s0 = s0.';
    %sound(s0, 44100, 16);

end

Saxphone

こちらにはLPFを適用しています。その分少しばかりコードが長いです。
ベースの時と同様、複数の波形を最終的に合成しています。
波形s1への(zerosの次の)代入の時に、mod_waveにものすごく周波数の低い(4Hz)Sine Waveが入っています。 これはビブラートを掛けるという意味があります。実際に波形を出力してみると、波形全体が上下しているのが非常に興味深いです。

function y = Saxphoneish(freq, length_of_s)

length_of_s=round(length_of_s);
s0=zeros(1,length_of_s);
s3=zeros(1,length_of_s);

s1 = pmApply(length_of_s, freq, createSine(4,length_of_s), 0.5); % vibrato
s1 = pmApply(length_of_s, freq, s1, 1);
s1 = pmApply(length_of_s, freq*4, s1, 1);
s1 = pmApply(length_of_s, freq, s1, 1);
env = ADSR(44100*0.08, length_of_s, 0.7, 44100*0.05, round(length_of_s*0.8), length_of_s);
s1=s1.*env;
fc= 2000/length_of_s;
Q=1/sqrt(2);
[b,a]=IIR_LPF(fc,Q);
I=2;
J=2;
for n=1:length_of_s,
    for m=1:J+1,
        if n-m+1>0
            s0(n)=s0(n)+b(m)*s1(n-m+1);
        end
    end
    for m=2:I+1,
        if n-m+1>0
            s0(n)=s0(n)-a(m)*s0(n-m+1);
        end
    end
end

s2 = pmApply(length_of_s, freq, createSine(4,length_of_s), 0.1);
env = ADSR(floor(44100*0.045), 44100*0.05, 0.5, 44100*0.05, round(length_of_s*0.9), length_of_s);
s2 = pmApply(length_of_s, freq, s2, 1);
s2 = pmApply(length_of_s, freq*4, s2, 3);
s2 = pmApply(length_of_s, freq, s2, 1);
s2=s2.*env;
fc= 18000/length_of_s;
Q=1/sqrt(2);
[b,a]=IIR_LPF(fc,Q);
I=2;
J=2;
for n=1:length_of_s,
    for m=1:J+1,
        if n-m+1>0
            s3(n)=s3(n)+b(m)*s2(n-m+1);
        end
    end
    for m=2:I+1,
        if n-m+1>0
            s3(n)=s3(n)-a(m)*s3(n-m+1);
        end
    end
end


s0 = s0*0.2+s2*0.8;

y = s0;

end

実用

これで数小節だけ書いてみた結果が以下のレポジトリにあります。これを作るだけでもかなり時間がかかりました。

https://github.com/Clpsplug/Octave_FM_Insts

さすがになっている音を音階や長さごとに1つ1つ音を作って最後に合計しようとすると組み合わせが爆発するので、 配列で楽譜をかけるようにしてみました。 bass_melody.mというファイルが音を生成しているのですが、最初に同じようなパターンが並びまくっているのはそれです。

Octaveでこれを実行する際、オーディオがならせる状態でビルド・インストールするようにしてください。homebrewでインストールするときは次のコマンドです。

brew install science/octave --with-portaudio

まとめ

やはりFM音源での音作りは難しいです。しかもOctaveはMATLABに比べて遅い3ので結果がなかなか出ません。
MIDIファイルを直接読んで生成しようかと思いましたが、まずは処理系の効率向上が先のようです。
それができたら一般化も行っていきたいところです。

明日の担当は@RIRA氏です!


  1. 北大ではCentOSとWindowsなどが入ったPCがあります ネットワークブートで選択できるようになっているようです 

  2. sensitivityを極端に高くした場合はこの限りではない。思いっきり発散し、いろんなところにスペクトルが見られるようになる 

  3. ちなみに、OctaveはDAWでもないのによく落ちます。 

Only registered users are allowed to comment on the entry. Please log in to comment.