サウンドプログラミング2

科目名：　メディアネットワーク実験IIA（2022年～）
対象：　メディアネットワークコース3年目
日時：　10月30日（水） - 31日（木）13:00～18:00
場所：　M棟1階計算機室
レポート提出締切：　11月6日（水）13:00
レポート提出先：　情報エレクトロニクス棟6階6-08
連絡先：　青木 直史（Tel: 706-6532）（E-mail: naofumi.aoki@ist.hokudai.ac.jp）

目的

　音はマルチメディアコンテンツを構成する重要な要素である．本実験は，Pythonによるプログラミングを通して，サウンド処理に対する理解を深めることを目的としている．

1．はじめに

　本実験は，Jupyter Notebookを利用し，ブラウザを使ってプログラムを実行しながら進めるものとする．なお，音を確認する場合は各自のイヤフォンまたはヘッドフォンを使うこと．

(1) Jupyter Notebookの起動

Anaconda Promptを開く．

jupyter notebook
と入力して実行．

(2) ブラウザ（Chrome推奨）を利用してプログラミングを行う．

「New」ボタンをクリックして「Python3」を選択．

音ファイルや画像ファイルを利用するときは「Upload」ボタンをクリックし，必要なファイルをJupyter Notebookにアップロードすること．

　サウンドプログラミングの前提として、以下の動画を参考に、各自、音響学について理解を深めること。

音響学講義「第1回 サイン波」

音響学講義「第2回 サイン波の重ね合わせ」

音響学講義「第3回 周波数分析」

音響学講義「第4回 音の性質」

音響学講義「第5回 音声」

サウンドプログラミング1（メディアネットワーク実験IA）の手引書

サウンドプログラミング1

2．自動演奏

　「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．カノンの演奏が聞こえてくることを確かめてください．

(1)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def sine_wave(fs, f, a, duration):
    length_of_s = int(fs * duration)
    s = np.zeros(length_of_s)
    for n in range(length_of_s):
        s[n] = np.sin(2 * np.pi * f * n / fs)

    for n in range(int(fs * 0.01)):
        s[n] *= n / (fs * 0.01)
        s[length_of_s - n - 1] *= n / (fs * 0.01)

    gain = a / np.max(np.abs(s))
    s *= gain
    return s

(3)

score = np.array([[1, 2, 659.26, 0.8, 1],
                  [1, 3, 587.33, 0.8, 1],
                  [1, 4, 523.25, 0.8, 1],
                  [1, 5, 493.88, 0.8, 1],
                  [1, 6, 440.00, 0.8, 1],
                  [1, 7, 392.00, 0.8, 1],
                  [1, 8, 440.00, 0.8, 1],
                  [1, 9, 493.88, 0.8, 1],
                  [2, 2, 261.63, 0.8, 1],
                  [2, 3, 196.00, 0.8, 1],
                  [2, 4, 220.00, 0.8, 1],
                  [2, 5, 164.81, 0.8, 1],
                  [2, 6, 174.61, 0.8, 1],
                  [2, 7, 130.81, 0.8, 1],
                  [2, 8, 174.61, 0.8, 1],
                  [2, 9, 196.00, 0.8, 1]])

number_of_track = 2
number_of_note = score.shape[0]

(4)

fs = 44100
length_of_s = int(fs * 12)
track = np.zeros((length_of_s, number_of_track))
s = np.zeros(length_of_s)

(5)

for i in range(number_of_note):
    j = int(score[i, 0] - 1)
    onset = score[i, 1]
    f = score[i, 2]
    a = score[i, 3]
    duration = score[i, 4]
    x = sine_wave(fs, f, a, duration)
    offset = int(fs * onset)
    length_of_x = len(x)
    for n in range(length_of_x):
        track[offset + n, j] += x[n]

(6)

for j in range(number_of_track):
    for n in range(length_of_s):
        s[n] += track[n, j]

(7)

master_volume = 0.5
s /= np.max(np.abs(s))
s *= master_volume

(8)

for n in range(length_of_s):
    s[n] = (s[n] + 1.0) / 2.0 * 65536.0
    if s[n] > 65535.0:
        s[n] = 65535.0
    elif s[n] < 0.0:
        s[n] = 0.0;
    s[n] = (s[n] + 0.5) - 32768

wavfile.write('p1.wav', fs, s.astype(np.int16))

(9)

Audio('p1.wav')

3．楽譜データをMIDIのパラメータで書き換える

　MIDIは、音の高さをノートナンバー，音の大きさをベロシティ，音の長さをゲートタイムによって定義している．「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．カノンの演奏が聞こえてくることを確かめてください．

(1)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def sine_wave(fs, note_number, velocity, gate):
    length_of_s = int(fs * gate)
    s = np.zeros(length_of_s)
    f = 440 * np.power(2, (note_number - 69) / 12)
    for n in range(length_of_s):
        s[n] = np.sin(2 * np.pi * f * n / fs)

    for n in range(int(fs * 0.01)):
        s[n] *= n / (fs * 0.01)
        s[length_of_s - n - 1] *= n / (fs * 0.01)

    gain = velocity / 127 / np.max(np.abs(s))
    s *= gain
    return s

(3)

score = np.array([[1, 1920, 76, 100, 960],
                  [1, 2880, 74, 100, 960],
                  [1, 3840, 72, 100, 960],
                  [1, 4800, 71, 100, 960],
                  [1, 5760, 69, 100, 960],
                  [1, 6720, 67, 100, 960],
                  [1, 7680, 69, 100, 960],
                  [1, 8640, 71, 100, 960],
                  [2, 1920, 60, 100, 960],
                  [2, 2880, 55, 100, 960],
                  [2, 3840, 57, 100, 960],
                  [2, 4800, 52, 100, 960],
                  [2, 5760, 53, 100, 960],
                  [2, 6720, 48, 100, 960],
                  [2, 7680, 53, 100, 960],
                  [2, 8640, 55, 100, 960]])

division = 480
tempo = 120
number_of_track = 2
end_of_track = 10
number_of_note = score.shape[0]

(4)

fs = 44100
length_of_s = int(fs * (end_of_track + 2))
track = np.zeros((length_of_s, number_of_track))
s = np.zeros(length_of_s)

(5)

for i in range(number_of_note):
    j = int(score[i, 0] - 1)
    onset = (score[i, 1] / division) * (60 / tempo)
    note_number = score[i, 2]
    velocity = score[i, 3]
    gate = (score[i, 4] / division) * (60 / tempo)
    x = sine_wave(fs, note_number, velocity, gate)
    offset = int(fs * onset)
    length_of_x = len(x)
    for n in range(length_of_x):
        track[offset + n, j] += x[n]

(6)

for j in range(number_of_track):
    for n in range(length_of_s):
        s[n] += track[n, j]

(7)

master_volume = 0.5
s /= np.max(np.abs(s))
s *= master_volume

(8)

for n in range(length_of_s):
    s[n] = (s[n] + 1.0) / 2.0 * 65536.0
    if s[n] > 65535.0:
        s[n] = 65535.0
    elif s[n] < 0.0:
        s[n] = 0.0;
    s[n] = (s[n] + 0.5) - 32768

wavfile.write('p2.wav', fs, s.astype(np.int16))

(9)

Audio('p2.wav')

4．加算合成1（鉄琴）

　同じ曲でも，楽器を取り替えると，雰囲気は一変する．「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．鉄琴の音色でカノンの演奏が聞こえてくることを確かめてください．音の時間変化をADSR関数を使って定義していることに注意してください．

(1)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def ADSR(fs, A, D, S, R, gate, duration):
    A = int(fs * A)
    D = int(fs * D)
    R = int(fs * R)
    gate = int(fs * gate)
    duration = int(fs * duration)
    e = np.zeros(duration)
    if A != 0:
        for n in range(A):
            e[n] = 1.0 - np.exp(-5.0 * n / A)

    if D != 0:
        for n in range(A, gate):
            e[n] = S + (1.0 - S) * np.exp(-5.0 * (n - A) / D)

    else:
        for n in range(A, gate):
            e[n] = S

    if R != 0:
        for n in range(gate, duration):
            e[n]= e[gate - 1] * np.exp(-5.0 * (n - gate + 1) / R)

    return e

(3)

def glockenspiel(fs, note_number, velocity, gate):
    f0 = 440 * np.power(2, (note_number - 69) / 12)

    number_of_partial = 5

    VCO_A = [0, 0, 0, 0, 0]
    VCO_D = [0, 0, 0, 0, 0]
    VCO_S = [1, 1, 1, 1, 1]
    VCO_R = [0, 0, 0, 0, 0]
    VCO_offset = [f0 * 1, f0 * 2.76, f0 * 5.40, f0 * 8.93, f0 * 13.32]
    VCO_depth = [0, 0, 0, 0, 0]

    VCA_A = [0, 0, 0, 0, 0]
    VCA_D = [2, 0.5, 0.2, 0.2, 0.1]
    VCA_S = [0, 0, 0, 0, 0]
    VCA_R = [2, 0.5, 0.2, 0.2, 0.1]
    VCA_offset = [0, 0, 0, 0, 0]
    VCA_depth = [1, 0.5, 0.4, 0.4, 0.2]

    duration = 2

    length_of_s = int(fs * duration)
    s = np.zeros(length_of_s)

    for i in range(number_of_partial):
        vco = ADSR(fs, VCO_A[i], VCO_D[i], VCO_S[i], VCO_R[i], gate, duration)
        for n in range(length_of_s):
            vco[n] = VCO_offset[i] + vco[n] * VCO_depth[i];

        if np.max(vco) < fs / 2:
            x = np.zeros(length_of_s)
            t = 0;
            for n in range(length_of_s):
                x[n] = np.sin(2 * np.pi * t)
                delta = vco[n] / fs
                t += delta
                if t >= 1:
                    t -= 1

            vca = ADSR(fs, VCA_A[i], VCA_D[i], VCA_S[i], VCA_R[i], gate, duration)
            for n in range(length_of_s):
                vca[n] = VCA_offset[i] + vca[n] * VCA_depth[i];

            for n in range(length_of_s):
                s[n] += x[n] * vca[n]

    gain = velocity / 127 / np.max(np.abs(s))
    s *= gain

    return s

(4)

score = np.array([[1, 1920, 76, 100, 960],
                  [1, 2880, 74, 100, 960],
                  [1, 3840, 72, 100, 960],
                  [1, 4800, 71, 100, 960],
                  [1, 5760, 69, 100, 960],
                  [1, 6720, 67, 100, 960],
                  [1, 7680, 69, 100, 960],
                  [1, 8640, 71, 100, 960],
                  [2, 1920, 60, 100, 960],
                  [2, 2880, 55, 100, 960],
                  [2, 3840, 57, 100, 960],
                  [2, 4800, 52, 100, 960],
                  [2, 5760, 53, 100, 960],
                  [2, 6720, 48, 100, 960],
                  [2, 7680, 53, 100, 960],
                  [2, 8640, 55, 100, 960]])

division = 480
tempo = 120
number_of_track = 2
end_of_track = 10
number_of_note = score.shape[0]

(5)

fs = 44100
length_of_s = int(fs * (end_of_track + 2))
track = np.zeros((length_of_s, number_of_track))
s = np.zeros(length_of_s)

(6)

for i in range(number_of_note):
    j = int(score[i, 0] - 1)
    onset = (score[i, 1] / division) * (60 / tempo)
    note_number = score[i, 2]
    velocity = score[i, 3]
    gate = (score[i, 4] / division) * (60 / tempo)
    x = glockenspiel(fs, note_number, velocity, gate)
    offset = int(fs * onset)
    length_of_x = len(x)
    for n in range(length_of_x):
        track[offset + n, j] += x[n]

(7)

for j in range(number_of_track):
    for n in range(length_of_s):
        s[n] += track[n, j]

(8)

master_volume = 0.5
s /= np.max(np.abs(s))
s *= master_volume

(9)

for n in range(length_of_s):
    s[n] = (s[n] + 1.0) / 2.0 * 65536.0
    if s[n] > 65535.0:
        s[n] = 65535.0
    elif s[n] < 0.0:
        s[n] = 0.0;
    s[n] = (s[n] + 0.5) - 32768

wavfile.write('p3.wav', fs, s.astype(np.int16))

(10)

Audio('p3.wav')

5．加算合成2（パイプオルガン）

　「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．パイプオルガンの音色でカノンの演奏が聞こえてくることを確かめてください．残響音をreverb関数を使って定義していることに注意してください．

(1)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def ADSR(fs, A, D, S, R, gate, duration):
    A = int(fs * A)
    D = int(fs * D)
    R = int(fs * R)
    gate = int(fs * gate)
    duration = int(fs * duration)
    e = np.zeros(duration)
    if A != 0:
        for n in range(A):
            e[n] = 1.0 - np.exp(-5.0 * n / A)

    if D != 0:
        for n in range(A, gate):
            e[n] = S + (1.0 - S) * np.exp(-5.0 * (n - A) / D)

    else:
        for n in range(A, gate):
            e[n] = S

    if R != 0:
        for n in range(gate, duration):
            e[n]= e[gate - 1] * np.exp(-5.0 * (n - gate + 1) / R)

    return e

(3)

def pipe_organ(fs, note_number, velocity, gate):
    f0 = 440 * np.power(2, (note_number - 69) / 12)

    number_of_partial = 16

    VCO_A = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    VCO_D = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    VCO_S = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    VCO_R = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    VCO_offset = [f0 * 1, f0 * 2, f0 * 3, f0 * 4, f0 * 5, f0 * 6, f0 * 7, f0 * 8, f0 * 9, f0 * 10, f0 * 11, f0 * 12, f0 * 13, f0 * 14, f0 * 15, f0 * 16]
    VCO_depth = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

    VCA_A = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
    VCA_D = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    VCA_S = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    VCA_R = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
    VCA_offset = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    VCA_depth = [1, 1, 1, 1, 0.8, 0.8, 0.8, 0.8, 0.5, 0.5, 0.5, 0.5, 0.3, 0.3, 0.3, 0.3]

    duration = gate + 0.1

    length_of_s = int(fs * duration)
    s = np.zeros(length_of_s)

    for i in range(number_of_partial):
        vco = ADSR(fs, VCO_A[i], VCO_D[i], VCO_S[i], VCO_R[i], gate, duration)
        for n in range(length_of_s):
            vco[n] = VCO_offset[i] + vco[n] * VCO_depth[i];

        if np.max(vco) < fs / 2:
            x = np.zeros(length_of_s)
            t = 0;
            for n in range(length_of_s):
                x[n] = np.sin(2 * np.pi * t)
                delta = vco[n] / fs
                t += delta
                if t >= 1:
                    t -= 1

            vca = ADSR(fs, VCA_A[i], VCA_D[i], VCA_S[i], VCA_R[i], gate, duration)
            for n in range(length_of_s):
                vca[n] = VCA_offset[i] + vca[n] * VCA_depth[i];

            for n in range(length_of_s):
                s[n] += x[n] * vca[n]

    gain = velocity / 127 / np.max(np.abs(s))
    s *= gain

    return s

(4)

def reverb(fs, x):
    length_of_x = len(x)

    d1 = int(fs * 0.03985)
    g1 = 0.871402
    u1 = np.zeros(length_of_x)
    for n in range(length_of_x):
        if n - d1 >= 0:
            u1[n] = x[n - d1] + g1 * u1[n - d1]

    d2 = int(fs * 0.03610)
    g2 = 0.882762
    u2 = np.zeros(length_of_x)
    for n in range(length_of_x):
        if n - d2 >= 0:
            u2[n] = x[n - d2] + g2 * u2[n - d2]

    d3 = int(fs * 0.03327)
    g3 = 0.891443
    u3 = np.zeros(length_of_x)
    for n in range(length_of_x):
        if n - d3 >= 0:
            u3[n] = x[n - d3] + g3 * u3[n - d3]

    d4 = int(fs * 0.03015)
    g4 = 0.901117
    u4 = np.zeros(length_of_x)
    for n in range(length_of_x):
        if n - d4 >= 0:
            u4[n] = x[n - d4] + g4 * u4[n - d4]

    v1 = np.zeros(length_of_x)
    for n in range(length_of_x):
        v1[n] = u1[n] + u2[n] + u3[n] + u4[n]

    d5 = int(fs * 0.005)
    g5 = 0.7
    u5 = np.zeros(length_of_x)
    v2 = np.zeros(length_of_x)
    for n in range(length_of_x):
        if n - d5 >= 0:
            u5[n] = v1[n - d5] + g5 * u5[n - d5]

        v2[n] = u5[n] - g5 * (v1[n] + g5 * u5[n])

    d6 = int(fs * 0.0017)
    g6 = 0.7
    u6 = np.zeros(length_of_x)
    y = np.zeros(length_of_x)
    for n in range(length_of_x):
        if n - d6 >= 0:
            u6[n] = v2[n - d6] + g6 * u6[n - d6]

        y[n] = u6[n] - g6 * (v2[n] + g6 * u6[n])

    return y

(5)

score = np.array([[1, 1920, 76, 100, 960],
                  [1, 2880, 74, 100, 960],
                  [1, 3840, 72, 100, 960],
                  [1, 4800, 71, 100, 960],
                  [1, 5760, 69, 100, 960],
                  [1, 6720, 67, 100, 960],
                  [1, 7680, 69, 100, 960],
                  [1, 8640, 71, 100, 960],
                  [2, 1920, 60, 100, 960],
                  [2, 2880, 55, 100, 960],
                  [2, 3840, 57, 100, 960],
                  [2, 4800, 52, 100, 960],
                  [2, 5760, 53, 100, 960],
                  [2, 6720, 48, 100, 960],
                  [2, 7680, 53, 100, 960],
                  [2, 8640, 55, 100, 960]])

division = 480
tempo = 120
number_of_track = 2
end_of_track = 10
number_of_note = score.shape[0]

(6)

fs = 44100
length_of_s = int(fs * (end_of_track + 2))
track = np.zeros((length_of_s, number_of_track))
s = np.zeros(length_of_s)

(7)

for i in range(number_of_note):
    j = int(score[i, 0] - 1)
    onset = (score[i, 1] / division) * (60 / tempo)
    note_number = score[i, 2]
    velocity = score[i, 3]
    gate = (score[i, 4] / division) * (60 / tempo)
    x = pipe_organ(fs, note_number, velocity, gate)
    offset = int(fs * onset)
    length_of_x = len(x)
    for n in range(length_of_x):
        track[offset + n, j] += x[n]

(8)

for j in range(number_of_track):
    for n in range(length_of_s):
        s[n] += track[n, j]

(9)

s = reverb(fs, s)

(10)

master_volume = 0.5
s /= np.max(np.abs(s))
s *= master_volume

(11)

for n in range(length_of_s):
    s[n] = (s[n] + 1.0) / 2.0 * 65536.0
    if s[n] > 65535.0:
        s[n] = 65535.0
    elif s[n] < 0.0:
        s[n] = 0.0;
    s[n] = (s[n] + 0.5) - 32768

wavfile.write('p4.wav', fs, s.astype(np.int16))

(12)

Audio('p4.wav')

6．フィルタ

　フィルタは、特定の帯域の周波数成分だけを選択的に通過させるふるいである．「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．高域の周波数成分がカットされることを確認しなさい．

(1)

%matplotlib inline
import numpy as np
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def LPF(fs, fc, Q):
    fc /= fs
    fc = np.tan(np.pi * fc) / (2.0 * np.pi)
    a = np.zeros(3)
    b = np.zeros(3)
    a[0] = 1.0 + 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc
    a[1] = (8.0 * np.pi * np.pi * fc * fc - 2.0) / a[0]
    a[2] = (1.0 - 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc) / a[0]
    b[0] = 4.0 * np.pi * np.pi * fc * fc / a[0]
    b[1] = 8.0 * np.pi * np.pi * fc * fc / a[0]
    b[2] = 4.0 * np.pi * np.pi * fc * fc / a[0]
    a[0] = 1.0
    return a, b

(3)

fs = 8000

(4)

length_of_s = int(fs * 1)
s0 = np.zeros(length_of_s)

(5)

for i in range(1, 9):
    for n in range(length_of_s):
        s0[n] += 1.0 * np.sin(2 * np.pi * 440 * i * n / fs)

gain = 0.5 / np.max(np.abs(s0))
s0 *= gain

(6)

for n in range(length_of_s):
    s0[n] = (s0[n] + 1.0) / 2.0 * 65536.0
    if s0[n] > 65535.0:
        s0[n] = 65535.0
    elif s0[n] < 0.0:
        s0[n] = 0.0;
    s0[n] = (s0[n] + 0.5) - 32768

wavfile.write('pulse_train.wav', fs, s0.astype(np.int16))

(7)

Audio('pulse_train.wav')

(8)

s1 = np.zeros(length_of_s)

(9)

fc = 880
Q = 1 / np.sqrt(2)
a, b = LPF(fs, fc, Q)
for n in range(length_of_s):
    for m in range(0, 3):
        if n - m >= 0:
            s1[n] += b[m] * s0[n - m]

    for m in range(1, 3):
        if n - m >= 0:
            s1[n] += -a[m] * s1[n - m]

(10)

master_volume = 0.5
s1 /= np.max(np.abs(s1))
s1 *= master_volume

(11)

for n in range(length_of_s):
    s1[n] = (s1[n] + 1.0) / 2.0 * 65536.0
    if s1[n] > 65535.0:
        s1[n] = 65535.0
    elif s1[n] < 0.0:
        s1[n] = 0.0;
    s1[n] = (s1[n] + 0.5) - 32768

wavfile.write('p5.wav', fs, s1.astype(np.int16))

(12)

Audio('p5.wav')

7．減算合成1（ハイハット）

　原音をフィルタで削ることで音色をつくる音響合成のテクニックを減算合成と呼ぶ。「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．白色雑音からハイハットの音をつくることができることを確認しなさい．

(1)

%matplotlib inline
import numpy as np
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def HPF(fs, fc, Q):
    fc /= fs
    fc = np.tan(np.pi * fc) / (2.0 * np.pi)
    a = np.zeros(3)
    b = np.zeros(3)
    a[0] = 1.0 + 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc;
    a[1] = (8.0 * np.pi * np.pi * fc * fc - 2.0) / a[0]
    a[2] = (1.0 - 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc) / a[0]
    b[0] = 1.0 / a[0]
    b[1] = -2.0 / a[0]
    b[2] = 1.0 / a[0]
    a[0] = 1.0
    return a, b

(3)

def ADSR(fs, A, D, S, R, gate, duration):
    A = int(fs * A)
    D = int(fs * D)
    R = int(fs * R)
    gate = int(fs * gate)
    duration = int(fs * duration)
    e = np.zeros(duration)
    if A != 0:
        for n in range(A):
            e[n] = 1.0 - np.exp(-5.0 * n / A)

    if D != 0:
        for n in range(A, gate):
            e[n] = S + (1.0 - S) * np.exp(-5.0 * (n - A) / D)

    else:
        for n in range(A, gate):
            e[n] = S

    if R != 0:
        for n in range(gate, duration):
            e[n]= e[gate - 1] * np.exp(-5.0 * (n - gate + 1) / R)

    return e

(4)

def hihat_cymbal_close(fs, velocity, gate):
    duration = 1

    length_of_s = int(fs * duration)
    s0 = np.zeros(length_of_s)

    np.random.seed(0)
    for n in range(length_of_s):
        s0[n] = (np.random.rand() * 2.0) - 1.0

    VCF_A = [0]
    VCF_D = [0]
    VCF_S = [1]
    VCF_R = [0]
    VCF_offset = [10000]
    VCF_depth=[0]

    vcf = ADSR(fs, VCF_A[0], VCF_D[0], VCF_S[0], VCF_R[0], gate, duration)
    for n in range(length_of_s):
        vcf[n] = VCF_offset[0] + vcf[n] * VCF_depth[0]

    s1 = np.zeros(length_of_s)
    Q = 1 / np.sqrt(2)
    for n in range(length_of_s):
        a, b = HPF(fs, vcf[n], Q)
        for m in range(0, 3):
            if n - m >= 0:
                s1[n] += b[m] * s0[n - m]

        for m in range(1, 3):
            if n - m >= 0:
                s1[n] += -a[m] * s1[n - m]

    VCA_A = [0]
    VCA_D = [0.1]
    VCA_S = [0]
    VCA_R = [0.1]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], gate, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        s1[n] *= vca[n]

    gain = velocity / 127 / np.max(np.abs(s1))
    s1 *= gain

    return s1

(5)

def percussion(fs, note_number, velocity, gate):
    if note_number == 36:
        s = bass_drum(fs, velocity, gate)
    elif note_number == 40:
        s = snare_drum(fs, velocity, gate)
    elif note_number == 42:
        s = hihat_cymbal_close(fs, velocity, gate)

    return s

(6)

fs = 44100

(7)

note_number = 42
velocity = 100
gate = 0.1

(8)

s = percussion(fs, note_number, velocity, gate)
length_of_s = len(s)

(9)

for n in range(length_of_s):
    s[n] = (s[n] + 1.0) / 2.0 * 65536.0
    if s[n] > 65535.0:
        s[n] = 65535.0
    elif s[n] < 0.0:
        s[n] = 0.0;
    s[n] = (s[n] + 0.5) - 32768

wavfile.write('hihat.wav', fs, s.astype(np.int16))

(10)

Audio('hihat.wav')

8．減算合成2（バスドラム）

　「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．サイン波と白色雑音からバスドラムの音をつくることができることを確認しなさい．

(1)

%matplotlib inline
import numpy as np
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def LPF(fs, fc, Q):
    fc /= fs
    fc = np.tan(np.pi * fc) / (2.0 * np.pi)
    a = np.zeros(3)
    b = np.zeros(3)
    a[0] = 1.0 + 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc
    a[1] = (8.0 * np.pi * np.pi * fc * fc - 2.0) / a[0]
    a[2] = (1.0 - 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc) / a[0]
    b[0] = 4.0 * np.pi * np.pi * fc * fc / a[0]
    b[1] = 8.0 * np.pi * np.pi * fc * fc / a[0]
    b[2] = 4.0 * np.pi * np.pi * fc * fc / a[0]
    a[0] = 1.0
    return a, b

(3)

def ADSR(fs, A, D, S, R, gate, duration):
    A = int(fs * A)
    D = int(fs * D)
    R = int(fs * R)
    gate = int(fs * gate)
    duration = int(fs * duration)
    e = np.zeros(duration)
    if A != 0:
        for n in range(A):
            e[n] = 1.0 - np.exp(-5.0 * n / A)

    if D != 0:
        for n in range(A, gate):
            e[n] = S + (1.0 - S) * np.exp(-5.0 * (n - A) / D)

    else:
        for n in range(A, gate):
            e[n] = S

    if R != 0:
        for n in range(gate, duration):
            e[n]= e[gate - 1] * np.exp(-5.0 * (n - gate + 1) / R)

    return e

(4)

def compressor(fs, x):
    length_of_s = len(x)

    gain = 1.0 / np.max(np.abs(x))
    x *= gain

    threshold = 0.6
    ratio = 1 / 8
    width = 0.2
    gain = 1 / (threshold + (1.0 - threshold) * ratio)

    for n in range(length_of_s):
        if x[n] < 0:
            sign_of_s = -1
        else:
            sign_of_s = 1

        abs_of_s = np.abs(x[n])

        if abs_of_s >= threshold - width / 2 and abs_of_s < threshold + width / 2:
            abs_of_s = abs_of_s + (ratio - 1) * (abs_of_s - threshold + width / 2)*(abs_of_s - threshold + width / 2) / (width * 2)
        elif abs_of_s >= threshold + width / 2:
            abs_of_s = threshold + (abs_of_s - threshold) * ratio

        x[n] = sign_of_s * abs_of_s * gain

    return x

(5)

def bass_drum(fs, velocity, gate):
    duration = 1

    length_of_s = int(fs * duration)
    sa0 = np.zeros(length_of_s)
    sb0 = np.zeros(length_of_s)

    VCO_A = [0]
    VCO_D = [0.1]
    VCO_S = [0]
    VCO_R = [0.1]
    VCO_offset = [55]
    VCO_depth = [50]

    vco = ADSR(fs, VCO_A[0], VCO_D[0], VCO_S[0], VCO_R[0], gate, duration)
    for n in range(length_of_s):
        vco[n] = VCO_offset[0] + vco[n] * VCO_depth[0]

    x = 0
    for n in range(length_of_s):
        sa0[n] = np.sin(2.0 * np.pi * x)
        delta = vco[n] / fs
        x += delta
        if x >= 1:
            x -= 1

    np.random.seed(0)
    for n in range(length_of_s):
        sb0[n] = (np.random.rand() * 2.0) - 1.0

    VCA_A = [0]
    VCA_D = [0.2]
    VCA_S = [0]
    VCA_R = [0.2]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], gate, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        sb0[n] *= vca[n]

    s0 = np.zeros(length_of_s)
    for n in range(length_of_s):
        s0[n] = sa0[n] * 0.6 + sb0[n] * 0.4

    VCF_A = [0]
    VCF_D = [0.04]
    VCF_S = [0]
    VCF_R = [0.04]
    VCF_offset = [200]
    VCF_depth = [4000]

    vcf = ADSR(fs, VCF_A[0], VCF_D[0], VCF_S[0], VCF_R[0], gate, duration)
    for n in range(length_of_s):
        vcf[n] = VCF_offset[0] + vcf[n] * VCF_depth[0]

    s1 = np.zeros(length_of_s)
    Q = 1 / np.sqrt(2)
    for n in range(length_of_s):
        a, b = LPF(fs, vcf[n], Q)
        for m in range(0, 3):
            if n - m >= 0:
                s1[n] += b[m] * s0[n - m]

        for m in range(1, 3):
            if n - m >= 0:
                s1[n] += -a[m] * s1[n - m]

    VCA_A = [0]
    VCA_D = [0.3]
    VCA_S = [0]
    VCA_R = [0.3]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], gate, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        s1[n] *= vca[n]

    s1 = compressor(fs, s1)

    gain = velocity / 127 / np.max(np.abs(s1))
    s1 *= gain

    return s1

(6)

def percussion(fs, note_number, velocity, gate):
    if note_number == 36:
        s = bass_drum(fs, velocity, gate)
    elif note_number == 40:
        s = snare_drum(fs, velocity, gate)
    elif note_number == 42:
        s = hihat_cymbal_close(fs, velocity, gate)

    return s

(7)

fs = 44100

(8)

note_number = 36
velocity = 100
gate = 0.1

(9)

s = percussion(fs, note_number, velocity, gate)
length_of_s = len(s)

(10)

for n in range(length_of_s):
    s[n] = (s[n] + 1.0) / 2.0 * 65536.0
    if s[n] > 65535.0:
        s[n] = 65535.0
    elif s[n] < 0.0:
        s[n] = 0.0;
    s[n] = (s[n] + 0.5) - 32768

wavfile.write('bass_drum.wav', fs, s.astype(np.int16))

(11)

Audio('bass_drum.wav')

9．減算合成3（スネアドラム）

　「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．矩形波と白色雑音からスネアドラムの音をつくることができることを確認しなさい．

(1)

%matplotlib inline
import numpy as np
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def LPF(fs, fc, Q):
    fc /= fs
    fc = np.tan(np.pi * fc) / (2.0 * np.pi)
    a = np.zeros(3)
    b = np.zeros(3)
    a[0] = 1.0 + 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc
    a[1] = (8.0 * np.pi * np.pi * fc * fc - 2.0) / a[0]
    a[2] = (1.0 - 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc) / a[0]
    b[0] = 4.0 * np.pi * np.pi * fc * fc / a[0]
    b[1] = 8.0 * np.pi * np.pi * fc * fc / a[0]
    b[2] = 4.0 * np.pi * np.pi * fc * fc / a[0]
    a[0] = 1.0
    return a, b

(3)

def ADSR(fs, A, D, S, R, gate, duration):
    A = int(fs * A)
    D = int(fs * D)
    R = int(fs * R)
    gate = int(fs * gate)
    duration = int(fs * duration)
    e = np.zeros(duration)
    if A != 0:
        for n in range(A):
            e[n] = 1.0 - np.exp(-5.0 * n / A)

    if D != 0:
        for n in range(A, gate):
            e[n] = S + (1.0 - S) * np.exp(-5.0 * (n - A) / D)

    else:
        for n in range(A, gate):
            e[n] = S

    if R != 0:
        for n in range(gate, duration):
            e[n]= e[gate - 1] * np.exp(-5.0 * (n - gate + 1) / R)

    return e

(4)

def compressor(fs, x):
    length_of_s = len(x)

    gain = 1.0 / np.max(np.abs(x))
    x *= gain

    threshold = 0.6
    ratio = 1 / 8
    width = 0.2
    gain = 1 / (threshold + (1.0 - threshold) * ratio)

    for n in range(length_of_s):
        if x[n] < 0:
            sign_of_s = -1
        else:
            sign_of_s = 1

        abs_of_s = np.abs(x[n])

        if abs_of_s >= threshold - width / 2 and abs_of_s < threshold + width / 2:
            abs_of_s = abs_of_s + (ratio - 1) * (abs_of_s - threshold + width / 2)*(abs_of_s - threshold + width / 2) / (width * 2)
        elif abs_of_s >= threshold + width / 2:
            abs_of_s = threshold + (abs_of_s - threshold) * ratio

        x[n] = sign_of_s * abs_of_s * gain

    return x

(5)

def snare_drum(fs, velocity, gate):
    duration = 1

    length_of_s = int(fs * duration)
    sa0 = np.zeros(length_of_s)
    sb0 = np.zeros(length_of_s)

    VCO_A = [0]
    VCO_D = [0]
    VCO_S = [1]
    VCO_R = [0]
    VCO_offset = [150]
    VCO_depth = [0]

    vco = ADSR(fs, VCO_A[0], VCO_D[0], VCO_S[0], VCO_R[0], duration, duration)
    for n in range(length_of_s):
        vco[n] = VCO_offset[0] + vco[n] * VCO_depth[0]

    x = 0
    for n in range(length_of_s):
        if x < 0.5:
            sa0[n] = 1
        else:
            sa0[n] = -1

        delta = vco[n] / fs

        if 1 - delta <= x and x < 1:
            t = (x - 1) / delta
            d = t * t + 2 * t + 1
            sa0[n] += d
        elif 0 <= x and x < delta:
            t = x / delta
            d = -t * t + 2 * t - 1
            sa0[n] += d

        if 0.5 - delta <= x and x < 0.5:
            t = (x - 0.5) / delta
            d = t * t + 2 * t + 1
            sa0[n] -= d
        elif 0.5 <= x and x < 0.5 + delta:
            t = (x - 0.5) / delta
            d = -t * t + 2 * t - 1
            sa0[n] -= d

        x += delta
        if x >= 1:
            x -= 1

    np.random.seed(0)
    for n in range(length_of_s):
        sb0[n] = (np.random.rand() * 2.0) - 1.0

    VCA_A = [0]
    VCA_D = [0]
    VCA_S = [1]
    VCA_R = [0]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], duration, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        sb0[n] *= vca[n]

    s0 = np.zeros(length_of_s)
    for n in range(length_of_s):
        s0[n] = sa0[n] * 0.3 + sb0[n] * 0.7

    VCF_A = [0]
    VCF_D = [0.1]
    VCF_S = [0]
    VCF_R = [0.1]
    VCF_offset = [8000]
    VCF_depth = [-7800]

    vcf = ADSR(fs, VCF_A[0], VCF_D[0], VCF_S[0], VCF_R[0], gate, duration)
    for n in range(length_of_s):
        vcf[n] = VCF_offset[0] + vcf[n] * VCF_depth[0]

    s1 = np.zeros(length_of_s)
    Q = 1 / np.sqrt(2)
    for n in range(length_of_s):
        a, b = LPF(fs, vcf[n], Q)
        for m in range(0, 3):
            if n - m >= 0:
                s1[n] += b[m] * s0[n - m]

        for m in range(1, 3):
            if n - m >= 0:
                s1[n] += -a[m] * s1[n - m]

    VCA_A = [0]
    VCA_D = [0.2]
    VCA_S = [0]
    VCA_R = [0.2]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], gate, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        s1[n] *= vca[n]

    s1 = compressor(fs, s1)

    gain = velocity / 127 / np.max(np.abs(s1))
    s1 *= gain

    return s1

(6)

def percussion(fs, note_number, velocity, gate):
    if note_number == 36:
        s = bass_drum(fs, velocity, gate)
    elif note_number == 40:
        s = snare_drum(fs, velocity, gate)
    elif note_number == 42:
        s = hihat_cymbal_close(fs, velocity, gate)

    return s

(7)

fs = 44100

(8)

note_number = 40
velocity = 100
gate = 0.1

(9)

s = percussion(fs, note_number, velocity, gate)
length_of_s = len(s)

(10)

for n in range(length_of_s):
    s[n] = (s[n] + 1.0) / 2.0 * 65536.0
    if s[n] > 65535.0:
        s[n] = 65535.0
    elif s[n] < 0.0:
        s[n] = 0.0;
    s[n] = (s[n] + 0.5) - 32768

wavfile.write('snare_drum.wav', fs, s.astype(np.int16))

(11)

Audio('snare_drum.wav')

10．自動演奏

　「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．ドラムセットの演奏が聞こえてくることを確かめてください．

(1)

%matplotlib inline
import numpy as np
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def LPF(fs, fc, Q):
    fc /= fs
    fc = np.tan(np.pi * fc) / (2.0 * np.pi)
    a = np.zeros(3)
    b = np.zeros(3)
    a[0] = 1.0 + 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc
    a[1] = (8.0 * np.pi * np.pi * fc * fc - 2.0) / a[0]
    a[2] = (1.0 - 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc) / a[0]
    b[0] = 4.0 * np.pi * np.pi * fc * fc / a[0]
    b[1] = 8.0 * np.pi * np.pi * fc * fc / a[0]
    b[2] = 4.0 * np.pi * np.pi * fc * fc / a[0]
    a[0] = 1.0
    return a, b

(3)

def HPF(fs, fc, Q):
    fc /= fs
    fc = np.tan(np.pi * fc) / (2.0 * np.pi)
    a = np.zeros(3)
    b = np.zeros(3)
    a[0] = 1.0 + 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc;
    a[1] = (8.0 * np.pi * np.pi * fc * fc - 2.0) / a[0]
    a[2] = (1.0 - 2.0 * np.pi * fc / Q + 4.0 * np.pi * np.pi * fc * fc) / a[0]
    b[0] = 1.0 / a[0]
    b[1] = -2.0 / a[0]
    b[2] = 1.0 / a[0]
    a[0] = 1.0
    return a, b

(4)

def ADSR(fs, A, D, S, R, gate, duration):
    A = int(fs * A)
    D = int(fs * D)
    R = int(fs * R)
    gate = int(fs * gate)
    duration = int(fs * duration)
    e = np.zeros(duration)
    if A != 0:
        for n in range(A):
            e[n] = 1.0 - np.exp(-5.0 * n / A)

    if D != 0:
        for n in range(A, gate):
            e[n] = S + (1.0 - S) * np.exp(-5.0 * (n - A) / D)

    else:
        for n in range(A, gate):
            e[n] = S

    if R != 0:
        for n in range(gate, duration):
            e[n]= e[gate - 1] * np.exp(-5.0 * (n - gate + 1) / R)

    return e

(5)

def compressor(fs, x):
    length_of_s = len(x)

    gain = 1.0 / np.max(np.abs(x))
    x *= gain

    threshold = 0.6
    ratio = 1 / 8
    width = 0.2
    gain = 1 / (threshold + (1.0 - threshold) * ratio)

    for n in range(length_of_s):
        if x[n] < 0:
            sign_of_s = -1
        else:
            sign_of_s = 1

        abs_of_s = np.abs(x[n])

        if abs_of_s >= threshold - width / 2 and abs_of_s < threshold + width / 2:
            abs_of_s = abs_of_s + (ratio - 1) * (abs_of_s - threshold + width / 2)*(abs_of_s - threshold + width / 2) / (width * 2)
        elif abs_of_s >= threshold + width / 2:
            abs_of_s = threshold + (abs_of_s - threshold) * ratio

        x[n] = sign_of_s * abs_of_s * gain

    return x

(6)

def hihat_cymbal_close(fs, velocity, gate):
    duration = 1

    length_of_s = int(fs * duration)
    s0 = np.zeros(length_of_s)

    np.random.seed(0)
    for n in range(length_of_s):
        s0[n] = (np.random.rand() * 2.0) - 1.0

    VCF_A = [0]
    VCF_D = [0]
    VCF_S = [1]
    VCF_R = [0]
    VCF_offset = [10000]
    VCF_depth=[0]

    vcf = ADSR(fs, VCF_A[0], VCF_D[0], VCF_S[0], VCF_R[0], gate, duration)
    for n in range(length_of_s):
        vcf[n] = VCF_offset[0] + vcf[n] * VCF_depth[0]

    s1 = np.zeros(length_of_s)
    Q = 1 / np.sqrt(2)
    for n in range(length_of_s):
        a, b = HPF(fs, vcf[n], Q)
        for m in range(0, 3):
            if n - m >= 0:
                s1[n] += b[m] * s0[n - m]

        for m in range(1, 3):
            if n - m >= 0:
                s1[n] += -a[m] * s1[n - m]

    VCA_A = [0]
    VCA_D = [0.1]
    VCA_S = [0]
    VCA_R = [0.1]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], gate, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        s1[n] *= vca[n]

    gain = velocity / 127 / np.max(np.abs(s1))
    s1 *= gain

    return s1

(7)

def bass_drum(fs, velocity, gate):
    duration = 1

    length_of_s = int(fs * duration)
    sa0 = np.zeros(length_of_s)
    sb0 = np.zeros(length_of_s)

    VCO_A = [0]
    VCO_D = [0.1]
    VCO_S = [0]
    VCO_R = [0.1]
    VCO_offset = [55]
    VCO_depth = [50]

    vco = ADSR(fs, VCO_A[0], VCO_D[0], VCO_S[0], VCO_R[0], gate, duration)
    for n in range(length_of_s):
        vco[n] = VCO_offset[0] + vco[n] * VCO_depth[0]

    x = 0
    for n in range(length_of_s):
        sa0[n] = np.sin(2.0 * np.pi * x)
        delta = vco[n] / fs
        x += delta
        if x >= 1:
            x -= 1

    np.random.seed(0)
    for n in range(length_of_s):
        sb0[n] = (np.random.rand() * 2.0) - 1.0

    VCA_A = [0]
    VCA_D = [0.2]
    VCA_S = [0]
    VCA_R = [0.2]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], gate, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        sb0[n] *= vca[n]

    s0 = np.zeros(length_of_s)
    for n in range(length_of_s):
        s0[n] = sa0[n] * 0.6 + sb0[n] * 0.4

    VCF_A = [0]
    VCF_D = [0.04]
    VCF_S = [0]
    VCF_R = [0.04]
    VCF_offset = [200]
    VCF_depth = [4000]

    vcf = ADSR(fs, VCF_A[0], VCF_D[0], VCF_S[0], VCF_R[0], gate, duration)
    for n in range(length_of_s):
        vcf[n] = VCF_offset[0] + vcf[n] * VCF_depth[0]

    s1 = np.zeros(length_of_s)
    Q = 1 / np.sqrt(2)
    for n in range(length_of_s):
        a, b = LPF(fs, vcf[n], Q)
        for m in range(0, 3):
            if n - m >= 0:
                s1[n] += b[m] * s0[n - m]

        for m in range(1, 3):
            if n - m >= 0:
                s1[n] += -a[m] * s1[n - m]

    VCA_A = [0]
    VCA_D = [0.3]
    VCA_S = [0]
    VCA_R = [0.3]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], gate, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        s1[n] *= vca[n]

    s1 = compressor(fs, s1)

    gain = velocity / 127 / np.max(np.abs(s1))
    s1 *= gain

    return s1

(8)

def snare_drum(fs, velocity, gate):
    duration = 1

    length_of_s = int(fs * duration)
    sa0 = np.zeros(length_of_s)
    sb0 = np.zeros(length_of_s)

    VCO_A = [0]
    VCO_D = [0]
    VCO_S = [1]
    VCO_R = [0]
    VCO_offset = [150]
    VCO_depth = [0]

    vco = ADSR(fs, VCO_A[0], VCO_D[0], VCO_S[0], VCO_R[0], duration, duration)
    for n in range(length_of_s):
        vco[n] = VCO_offset[0] + vco[n] * VCO_depth[0]

    x = 0
    for n in range(length_of_s):
        if x < 0.5:
            sa0[n] = 1
        else:
            sa0[n] = -1

        delta = vco[n] / fs

        if 1 - delta <= x and x < 1:
            t = (x - 1) / delta
            d = t * t + 2 * t + 1
            sa0[n] += d
        elif 0 <= x and x < delta:
            t = x / delta
            d = -t * t + 2 * t - 1
            sa0[n] += d

        if 0.5 - delta <= x and x < 0.5:
            t = (x - 0.5) / delta
            d = t * t + 2 * t + 1
            sa0[n] -= d
        elif 0.5 <= x and x < 0.5 + delta:
            t = (x - 0.5) / delta
            d = -t * t + 2 * t - 1
            sa0[n] -= d

        x += delta
        if x >= 1:
            x -= 1

    np.random.seed(0)
    for n in range(length_of_s):
        sb0[n] = (np.random.rand() * 2.0) - 1.0

    VCA_A = [0]
    VCA_D = [0]
    VCA_S = [1]
    VCA_R = [0]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], duration, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        sb0[n] *= vca[n]

    s0 = np.zeros(length_of_s)
    for n in range(length_of_s):
        s0[n] = sa0[n] * 0.3 + sb0[n] * 0.7

    VCF_A = [0]
    VCF_D = [0.1]
    VCF_S = [0]
    VCF_R = [0.1]
    VCF_offset = [8000]
    VCF_depth = [-7800]

    vcf = ADSR(fs, VCF_A[0], VCF_D[0], VCF_S[0], VCF_R[0], gate, duration)
    for n in range(length_of_s):
        vcf[n] = VCF_offset[0] + vcf[n] * VCF_depth[0]

    s1 = np.zeros(length_of_s)
    Q = 1 / np.sqrt(2)
    for n in range(length_of_s):
        a, b = LPF(fs, vcf[n], Q)
        for m in range(0, 3):
            if n - m >= 0:
                s1[n] += b[m] * s0[n - m]

        for m in range(1, 3):
            if n - m >= 0:
                s1[n] += -a[m] * s1[n - m]

    VCA_A = [0]
    VCA_D = [0.2]
    VCA_S = [0]
    VCA_R = [0.2]
    VCA_offset = [0]
    VCA_depth = [1]

    vca = ADSR(fs, VCA_A[0], VCA_D[0], VCA_S[0], VCA_R[0], gate, duration)
    for n in range(length_of_s):
        vca[n] = VCA_offset[0] + vca[n] * VCA_depth[0]

    for n in range(length_of_s):
        s1[n] *= vca[n]

    s1 = compressor(fs, s1)

    gain = velocity / 127 / np.max(np.abs(s1))
    s1 *= gain

    return s1

(9)

def percussion(fs, note_number, velocity, gate):
    if note_number == 36:
        s = bass_drum(fs, velocity, gate)
    elif note_number == 40:
        s = snare_drum(fs, velocity, gate)
    elif note_number == 42:
        s = hihat_cymbal_close(fs, velocity, gate)

    return s

(10)

score = np.array([[1, 1920, 42, 100, 96],
                  [1, 2160, 42, 100, 96],
                  [1, 2400, 42, 100, 96],
                  [1, 2640, 42, 100, 96],
                  [1, 2880, 42, 100, 96],
                  [1, 3120, 42, 100, 96],
                  [1, 3360, 42, 100, 96],
                  [1, 3600, 42, 100, 96],
                  [1, 3840, 42, 100, 96],
                  [1, 4080, 42, 100, 96],
                  [1, 4320, 42, 100, 96],
                  [1, 4560, 42, 100, 96],
                  [1, 4800, 42, 100, 96],
                  [1, 5040, 42, 100, 96],
                  [1, 5280, 42, 100, 96],
                  [1, 5520, 42, 100, 96],
                  [2, 1920, 36, 100, 96],
                  [2, 2880, 36, 100, 96],
                  [2, 3840, 36, 100, 96],
                  [2, 4800, 36, 100, 96],
                  [3, 2400, 40, 100, 96],
                  [3, 3360, 40, 100, 96],
                  [3, 4320, 40, 100, 96],
                  [3, 5280, 40, 100, 96]])

division = 480
tempo = 120
number_of_track = 3
end_of_track = 6
number_of_note = score.shape[0]

(11)

fs = 44100
length_of_s = int(fs * (end_of_track + 2))
track = np.zeros((length_of_s, number_of_track))
s = np.zeros(length_of_s)

(12)

for i in range(number_of_note):
    j = int(score[i, 0] - 1)
    onset = (score[i, 1] / division) * (60 / tempo)
    note_number = score[i, 2]
    velocity = score[i, 3]
    gate = (score[i, 4] / division) * (60 / tempo)
    x = percussion(fs, note_number, velocity, gate)
    offset = int(fs * onset)
    length_of_x = len(x)
    for n in range(length_of_x):
        track[offset + n, j] += x[n]

(13)

for j in range(number_of_track):
    for n in range(length_of_s):
        s[n] += track[n, j]

(14)

master_volume = 0.5
s /= np.max(np.abs(s))
s *= master_volume

(15)

for n in range(length_of_s):
    s[n] = (s[n] + 1.0) / 2.0 * 65536.0
    if s[n] > 65535.0:
        s[n] = 65535.0
    elif s[n] < 0.0:
        s[n] = 0.0;
    s[n] = (s[n] + 0.5) - 32768

wavfile.write('p6.wav', fs, s.astype(np.int16))

(16)

Audio('p6.wav')

11．レポートについて

　下記の課題について，レポートを作成しなさい．

[1] ウェブサイトなどから楽器音の音データを入手し，波形，周波数特性，スペクトログラムを表示しなさい．これらをお手本とすることで，楽器音を合成するプログラムをつくりなさい．どのような手順で音をつくったのか，自分のプログラムについて説明してください．

[2] 自分でつくった楽器音を使ってカノンの演奏を確認し、そのWAVEファイルをメールに添付して提出してください．

12．創成課題

[1] 自分でつくった楽器音を使ってカノン以外の音楽データを作成しなさい．少なくともふたつ以上の楽器音をつくり，音楽データを作成すること．MIDIファイルから楽譜データを読み取るプログラムをつくる場合は，http://floor13.sakura.ne.jp/book09/book09.htmlからprogram.zipをダウンロードし、そのなかにあるp4_2.pyを参考にしてください．

[2] 創成課題の発表会に向けて，プレゼンテーションのための資料を作成しなさい．かならず，音楽データのデモをまじえること．どのような手順で音をつくったのか，具体的な説明があるものを高く評価します．

創成課題の作品

2022年度
2023年度

参考文献

青木 直史, ``Pythonではじめる 音のプログラミング　～ コンピュータミュージックの信号処理 ～ ,'' オーム社, 2022.

青木 直史, ``Pythonではじめる 音のプログラミング　～ コンピュータミュージックの信号処理 ～ ,'' Googleブックス, 2022.

Last Modified: November 5 12:00 JST 2022 by Naofumi Aoki
E-mail: aoki@ime.ist.hokudai.ac.jp