リアルタイム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)が増加します。しかし基本的にリアルタイムで処理する際、未来の波形は取得できません。取れてもサンプリングフレーム分だけです。そのため、リアルタイムで処理するためにはあらかじめ大量のバッファを確保して、基本的に過去の波形を参照し続けるという方式でないと理論上不可能です。
まあ、おそらく大量のバッファには適当な波形を起動時に生成しているかゼロ・パディングでバッファをクリアしているのでしょう。そこらへんは、実際に実装するときに考えましょう。
以上。
AudioUnitでFuzzFaceを作る。その１：ダイオードシミュレーション
FuzzFaceって知っていますか？
↓これですね
モリダイラ楽器から引用
Jimi Hendrixも使っていたりする超有名なエフェクターですね。
これが欲しいのですが、お金がないです。
自作するにも環境がないです。あとパーツを買う元気もないです。
そこでAudioUnitでこれを実装しようという話です。プログラム上で組めればタダですから。
まあニッチな内容なので誰も興味ないでしょうが、自分用のメモとしてここに残しておきます。
cent ↔ Hzの変換式
最近、論文を読んでたら結構な割合で周波数軸をcent表記している人が多い。
よくよく考えたらチューナーなどでも合わせる対象の周波数に対してどれだけずれているかをcent表記で示していることが多い。(代表例:GarageBandのチューナーの単位もcent)
"ではcentとはなんだ？"といますと、wikipediaでは以下のように書かれています。
セントは2つの周波数の比である。平均律の半音の間隔は100セントと定義される。
https://ja.wikipedia.org/wiki/%E3%82%BB%E3%83%B3%E3%83%88_%28%E9%9F%B3%E6%A5%BD%29
オクターヴ（2つの音の周波数比が2:1）は12半音であり、1200セントである。
なんとなく分かるけど実際にHzからcentに変換するにはどう計算するのか分からなかったので、今回それについて考えてみます。
続きを読むAudioUnits パラメータの単位設定
AudioUnitsを何個か作ったなかで情報が少なくて面倒臭くなることが幾つかあったので、とりあえず記事にまとめてみます。
まあ、これぐらいはすぐに調べたら分かるような内容ですが一箇所にあるほうがよさそうなので。
続きを読む
AudioUnitEffectExample読み解き #その1
こういうソース解説の記事って自分がよくお世話になったりしたので自分で書いてみるけど、基本は自分の自己満足兼備忘録。
今回の対象は、以前の記事で落としたであろうAudioUnitのサンプルコード群の中にAudioUnitEffectExample。
本来の目的はエフェクタな訳で、他のAudioUnitGeneratorExample(PinkNoise生成)やAudioUnitInstrumentExample(おそらくMIDI制御のsin波シンセ)には目もくれず、これを弄ろう。
最初、StarterAudioUnitExample(Tremolo)を弄ろうと思ったんだけど、何かぱっとみ書き方がややこしい
(読むのが面倒くせえ)ので、見た感じがすっきりしているFilterDemoさんを見ていこうという感じ。
あとCocoaGUIが用意されているので要らない場合はCocoaGUI関係のソースコードを消せば、StarterAudioUnitExample(Tremolo)のような初期のパラメータ設定ウィンドウが表示されるので、今回はそうします。
そのため、実はFilter.hをゴミ箱にポイーしても大丈夫です。ただでさえ中身が少ないのによく見たら全部CocoaGUIに関するものだったので。
ここから内部処理に触れていきますが、著者はC++を触ったことある程度の人間ですので何か間違いや勘違いなどがあると思われますが、その時はコメントか何かで報告を……。
開発環境 : Yosemite(10.10.3), Xcode 6.3
PyAudio環境構築
将来的に音響解析などに手を出したいと思ったのは良かったのだが、それを行う言語をどうしようと悩みあげた挙句、pythonさんを利用することにした。理由はmatlabとかお金かかるし、最近はpythonさんが人気あるとか、そんな感じ。
ただ決めたのは良かったのだけど、PyAudioの環境構築に中々手間取ったので備忘録。
追記(15/05/25):
Anaconda 2.2.0の中に既にpyaudioがあるもよう。そのため色々と中身を変更しました。
開発環境 : Yosemite(10.10.3), Python 2.7.8 :: Anaconda 2.1.0 → Python 2.7.9 :: Anaconda 2.2.0 (x86_64)
