Help us understand the problem. What is going on with this article?

FM音源について

More than 3 years have passed since last update.

FM音源についてです。

2オペレータのFM合成

練習

 簡単な2オペレータのFM合成について考えてみましょう。今回はNative InstrumentsのFM8を使用します。最初のお題はコレ。

fm8_test_1_1.png

 上のように設定すると"Waveform"に次のようなグラフが表示されます。

fm8_test_1_2.png

 このグラフの波形が本当に出力されているかと言われると微妙ですが、このグラフの波形を自前で描いてみましょう。数式で書くと次のようなイメージだと思います。

g_2(t_n) = \sin( \, 2\pi ( \, (\alpha_2 f + \beta_2 )t_n 
+ m_1 \cdot \sin( \, 2\pi ( \alpha_1 f + \beta_1 )t_n ) \, ) \, )

 ここで $f$は基本音の周波数、$ \{ \alpha_1, \alpha_2 \} $はレシオ、$ \{ \beta_1, \beta_2 \} $は オフセット周波数 、$m_1$は モジュレータ振幅 としました。更に次のように変形してみます。

\begin{eqnarray}
c_1 &=& \alpha_1 f + \beta_1  \\
c_2 &=& \alpha_2 f + \beta_2  \\
g_1(t_n) &=& \sin( \, 2\pi c_1 t_n  \, ) \\
g_2(t_n) &=& \sin( \, 2\pi ( c_2 t_n 
+ m_1 \cdot g_1(t_n) ) \, )
\end{eqnarray}

 $\{ c_1, c_2 \}$を キャリア周波数 、更に変調に係わる部分を$g_1$と置きました。揺らす方の$g_1$は モジュレータ で揺らされる方の$g_2$は キャリア と言います。また、モジュレータとキャリアのことをまとめて オペレータ と言ったりします。では実際に計算・・・と行きたいところですが、FMマトリクスに表示されている値と実際の値は異なるようです。調べてみたところ次のような感じでした。

fm_param.png

\begin{eqnarray}
m_{ij} &=&
\begin{cases}
u_{ij}^2 \, / \, 9600 \quad (i=j) \\
u_{ij}^2 \, / \, 4800 \quad (i \neq j )
\end{cases} \\

h_j &=& w_j^2 \, / \, 10000
\end{eqnarray}

 間違っている可能性大ですが、深入りはしないで先に進みましょう。ではGnuplot用のdatファイルを作成してみます。

fm_test_1.rb
# coding: utf-8

# モジュレータ振幅の変換
def convert_modulation_amplitude(x)
    return x*x/4800.0
end

f  = 260.0   # 基本音の周波数[Hz]
fs = 44100.0 # サンプリング周波数[Hz]

a1 = 1.45  # Ratio(オペレータ1)
a2 = 1.0   # Ratio(オペレータ2)
b1 = 200.0 # Offset[Hz](オペレータ1)
b2 = 0.0   # Offset[Hz](オペレータ2)

m12 = convert_modulation_amplitude(55.0) # モジュレータ振幅

c1 = a1*f + b1 # キャリア周波数[Hz](オペレータ1)
c2 = a2*f + b2 # キャリア周波数[Hz](オペレータ2)

dt1 = c1/fs # 刻み間隔[s](オペレータ1)
dt2 = c2/fs # 刻み間隔[s](オペレータ2)

t1 = 0.0 # 時刻(オペレータ1)
t2 = 0.0 # 時刻(オペレータ2)

(2*fs/f).truncate.times do
    # update operator
    g1 = Math.sin(2*Math::PI*(t1))
    g2 = Math.sin(2*Math::PI*(t2 + m12*g1))

    # output
    p g2

    # update time
    t1 += dt1
    t2 += dt2
end

 実際にVSTプラグインなどのリアルタイム処理では負荷を軽くするためにsin関数を使わずにウェーブテーブルを使用するなど工夫しますがここでは頑張らないことにします。g2の値をグラフに表示してみましょう。

Screenshot from 2014-05-23 12_12_10.png

 おお。大体あってそうです。

フィードバッグ

 FM合成でよく フィードバック という単語を聞きますが、具体的にはどのような処理を行っているのでしょうか? 次はコレを考えてみましょう。

fm8_test_2_1.png

fm8_test_2_2.png

 $g_1(t_n)$はこれから計算しようとしている値なので単純にフィードバックとして使用するのは難しそうです。そうすると1サンプル前の$g_1(t_{n-1})$が使用するのが手っ取り早そうです。

\begin{eqnarray}
g_1(t_n) &=& \sin( \, 2\pi ( c_1 t_n
&+& m_{11} \cdot g_1(t_{n-1} ) ) \, )  \\
g_2(t_n) &=& \sin( \, 2\pi( c_2 t_n 
&+& m_{21} \cdot g_1(t_n) ) \, )
\end{eqnarray}

 ではdatファイル作成します。

fm_test_2.rb
# coding: utf-8

# モジュレータ振幅の変換
def convert_modulation_amplitude(x)
    return x*x/4800.0
end

# モジュレータ振幅の変換(対角成分)
def convert_modulation_amplitude_diag(x)
    return x*x/9600.0
end

f  = 260.0   # 基本音の周波数[Hz]
fs = 44100.0 # サンプリング周波数[Hz]

a1 = 1.45  # Ratio(オペレータ1)
a2 = 1.0   # Ratio(オペレータ2)
b1 = 200.0 # Offset[Hz](オペレータ1)
b2 = 0.0   # Offset[Hz](オペレータ2)

m11 = convert_modulation_amplitude_diag(40.0) # モジュレータ振幅
m12 = convert_modulation_amplitude(55.0)      # モジュレータ振幅

c1 = a1*f + b1 # キャリア周波数[Hz](オペレータ1)
c2 = a2*f + b2 # キャリア周波数[Hz](オペレータ2)

dt1 = c1/fs # 刻み間隔[s](オペレータ1)
dt2 = c2/fs # 刻み間隔[s](オペレータ2)

g1 = 0.0 # オペレータ1
t1 = 0.0 # 時刻(オペレータ1)
t2 = 0.0 # 時刻(オペレータ2)

(2*fs/f).truncate.times do
    # update operator
    g1 = Math.sin(2*Math::PI*(t1 + m11*g1))
    g2 = Math.sin(2*Math::PI*(t2 + m12*g1))

    # output
    p g2

    # update time
    t1 += dt1
    t2 += dt2
end

 $g_2$の値をグラフに表示してみましょう。

Screenshot from 2014-05-23 12_12_21.png

 今回もまあまあ大丈夫そうですね。

拡張

 FM合成の式の拡張を考えてみましょう。

他のオペレータへのフィードバック

 フィードバックについてですが、$g_2$のフィードバックというのも考えられそうです。$g_1$のときと同じように1サンプル前の値 $g_2(t_{n-1})$ を使用します。

\begin{eqnarray}
g_1(t_n) &=& \sin( \, 2\pi ( c_1 t_n 
&+& m_{11} \cdot g_1(t_{n-1}) ) \, )  \\
g_2(t_n) &=& \sin( \, 2\pi ( c_2 t_n 
&+& m_{21} \cdot g_1(t_n) + m_{22} \cdot g_2(t_{n-1}) ) \, )
\end{eqnarray}

 ついでに$g_2$ から$g_1$へのフィードバックというのがあっても良さそうですね。

\begin{eqnarray}
g_1(t_n) &=& \sin( \, 2\pi ( c_1 t_n 
&+& m_{11} \cdot g_1(t_{n-1} ) &+& m_{12} \cdot g_2(t_{n-1}) ) \, )  \\
g_2(t_n) &=& \sin( \, 2\pi ( c_2 t_n 
&+& m_{21} \cdot g_1(t_n) &+& m_{22} \cdot g_2(t_{n-1}) ) \, )
\end{eqnarray}

他のオペレータ値の出力

 $g_2$だけでなく$g_1$も出力に使用してみるのもアリかもしれません。$\{ h_1, h_2 \}$をボリューム、$v_n$を出力値とすると次のように表現出来ると思います。

v_n = h_1 \cdot g_1(t_n) + h_2 \cdot g_2(t_n)

FM合成の行列表現

モジュレータ項の統一

 $g_1, g_2$ 式をもう一度見てみましょう。

\begin{eqnarray}
g_1(t_n) &=& \sin( \, 2\pi ( c_1 t_n 
&+& m_{11} \cdot g_1(t_{n-1} ) &+& m_{12} \cdot g_2(t_{n-1}) ) \, )  \\
g_2(t_n) &=& \sin( \, 2\pi ( c_2 t_n 
&+& m_{21} \cdot g_1(t_n) &+& m_{22} \cdot g_2(t_{n-1}) ) \, )
\end{eqnarray}

 モジュレータ項の$m_{21} \cdot g_1(t_n)$が仲間外れのような気がします。この時点で$g_1(t_n)$の値は求まっていますので使用して構わないのですが、変調に関しては1サンプル前で計算済みの値を使用するように統一してみます。

\begin{eqnarray}
g_1(t_n) &=& \sin( \, 2\pi ( c_1 t_n 
&+& m_{11} \cdot g_1(t_{n-1} ) &+& m_{12} \cdot g_2(t_{n-1}) ) \, ) \\
g_2(t_n) &=& \sin( \, 2\pi ( c_2 t_n 
&+& m_{21} \cdot g_1(t_{n-1}) &+& m_{22} \cdot g_2(t_{n-1}) ) \, )
\end{eqnarray}

 キャリアとモジュレータのことをオペレータと言いますが、どちらかと言うとオペレータ自身がキャリアになったりモジュレータになったりという表現の方がしっくりくるかもしれません。

行列表現

 結局何がしたいのかと言うと、上のように置くことによって次のように変形することが出来ます。
キャリア周波数を・・・

c =
\begin{pmatrix}
c_1 \\
c_2
\end{pmatrix}
=
\begin{pmatrix}
\alpha_1 \\
\alpha_2
\end{pmatrix}
f +
\begin{pmatrix}
\beta_1 \\
\beta_2
\end{pmatrix}

キャリア波形を・・・

\psi(x) = \psi
\begin{pmatrix}
x_1 \\ 
x_2 
\end{pmatrix}
=
\begin{pmatrix}
\sin(2 \pi x_1) \\
\sin(2 \pi x_2)
\end{pmatrix}

モジュレータ振幅を・・・

M =
\begin{pmatrix}
m_{11} & m_{12} \\
m_{21} & m_{22}
\end{pmatrix}

オペレータを・・・

g(t) =
\begin{pmatrix}
g_1(t) \\
g_2(t) \\
\end{pmatrix}

ボリュームを・・・

h =
\begin{pmatrix}
h_1 \\
h_2
\end{pmatrix}

と置けば

\begin{eqnarray}
g(t_n) &=& \psi( \, c t_n + M g(t_{n-1}) \, ) \\
v_n &=& ( \, g(t_n), \, h \, )
\end{eqnarray}

 おー、行列とベクトルの計算で表現することが出来ました。FMマトリクスがなぜ"FMマトリクス"と呼ばれている理由が解った気がします。

6オペレータのFM合成

任意のアルゴリズムに対するFM合成

 2オペレータで考えてきましたが6オペレータの場合(それ以上でもOK)について考えてみましょう。FM合成の線の繋ぎ方のことを アルゴリズム と言います。アルゴリズムが固定の場合はそれ用に処理を実装するのも良いですが、任意のFMマトリクスに対応するにはどのような実装にすれば良いでしょうか? 先程、FM合成を行列を使って表現しましたが、行列のサイズを2から6にするだけで行けそうな気がします。では適当の値で確認してみましょう。最後のお題はコレ。

fm8_test3_1.png

fm8_test3_2.png

 実装ですが、Rubyは標準でMatrixクラスとVectorクラスが用意されています。折角なので使ってみます。

fm_test_3.rb
# coding: utf-8
require "matrix"

# モジュレータ振幅の変換
def convert_modulation_amplitude(arr)
    arr.map.with_index do |row,i|
        row.map.with_index do |x,j|
            i==j ? (x*x/9600.0) : (x*x/4800.0)
        end
    end
end

# ボリュームの変換
def convert_volume(arr)
    arr.map do |x|
        x*x/10000.0
    end
end

f  = 260.0   # 基本音の周波数[Hz]
fs = 44100.0 # サンプリング周波数[Hz]

# Ratio
a = Vector[ 2.0000,
            0.9995,
            1.0005,
            0.9995,
            1.0005,
           11.0000]

# Offset[Hz]
b = Vector[200.00,
             2.00,
            -2.00,
            -2.00,
             2.00,
             0.00]

# モジュレータ振幅
diplay_fm_mat = [[ 0,20, 0, 0, 0, 0],
                 [ 0, 0, 0, 0,16, 0],
                 [ 0, 0,33,44, 0, 0],
                 [ 0,61, 0, 0, 0, 0],
                 [38, 0, 0, 0, 0, 0],
                 [ 0, 0, 0, 0,64,28]]
fm_mat = Matrix.rows(convert_modulation_amplitude(diplay_fm_mat))

# ボリューム
display_volume = [ 0, 0,85, 0, 0,50]
volume = Vector.elements(convert_volume(display_volume))

# キャリア周波数[Hz]
c = a*f + b

# キャリア波形
carrier_waveform = ->(vec){
    Vector[Math.sin(2*Math::PI*vec[0]),
           Math.sin(2*Math::PI*vec[1]),
           Math.sin(2*Math::PI*vec[2]),
           Math.sin(2*Math::PI*vec[3]),
           Math.sin(2*Math::PI*vec[4]),
           Math.sin(2*Math::PI*vec[5])]
}

dt = c/fs                    # 刻み間隔[s]
g = Vector[0, 0, 0, 0, 0, 0] # オペレータ
t = Vector[0, 0, 0, 0, 0, 0] # 時刻

(2*fs/f).truncate.times do
    # update operator
    g = carrier_waveform[t + fm_mat*g]

    # output
    p g.inner_product(volume)

    # update time
    t += dt
end

 で結果は・・・

Screenshot from 2014-05-23 13_29_41.png

 今回もまあまあ大丈夫ということで。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした