何かあれば GitHub のリポジトリに issue を作るか ryukau@gmail.com までお気軽にどうぞ。


インデックスに戻る

Update: 2019-01-09

Table of Contents

Inviscid な Burgers 方程式

Landajuela さんの BURGERS EUQATION で紹介されていた Godunov’s method で変形した Burgers 方程式を実装して音を入力します。

コードについて

次のリンクからコードをまとめて読むことができます。

リンク先の propagation.pySoX を使っています。

今回使った Python3 のライブラリです。

シミュレーション

ここでシミュレーションするのは inviscid な Burgers 方程式です。

\[ u_t + u u_x = 0 \]

\(u_t, u_x\)偏微分の表記です。

Landajuela さんの BURGERS EUQATION で紹介されていた Godunov’s method で変形した Burgers 方程式です。

\[ U_j^{n+1} = U_j^n - \frac{k}{h} \left( F(U_j^n, U_{j+1}^n) - F(U_{j-1}^n, U_j^n) \right) \]

\(F(U, V)\) は次のように定義されています。

\[ \begin{aligned} F(U, V) &= \frac{(u^*)^2}{2}\\ \text{If}\;U \geq V\;\text{then}\qquad u^* &= \begin{cases} U, &\text{if}\;\frac{U+V}{2} > 0,\\ V, &\text{in other case.} \end{cases}\\ \text{If}\;U < V\;\text{then}\qquad u^* &= \begin{cases} U, &\text{if}\;U > 0,\\ V, &\text{if}\;V > 0,\\ 0, &\text{if}\;U \leq 0 \leq V. \end{cases} \end{aligned} \]

実装します。

import matplotlib.pyplot as pyplot
import numpy

class Burgers1D:
    """
    dx = 1, dt = 1 で固定した inviscid Burgers' equation 。
    """

    def __init__(self, length):
        """length はシミュレーションする波の配列の長さ。"""
        self.wave = numpy.zeros((2, length))
        self.reset()

    def u_star(self, u, v):
        return numpy.where(
            u >= v,
            numpy.where((u + v) / 2 > 0, u, v),
            numpy.where(u > 0, u, numpy.where(v < 0, v, 0)),
        )

    def step(self):
        self.wave = numpy.roll(self.wave, 1, axis=0)

        last = self.wave.shape[1] - 1

        wave_l = numpy.roll(self.wave[1], 1)
        wave_r = numpy.roll(self.wave[1], -1)
        u_star_l = self.u_star(self.wave[1], wave_r)
        u_star_r = self.u_star(wave_l, self.wave[1])
        self.wave[0] = self.wave[1] - (u_star_l * u_star_l - u_star_r * u_star_r) / 2

        self.wave[0][0] = 0
        self.wave[0][last] = 0

        if self.pick_y != 0:
            self.wave[0][self.pick_x] = self.pick_y

    def reset(self):
        self.wave.fill(0)

    def pick(self, x, y):
        """
        x の範囲は [0, 1] 。
        y の範囲は [-1, 1] 。 |y| > 1 で発散する。
        """
        x = max(0, min(x, 1))
        self.pick_x = numpy.int32(x * self.wave.shape[1])
        self.pick_y = y

    def state(self):
        return self.wave[0]

length = 512
burgers = Burgers1D(length)

signal = numpy.sin(numpy.linspace(0, 8 * numpy.pi, length))

for i, pick_y in enumerate(signal):
    burgers.pick(0, pick_y)
    burgers.step()

pyplot.grid()
pyplot.plot(burgers.state())
pyplot.show()

コードを実行すると次のようなプロットが出力されます。

Image of a result of simulation of inviscid Burgers' equation.

デモ

キャンバスをクリックすると波が起こります。

中央の黒い線は \(U = 0\) を表しています。上半分は \(U > 0\) 、下半分は \(U < 0\) となります。

今回のシミュレーションでは入力の符号によって波が進む方向が決まるようです。

音を入力する

波の進む方向が一方向に固定されるので、入力信号がすべて正の値になるように整形します。またシミュレーションが発散しないように入力信号の最大値が 1 以下となるようにします。

次の図は dc = 0.6amp = 0.8 として整形した入力信号の例です。

Image of input signal.

上の図の入力信号から得られた inviscid な Burgers 方程式のシミュレーション結果です。

Image of raw output signal.

シミュレーションから得られた信号を wav ファイルとして保存するときに直流を切りたいのですが、このままハイパスフィルタをかけると音の始まりに大きなポップノイズが乗ってしまいます。今回は信号の numpy.median の値をしきい値として、信号の開始からしきい値を超えるまでの区間を numpy.median の値に置き換えることでポップノイズを抑えました。

実装は次のようになります。

import numpy
import scipy.signal
import soundfile
from burgers import Burgers1D

samplerate = 44100
duration = 0.4
frequency = 60
amp = 0.8
dc = 0.6

phase = numpy.linspace(
    0,
    2 * numpy.pi * frequency * duration,
    int(samplerate * duration),
)
input_signal = numpy.sin(phase)

# 入力信号を整形。
amp = dc * amp if dc < 0.5 else (1 - dc) * amp
signal = (dc - amp) + amp * (input_signal + 1)

input_signal = numpy.copy(signal)

# シミュレーション。
burgers = Burgers1D(257)
read_index = burgers.wave.shape[1] - 2
for i, y in enumerate(signal):
    burgers.pick(0, y)
    burgers.step()
    signal[i] = burgers.state()[read_index]

raw_signal = numpy.copy(signal)

# ポップノイズを防ぐために信号の開始からしきい値を超えるまでの値を置き換え。
threshold = numpy.median(signal)
i = 0
while signal[i] < threshold:
    i += 1
signal[:i] = threshold
signal -= threshold

# 直流の除去。
hp_sos = scipy.signal.butter(4, 2 * 20 / samplerate, btype="highpass", output="sos")
signal = scipy.signal.sosfilt(hp_sos, signal)

soundfile.write("output.wav", signal, samplerate)

直流を除去した出力信号です。

Image of processed output signal.

音のサンプル

ここまでのコードを propagation.py にまとめてコマンドラインから使えるようにしました。声のサンプルのラベルに書いている -d-lpropagation.py のオプションです。 -d は DC オフセット、 -lBurgers1D のインスタンスを生成するときに渡す length です。

1000Hz のサイン波の入力信号と Burgers1D に通した結果です。

ブラウンノイズの入力信号と Burgers1D に通した結果です。

サイン波とブラウンノイズは SoX で生成しました。

sox -n sin.wav synth 1.0 sine 1000.0
sox -n noise.wav synth 0.5 brownnoise

Freesound で見つけた AlienAudio さんによる人の声のサンプルBurgers1D に通しました。