正弦波の合成 - 人工知能に関する断創録

Pythonで音声信号処理（2011/05/14）

今までは、既存の音声ファイルを再生したり、波形を見たりしてきましたが、今回は、音の波の基本となる正弦波や正弦波を合成して作れる三角波、矩形波、ノコギリ波などを自前で作って音を鳴らしてみたいと思います。正弦波の基本式は、

$s(t) = A \sin(2 \pi f_0 t)$

です。ここで、Aは振幅、f0は基本周波数、tは時間です。今回は、コンピュータで音を扱うので時間を離散化します。

$t = n / f_s$

ここで、fsはサンプリング周波数（1秒間のサンプル数）、nはサンプルのインデックスです。これを先の正弦波の式に代入すると

$s(n) = A \sin(2 \pi f_0 n / fs)$

となって、tの式をnの式に変換できます。この式を使うと、nサンプル目の波の値が計算できます。

正弦波を作る

まずは、ドレミファソラシドの各音階の正弦波を作って鳴らしてみます。正弦波はこんな形。

f:id:aidiary:20110607204946p:plain

音だと

ドは262Hz、ラは440Hzというように周波数が決められています。周波数の間には、十二平均律という数学的に興味深い構造があります。「虚数の情緒」っていう本で初めて知りました。音楽の授業でもやったことあるかな？

#coding: utf-8
import wave
import struct
import numpy as np
from pylab import *

def createSineWave (A, f0, fs, length):
    """振幅A、基本周波数f0、サンプリング周波数 fs、
    長さlength秒の正弦波を作成して返す"""
    data = []
    # [-1.0, 1.0]の小数値が入った波を作成
    for n in arange(length * fs):  # nはサンプルインデックス
        s = A * np.sin(2 * np.pi * f0 * n / fs)
        # 振幅が大きい時はクリッピング
        if s > 1.0:  s = 1.0
        if s < -1.0: s = -1.0
        data.append(s)
    # [-32768, 32767]の整数値に変換
    data = [int(x * 32767.0) for x in data]
#    plot(data[0:100]); show()
    # バイナリに変換
    data = struct.pack("h" * len(data), *data)  # listに*をつけると引数展開される
    return data

def play (data, fs, bit):
    import pyaudio
    # ストリームを開く
    p = pyaudio.PyAudio()
    stream = p.open(format=pyaudio.paInt16,
                    channels=1,
                    rate=int(fs),
                    output= True)
    # チャンク単位でストリームに出力し音声を再生
    chunk = 1024
    sp = 0  # 再生位置ポインタ
    buffer = data[sp:sp+chunk]
    while buffer != '':
        stream.write(buffer)
        sp = sp + chunk
        buffer = data[sp:sp+chunk]
    stream.close()
    p.terminate()

if __name__ == "__main__" :
    freqList = [262, 294, 330, 349, 392, 440, 494, 523]  # ドレミファソラシド
    for f in freqList:
        data = createSineWave(1.0, f, 8000.0, 1.0)
        play(data, 8000, 16)

createSineWave()は、振幅、基本周波数、サンプリング周波数、長さの4つの引数をとり、正弦波のデータを作成します。ここで作られるデータは、WAVEファイルと同じ形式のバイナリデータです。なので、WAVEファイルの再生（2011/5/15）で使ったpyaudioを使って音を鳴らすことができます。以前は、WAVEファイルのポインタからフレームを読み取っていましたが、今回は、createSineWave()で返されたバイナリデータを読み取ってストリームに流しています。メイン関数では、ドレミファソラシドの各周波数の波形をcreateSineWave()で作ってplay()で再生しています。dataをplotすれば波の形が実際に見られます。

合成波を作る

次に、正弦波を組み合わせた合成波を作ってみます。音楽には疎いんですが、和音（コード）というのがあるようなのでその音を作ってみました。作曲の仕方、最低限覚えることを読みました。どうやら和音ってのはドレミの中から三つの音を同時に鳴らした音のようです。というわけで正弦波を三つ重ね合わせて鳴らしてみます。波の形は下のような形です。

f:id:aidiary:20110607204947p:plain

音だと

#coding: utf-8
import wave
import struct
import numpy as np
from pylab import *

def createCombinedWave (A, freqList, fs, length):
    """freqListの正弦波を合成した波を返す"""
    data = []
    amp = float(A) / len(freqList)
    # [-1.0, 1.0]の小数値が入った波を作成
    for n in arange(length * fs):  # nはサンプルインデックス
        s = 0.0
        for f in freqList:
            s += amp * np.sin(2 * np.pi * f * n / fs)
        # 振幅が大きい時はクリッピング
        if s > 1.0:  s = 1.0
        if s < -1.0: s = -1.0
        data.append(s)
    # [-32768, 32767]の整数値に変換
    data = [int(x * 32767.0) for x in data]
    # バイナリに変換
    data = struct.pack("h" * len(data), *data)  # listに*をつけると引数展開される
    return data

def play (data, fs, bit):
    import pyaudio
    # ストリームを開く
    p = pyaudio.PyAudio()
    stream = p.open(format=pyaudio.paInt16,
                    channels=1,
                    rate=int(fs),
                    output= True)
    # チャンク単位でストリームに出力し音声を再生
    chunk = 1024
    sp = 0  # 再生位置ポインタ
    buffer = data[sp:sp+chunk]
    while buffer != '':
        stream.write(buffer)
        sp = sp + chunk
        buffer = data[sp:sp+chunk]
    stream.close()
    p.terminate()

if __name__ == "__main__" :
    # 和音（あってる？）
    chordList = [(262, 330, 392),  # C（ドミソ）
                 (294, 370, 440),  # D（レファ#ラ）
                 (330, 415, 494),  # E（ミソ#シ）
                 (349, 440, 523),  # F（ファラド）
                 (392, 494, 587),  # G（ソシレ）
                 (440, 554, 659),  # A（ラド#ミ）
                 (494, 622, 740)]  # B（シレ#ファ#）
    for freqList in chordList:
        data = createCombinedWave(1.0, freqList, 8000.0, 1.0)
        play(data, 8000, 16)

振幅は、[-1.0, 1.0]の範囲を超えるとクリッピングされてしまうので、合成波の各波の振幅は下げています。

三角波

次に、三角波っていうぎざぎざした波を作ってみます。こんな形です。

f:id:aidiary:20110607202717p:plain

この三角波は、先ほどやったみたいに多数の正弦波を重ね合わせることで作れます。Wikipediaの三角波（triangle wave）の項目を見ると、三角波の式は、

$x(t) = \frac{8}{\pi^2} \sum_{k=0}^\infty (-1)^k \frac{\sin((2k+1) \omega t)}{(2k+1)^2}$

ここで、 $\omega = 2 \pi f_0$ です。三角波は、無限個の正弦波の組み合わせであることがわかりますが、無限個は無理なので10個までにします。以下、三角波でドレミを鳴らすプログラムです。

#coding: utf-8
import wave
import struct
import numpy as np
from pylab import *

def createTriangleWave (A, f0, fs, length):
    data = []
    # [-1.0, 1.0]の小数値が入った波を作成
    for n in arange(length * fs):  # nはサンプルインデックス
        s = 0.0
        for k in range(0, 10):  # サンプルごとに10個のサイン波を重ね合わせ
            s += (-1)**k * (A / (2*k+1)**2) * np.sin((2*k+1) * 2 * np.pi * f0 * n / fs)
        # 振幅が大きい時はクリッピング
        if s > 1.0:  s = 1.0
        if s < -1.0: s = -1.0
        data.append(s)
    # [-32768, 32767]の整数値に変換
    data = [int(x * 32767.0) for x in data]
#    plot(data[0:100]); show()  # 波を描画
    # バイナリに変換
    data = struct.pack("h" * len(data), *data)  # listに*をつけると引数展開される
    return data

def play (data, fs, bit):
    import pyaudio
    # ストリームを開く
    p = pyaudio.PyAudio()
    stream = p.open(format=pyaudio.paInt16,
                    channels=1,
                    rate=int(fs),
                    output= True)
    # チャンク単位でストリームに出力し音声を再生
    chunk = 1024
    sp = 0  # 再生位置ポインタ
    buffer = data[sp:sp+chunk]
    while buffer != '':
        stream.write(buffer)
        sp = sp + chunk
        buffer = data[sp:sp+chunk]
    stream.close()
    p.terminate()

if __name__ == "__main__" :
    freqList = [262, 294, 330, 349, 392, 440, 494, 523]  # ドレミファソラシド
    for f in freqList:
        data = createTriangleWave(0.5, f, 8000.0, 1.0)
        play(data, 8000, 16)

振幅は、Aで指定できるようにしています。

矩形波

次は、矩形波、四角い波です。

f:id:aidiary:20110607204547p:plain

三角波と同じく無限個の正弦波の重ね合わせで作れますが、式がちょっと違います。同じくWikipediaの矩形波（Square wave）の項目を見ると、矩形波の式は、

$x(t) = \frac{4}{\pi} \sum_{k=1}^\infty \frac{\sin((2k-1) \pi f t)}{(2k-1)}$

同じく矩形波でドレミを鳴らしてみます。

#coding: utf-8
import wave
import struct
import numpy as np
from pylab import *

def createSquareWave (A, f0, fs, length):
    data = []
    # [-1.0, 1.0]の小数値が入った波を作成
    for n in arange(length * fs):  # nはサンプルインデックス
        s = 0.0
        for k in range(1, 10):
            s += (A / (2*k-1)) * np.sin((2*k-1) * 2 * np.pi * f0 * n / fs)
        # 振幅が大きい時はクリッピング
        if s > 1.0:  s = 1.0
        if s < -1.0: s = -1.0
        data.append(s)
    # [-32768, 32767]の整数値に変換
    data = [int(x * 32767.0) for x in data]
    # バイナリに変換
    data = struct.pack("h" * len(data), *data)  # listに*をつけると引数展開される
    return data

def play (data, fs, bit):
    import pyaudio
    # ストリームを開く
    p = pyaudio.PyAudio()
    stream = p.open(format=pyaudio.paInt16,
                    channels=1,
                    rate=int(fs),
                    output= True)
    # チャンク単位でストリームに出力し音声を再生
    chunk = 1024
    sp = 0  # 再生位置ポインタ
    buffer = data[sp:sp+chunk]
    while buffer != '':
        stream.write(buffer)
        sp = sp + chunk
        buffer = data[sp:sp+chunk]
    stream.close()
    p.terminate()

if __name__ == "__main__" :
    freqList = [262, 294, 330, 349, 392, 440, 494, 523]  # ドレミファソラシド
    for f in freqList:
        data = createSquareWave(0.5, f, 8000.0, 1.0)
        play(data, 8000, 16)

先の図のように四角い波を作るには100個ぐらい重ね合わせています。10個程度ではちょっとぎざぎざが残って完全な四角になりませんでした。

のこぎり波

最後は、のこぎり波です。こんな形。

f:id:aidiary:20110607204948p:plain

これも三角波、矩形波と同じく無限の正弦波の重ね合わせで作れます。Wikipediaののこぎり波（Sawtooth wave）の項目を見ると、のこぎり波の式は、

$x(t) = \frac{2}{\pi} \sum_{k=1}^\infty \frac{\sin(2 \pi kft)}{k}$

です。こっちの方が簡単ですね。同じくのこぎり波でドレミを鳴らしてみます。

#coding: utf-8
import wave
import struct
import numpy as np
from pylab import *

def createSawtoothWave (A, f0, fs, length):
    data = []
    # [-1.0, 1.0]の小数値が入った波を作成
    for n in arange(length * fs):  # nはサンプルインデックス
        s = 0.0
        for k in range(1, 10):
            s += (A / k) * np.sin(2 * np.pi * k * f0 * n / fs)
        # 振幅が大きい時はクリッピング
        if s > 1.0:  s = 1.0
        if s < -1.0: s = -1.0
        data.append(s)
    # [-32768, 32767]の整数値に変換
    data = [int(x * 32767.0) for x in data]
    plot(data[0:100]); show()
    # バイナリに変換
    data = struct.pack("h" * len(data), *data)  # listに*をつけると引数展開される
    return data


def play (data, fs, bit):
    import pyaudio
    # ストリームを開く
    p = pyaudio.PyAudio()
    stream = p.open(format=pyaudio.paInt16,
                    channels=1,
                    rate=int(fs),
                    output= True)
    # チャンク単位でストリームに出力し音声を再生
    chunk = 1024
    sp = 0  # 再生位置ポインタ
    buffer = data[sp:sp+chunk]
    while buffer != '':
        stream.write(buffer)
        sp = sp + chunk
        buffer = data[sp:sp+chunk]
    stream.close()
    p.terminate()

if __name__ == "__main__" :
    freqList = [262, 294, 330, 349, 392, 440, 494, 523]  # ドレミファソラシド
    for f in freqList:
        data = createSawtoothWave(0.5, f, 8000.0, 1.0)
        play(data, 8000, 16)

ちょっとブザーっぽい感じがします。

WAVEファイルへ出力

最後に、作成したdataをWAVEファイルへ出力する関数を作って終わりにします。これは、waveモジュールを使うと簡単にできます。

def save(data, fs, bit, filename):
    """波形データをWAVEファイルへ出力"""
    wf = wave.open(filename, "w")
    wf.setnchannels(1)
    wf.setsampwidth(bit / 8)
    wf.setframerate(fs)
    wf.writeframes(data)
    wf.close()

下みたいな感じで使います。dataはバイナリですが、Pythonのおかげで文字列データと同じように扱えるので"+"で連結できたりします。今回、アップロードしたWAVEファイルは下のような感じで保存しました。

if __name__ == "__main__" :
    allData = ""
    freqList = [262, 294, 330, 349, 392, 440, 494, 523]  # ドレミファソラシド
    for f in freqList:
        data = createSineWave(1.0, f, 8000.0, 1.0)
        play(data, 8000, 16)
        allData += data
    save(allData, 8000, 16, "sine.wav")

そうそう、昔懐かしいファミコンのピコピコ音は、三角波と矩形波だけで作られていたそうです。限界を超えるテクニック（ファミコンサウンド編）。何か面白いですねー。