リアルタイムFM合成を考える。

最近になってようやくFM合成をどのように作ればいいのかが分かったので、リアルタイムで合成する時はどのように実装するのかを考える。
とりあえず目標はMAX/MSP上での実装。だって一番楽だもの。

では実装する前にFM音源について学びましょう。
まあ大体wikiを見たら分かるのですが、重要なところを抜き出したものが以下の通り。

• John Chowning氏が開発したFrequency Modulationを応用する音色合成方式で、現在はYAMAHAがライセンスを受け実用化(DX7, reface DXなど)
• 倍音減算方式のアナログシンセサイザーにはない複雑な倍音成分を持つ波形を生成することが可能
• 挙動はカオスであるため、パラメータ値の変動による倍音変化は予測し難い

こんな感じなのですが、どういった音色なのかは実際にYoutubeで「DX7」や「reface DX」などで検索して確認してみましょう。きっと氏家さんのプレイに魅了されることでしょう。

さて、では実際にDXなどで動いている数式は以下のとおりです。

数式にあるとおり基本形はsin波をsin波で変調しています。しかしこれをもっとsin波に限らず一般的な波形に当てはめるとこのような感じに。

Fc:キャリア波形, Fm:モジュレータ波形, A:キャリア振幅, C:キャリア周波数, β:変調指数, M:モジュレータ周波数, t:時刻

この数式が示しているのは、ある時刻tのモジュレータ波形の値だけずらした時刻t+Fm(t)のキャリア周波数を出力として出す、ということがわかります。
では簡単にPythonで組んでみて、波形を確認しましょう。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

n = 10000
x = np.linspace(0, 4*np.pi, num=n)

## generate waveforms
y1 = np.sin(x + np.sin(x))
y2 = np.sin(x + 5*np.sin(x))
y3 = np.sin(x + np.sin(5*x))

## plot
plt.subplot(311)
plt.plot(x, y1, label="y = sin(x + sin(x))")
plt.xlim(min(x), max(x))
plt.legend(loc='upper right')

plt.subplot(312)
plt.plot(x, y2, label="y = sin(x + 5sin(x))")
plt.xlim(min(x), max(x))
plt.legend(loc='upper right')

plt.subplot(313)
plt.plot(x, y3, label="y = sin(x + sin(5x))")
plt.xlim(min(x), max(x))
plt.legend(loc='upper right')

plt.show()

こういった感じで変調指数やモジュレータ周波数などを変えることで複雑な波形の合成が可能となります。それと同時に波形の変化の予想が難しく音作りが困難になるといったデメリットも生じています。これがFM音源の奥深さというやつですね。

さて、前置きはここまでにして重要なのはこれをリアルタイムで合成するという点です。
Max/MSPのオブジェクトとして用意することを考えるとこんな感じでしょうか？

このときcarrier, modulatorには任意の波形、modulation indexには任意の値(float)が与えられます。となると今まではsin波の計算だけで済んでいたものがそうともいかなくなりました。ではどうすればいいでしょう？
そこで先ほど、FM合成はある時刻tのモジュレータ波形の値だけずらした時刻t+Fm(t)のキャリア周波数を出力として出す、ということが分かりました。
つまり、モジュレータから受け取った値を無理やり時間軸に直して、キャリアにおけるその時刻の振幅を得ればいいわけです。

実際に検算してみましょう。

まずサンプリング周波数を44100[Hz]し、t1=5/44100と置きます。
このとき, sin(t1)=0.000113381255811, sin(t1 + sin(t1))=0.000226762509435となります。
ではモジュレータから受け取ったsin(t1)を時間軸に直したいのですが、波形は配列でインデックスでしか参照できません。なのでサンプリング周期で除算することでインデックスに変換します。よって、

sin(t1)*44100/2π = 4.99999998929

ここでインデックスが小数点になるので線形補完をつかっておおよその値を求めます。なのでここで得られた値をinteger=4, decimal=0.99999998929と置きます。
そこから、t2=t1+(integer/44100) = 9/44100 とまたその次のt3=t1+( (integer+1) /44100) = 10/44100におけるsinの値は、

sin(t2) = 0.000204086259043
sin(t3) = 0.000226762509678

ここから線形補完をつかってsin(t1 + sin(t1))の値を近似します。

sin(t3)-sin(t2) = 2.26762506355e-05
(sin(t3)-sin(t2))*decimal = 2.26762503926e-05
sin(t2) + (sin(t3)-sin(t2))*decimal = 0.000226762509435
sin(t1 + sin(t1))=0.000226762509435

よってうまく近似できていることが分かりました。実際にpythonで計算したところ誤差は-2.71050543121e-20だそうです。やったね。
ただサンプリング周波数を下げたり高周波成分になる程精度は下がりますが、そこは仕方ないので割り切り。線形補完の限界です。

では実際にこの方法で波形を近似していきます。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

n = 44100
x = np.linspace(0, 4*np.pi, num=2*n)

gap = x[1]

carrier = np.sin(x)
modulator = np.sin(x)

y1 = np.sin(x + np.sin(x))
y2 = np.zeros(n) ## output

## genarate wave
for i in range(n):
	fmt = modulator[i]/gap
	decimal, integer = np.modf(fmt)
	integer = int(integer)
	y2[i] = carrier[i+integer] + (carrier[i+1+integer]-carrier[i+integer])*decimal


plt.title("y = sin(x + sin(x))")
plt.plot(x,y1,label="correct")
plt.plot(x[:n],y2,'--',label="approximation")
plt.xlim(min(x),max(x)/2)
plt.legend(loc='upper right')
plt.show()

なかなかうまく近似できました。
当然ですが、モジュレーション関数を書き換えても同様の結果が得られます。

しかしお気付きの方もいると思いますが、この状態ではまだ大きな欠陥が存在します。
モジュレーション関数であるsin(x)が取る最大値は1です。このときFm(t)=44100になります。また振幅値を上げたらその分だけFm(t)が増加します。しかし基本的にリアルタイムで処理する際、未来の波形は取得できません。取れてもサンプリングフレーム分だけです。そのため、リアルタイムで処理するためにはあらかじめ大量のバッファを確保して、基本的に過去の波形を参照し続けるという方式でないと理論上不可能です。
まあ、おそらく大量のバッファには適当な波形を起動時に生成しているかゼロ・パディングでバッファをクリアしているのでしょう。そこらへんは、実際に実装するときに考えましょう。

以上。