目次

1.はじめに
 2.問題の原因調査
 3.位相について
 4.位相を理解するためのポイント
 5.アプローチ
 6.実装
 7.検証
 8.まとめ
 9.今後の検討
 10.所感
 11.参考

はじめに

Python+PyAudioを活用し、正弦波（sin波）を生成して音鳴らす場合、音を切り替える際に「プチッ」というノイズが聞こえることがあり、不快感を生じさせる恐れがあります。この記事では、そのようなノイズの発生を抑える方法について紹介します。

より具体的に理解して頂くために、デモ動画を紹介します。以下のデモでは、音が「シ→ド→シ」ととそれぞれ2秒ずつ切り替わります。

「ノイズあり」デモ

「ノイズ対処済み」デモ

「ノイズあり」のデモでは、音の切り替わり時に「プチッ」としたノイズが発生します。一方、「ノイズ対処済み」のデモでは、音と音のつなぎ目がより滑らかに聞こえるかと思います。

この記事では、音の切り替わりをより滑らかにする方法を解説します。具体的な実装方法に興味のある方は、「実装」セクションをご覧ください。

問題の原因調査

問題が生じていそうな部分

「ノイズあり」の音を確認すると、音の切り替わり部分で「プチッ」というノイズが発生しているように聞こえます。そこで、音が切り替わる部分（シ→ドの切り替わり）について、波の波形を確認することで、何が起きているのかを調査したいと思います。

なお、今回は音を切り替えたときのみに焦点を当て、音の再生開始時や終了時のノイズには対処していません。将来的にはこれらの状況にも対処できるようにしたいと考えています。

「ノイズあり」のコードについて

音の切り替わりを確認する前に、まずは「ノイズあり」の実装について記載します。このコードを出発点として、問題の原因特定をしたいと思います。

なお、本記事は先日公開した「PythonとPyAudioを使ってsin波を作成して音を鳴らす方法」のコードを拡張して解説します。基本的なコードの説明・動作確認方法はそちらに記載していますので、適宜ご参照ください。

「ノイズあり」のコード

import numpy as np
import pyaudio

def sound_sin(A, f, t):
    return A * np.sin(2 * np.pi * f * t)

A = 0.5
f = [493.88, 523.23]
play_length = 2
fs = 44100

p = pyaudio.PyAudio()

stream = p.open(format=pyaudio.paFloat32,
                channels=1,
                rate=fs,
                output=True)

samples = np.arange(fs * play_length) / fs

y1 = sound_sin(A, f[0], samples)
stream.write(y1.astype(np.float32).tobytes())

y2 = sound_sin(A, f[1], samples)
stream.write(y2.astype(np.float32).tobytes())

stream.write(y1.astype(np.float32).tobytes())

stream.stop_stream()
stream.close()
p.terminate()

コードの説明

上記のコードを簡単に説明します。

パラメータ：

振幅（A）: 音の大きさを決定します。値が大きいほど、出力される音は大きくなります
周波数（f）: 音の高さを決定します。このコードでは、周波数をリスト化しています
音の長さ（play_length）: 出力する音の秒数を指定します
サンプリングレート（fs）: サンプリングレートです。サンプリングレートは、1秒間を何個の点に分割するかを決めます
時刻配列（samples）:音を生成するための時刻を表します

処理の流れ:

sound_sin関数は、振幅A、周波数f、および時刻配列samplesを引数として受け取り、sin波の配列を返します。これにより、音波の形状が決定されます
生成された音波データはPyAudioを通じて音として出力されます
for文の中で周波数fを変えることで、異なる音を生成することができます

音の切り替わり部分の確認

音の切り替わり部分でノイズが発生しているので、y1とy2の切り替わり部分を確認します。

以下は、「y1からy2へ切り替わったとき」と、「y2からy1へ切り替わったとき」について図示した結果になります。y1とy2のつながりをわかりやすく可視化するため、あえて二つの波を一つのグラフに乗せ、切り替わりの前後それぞれ0.004秒ほどの区間を、拡大して図示しています。y1は青色で、y2は赤色です。横軸は時刻samplesの値（秒）であり、縦軸はy1, y2の値です。

図1を見ると、y1やy2単体では波が滑らかに連続しているのに対し、音の切り替わり部分である2.00秒付近は、波が不連続となっていることがわかります。また、図2を見ると、同様に切り替わり部分（4.00秒付近）で波が不連続となっていることがわかります。

図1: y1からy2へ切り替わったとき

図2: y2からy1へ切り替わったとき

ノイズが生じる原因の仮説

波が連続している部分（「シ」「ド」それぞれの音が鳴っている部分）ではノイズが聞こえないのに対し、切り替わり部分だけノイズが発生していることから、切り替わり部分で滑らかにつながってないことによって、「プチッ」としたノイズが発生していることが考えられます。

そこで本記事では、この音の切り替わり部分を滑らかにつなげることで、このノイズが低減されるかを確認します。そのための手段として、初期位相を活用します。初期位相により、以前の波の終わり位置情報を初期位相として保持し、波の連続性を保つことができるようになります。

位相について

位相や初期位相の具体的な説明に入る前に、まずは正弦波、位相、初期位相の式について記載します。これらの式は、2つの波に連続性を持たせるために利用できます。詳しい使い方は、後述するアプローチ章で説明しますので、まずは前提としてご理解ください。なお、以降は独学で理解した内容が続くので、間違いがあればご指摘ください。

正弦波（sin波）の式（参考：Wikipedia）

正弦波の式は以下です。ここではこちらを参照して、角周波数ωを2πfに式変形しています。

\displaylines{
y(t) = A\sin(\omega t + \varphi) = A\sin(2\pi f t + \varphi)
}

各パラメータは、

振幅（A）：音の大きさを決定します。値が大きいほど音は大きくなります
角周波数（ω）：一秒間に進む角度を表します。単位はrad（ラジアン※）/sです。
周波数（f）：音の高さを決定します。値が大きいほど音は高くなります
時刻（t）：音を生成する時刻を表します
初期位相（φ）：波が、周期の中のどの位置にいるかを表しています。角度であり、単位はradです

※ラジアンとは、単位円（半径１の円）で，長さが１の弧に対する中心角のことを指します。例えば360°は2π[rad], 90°はπ/2[rad]になります。詳しい説明は、様々なサイトで行われているので、今回は説明は省略します。こちらのサイトなどは大変参考になりました。

位相の式（参考：高校生から味わう理論物理入門）

位相の式は以下になります。位相は角度であり、波yの大きさを決定づけます。なお、各パラメータは正弦波の式で説明したものと同一です。

\theta(t) = 2\pi f t + \varphi

初期位相の式

初期位相は、時刻t=0の時の位相です。ゆえに以下で計算されます、

\theta(0) = \varphi

正弦波と位相を組み合わせる

正弦波の式は、位相を用いることでよりわかりやすく書き変えることができます。

\displaylines{
y(t) = A\sin(θ(t))\\
θ(t) = 2\pi f t + \varphi
}

ここでAは定数であり、sin関数はθ(t)の値によって定まるので、y(t)の値は位相θ(t)によって決まることが分かります。

位相を理解するためのポイント

位相を理解するためのポイントは、大きく分けて2つあります。

ポイント1: （初期）位相とは角度である

（初期）位相について理解すべき最初の点は、（初期）位相とは「角度」（単位はラジアン※）であるということです。

ポイント2:　ある波の最後の位相は、連続する次の波の初期位相となる

ある波の最後の位相は、連続する次の波初期位相になります。例えば時刻2.0秒におけるy1の位相が、y2の初期位相になります。

アプローチ

2つの波を滑らかにつなぐことを考えると、波y1の末端の点（=終点）におけるy1の値と、y2の始点におけるy2の値が等しくなっていればつながると考えられます。例えば、冒頭で説明した図1が、図3のようになれば、滑らかに聞こえる可能性があります。

（再掲）図1: y1からy2へ切り替わったとき

図3. y1からy2へ切り替わったとき

ここで、y1, y2の値は、前述の式の通り、位相θによって決まります。すなわち、y1の終点における位相と、y2の初期位相が一致していれば滑らかに繋がりそうです。
（再掲）

\displaylines{
y(t) = A\sin(θ(t))\\
θ(t) = 2\pi f t + \varphi
}

そこで、ここからは具体例を用いて、y1の終点における位相と、y2初期位相を3つのステップで一致させます。

設定

まず、具体例で使用する設定について説明します。

波y1の設定
- 初期位相は0[rad]とする
- 時刻0秒から、2秒間の波である
- 周波数は493.88Hzとする
- 振幅は0.5とする
波y2の設定
- 初期位相は波y1の最後の位相とする
- 時刻2秒から、2秒間の波である
- 周波数は523.23Hzとする
- 振幅は0.5とする

ここで大事な点として、y1の初期位相は角度（[rad]）である必要があります。今回であれば0[rad]としています。

ステップ1. 波y1の、終点における位相を計算する

まず波y1の時刻tにおける値は、正弦波の式によって以下の通りでした。

\displaylines{
y(t) = A\sin(θ(t))\\
θ(t) = 2\pi f t + \varphi
}

そこで、まずはy1の終点、すなわちt=2における位相を計算します。具体的には波y1に関する設定を、上記の式に代入します。

初期位相はゼロのため、φ=0
時刻tは、y1の終点の時刻 t=2.0(s)
周波数f=493.88Hz

すると、y1の終点（t=2秒）での位相$\theta_1(2.0)$は次のように計算できます：

\theta_1(2.0)[rad] = 2\pi[rad] * 493.88[Hz] * 2.0[s] + 0[rad]

結果として、t=2.0における位相θは約6203[rad]になります。なお、[]内は単位を表します。

以上より、y1の終点（時刻t=2.0）における位相を求めることができました。

ステップ2. 波y1の終点における位相と、波y2の初期位相を一致させる

次に、波y2の初期位相を、波y1の終点の位相と一致させます。まず、正弦波と位相の式は以下の通りでした。

\displaylines{
y(t) = A\sin(θ(t))\\
θ(t) = 2\pi f t + \varphi
}

この$\varphi$の部分に、ステップ1で求めた$θ_1(2.0)$の値を代入することで、y2の初期位相を設定することができます。

具体的には、正弦波の式は以下になります。

\displaylines{
y_2(t) = A\sin(2 \pi f_2 t + θ_1(2.0)) \\
}

この式に、波y2に関する設定を代入します。

Aは振幅であり0.5
$f_2$は$y_2$の周波数であり、523.23Hz
tは時刻
$θ_1(2.0)$ は波y1の時刻t=2.0における位相

以上のように、ステップ1とステップ2によって、y1の終点における位相とy2の初期位相を一致させることで、2つの異なる音波が接続された際に波形が自然につながり、ノイズ低減が期待されます。

ステップ3. θの大きさをスケール

このステップでは、φの値が大きくなるという問題に取り組みます。φの値が大きいと、数値計算上の誤差が生じやすくなり、プログラムの実行効率も低下する可能性があります。例えば、波y1の終点での位相は約6203ラジアンであり、この大きな値を扱うのは非効率的です。

解決策として、sin関数が周期関数であることを利用します。周期関数とは、Wikipediaによれば、一定の間隔ごとに取る値を繰り返す関数のことを言います。sin関数の場合、周期は2πラジアンです。

周期関数であるsin(θ(t))は、以下のような式が成り立ちます：

sin(θ(t)) = sin(θ(t) + n*2π)

つまり、sin関数においては、ある位相θ(t)と、その位相に2πを整数倍した値を加えたθ(t) + n*2π（ここでnは整数）は、同じ値を返します。

この性質を利用します。具体的には、φには位相の値を2πで割った余りを使用します。これにより、位相がどれだけ大きくなっても、φは常に0から2πの間に収まり、かつsin関数の結果は変わりません。

実装

以上より、位相及び初期位相の計算方法が分かりました。ここからは、初期位相を考慮した実装と、実際に波形がなめらかに接続できるかを検証したいと思います。

以下に、PythonとPyAudioを使用した簡単な音の生成と再生のコードを示します。このコードでは、異なる周波数の音を順番に再生します。連続する次の波の初期位相に、前の波の終点の位相を与えることで、音の切り替え時のノイズを最小限に抑えています。コードを実行すると、音と音の切り替わりが滑らかに聞こえることが確認できるかと思います。

動作環境

Windows 11（Windows Powershell）
Poetry: 1.8.2
python: 3.10
numpy: 1.26.4
PyAudio: 0.2.14

コード

import numpy as np
import pyaudio


A = 0.5  # 振幅
f_list = [493.88, 523.23, 493.88] # 周波数のリスト
play_length = 2 # 音を再生する長さ
fs = 44100      # サンプリングレート

p = pyaudio.PyAudio()

stream = p.open(format=pyaudio.paFloat32,
                channels=1,
                rate=fs,
                output=True)

φ = 0 # 初期位相
t = np.arange(fs * play_length) / fs

for f in f_list:
    θ = 2 * np.pi * f * t + φ  # 位相の計算
    y = A * np.sin(θ)

    stream.write(y.astype(np.float32).tobytes())

    φ = θ[-1] % (2*np.pi) # 初期位相の更新


# 終了処理
stream.stop_stream()
stream.close()
p.terminate()

コードの説明

上記のコードについて簡単に説明します。
まず、φ=0とすることで、最初の波の初期位相を0[rad]に設定しています。

t = np.arange(fs * play_length) / fs

では、tを時刻配列とすることで、後段の処理を効率的に行なえるようにしています。
例えばtの値は[0, 1/fs, 2/fs, ...., fs*play_length/fs]となります。

for f in f_list:
    θ = 2 * np.pi * f * t + φ  # 位相の計算
    y = A * np.sin(θ)

    stream.write(y.astype(np.float32).tobytes())

    φ = θ[-1] % (2*np.pi) # 初期位相の更新

続けてメインの部分です。θでは、先ほど求めた時刻配列tの各要素に対し、位相を計算します。yは、位相（配列）から生成されるsin波であり、生成された波はPyAudioによって再生されます。

この処理を周波数配列の各要素ごとに行うことで、音の高さを切り替えて再生することができます。

ここで、φは一つ前の波の終点における位相を保持しており、次の波の初期位相として設定しています。2πで割った余りを使用しているのは、ステップ3でも説明したように、sin関数の周期性を利用してφの大きさが0~2πに収まるようにするためです。

検証

実際に波が滑らかにつながっているかを確認するために、y1からy2への切り替わりと、y2からy1への切り替わり部分について図示して確認しようと思います。

図3は「y1からy2へ切り替わったとき」と、図4は「y2からy1へ切り替わったとき」について図示した結果になります。y1とy2のつながりをわかりやすく可視化するため、あえて二つの波を一つのグラフに乗せ、切り替わりの前後それぞれ0.004秒ほどの区間を、拡大して図示しています。y1は青色で、y2は赤色です。横軸は時刻samples（秒）であり、縦軸はyの値です。

図3を見ると、音の切り替わり部分である2.00秒付近では、図1とは異なり波が連続につながっているように見えます。また、図4でも、図2と比較するとおおむね波が連続しているように見えます。ただし、図4については、切り替わり部分（t=4.00秒付近）のところで滑らかに繋がっていない部分がわずかにあります。

この原因として考えられることは、初期位相を用いる手法の限界です。例として、y1の波の始まり時刻t = 4.00+1/fs（fs=44100）秒の時を考えます。この時、t=4.00秒以降もy2が続くと仮定すると、t=4.00+1/fs時点のy2の値は、t=4.00時点のy2の値よりも大きくなります。しかし今回の手法では、y1の初期位相をy2の終点の位相に設定しているため、時刻t=4.00+fsにおけるy1の値は、y2の波がそのまま続いた場合よりもわずかに小さくなり、結果として、このようなずれが生じたのだと考えられます。

ただし、1/fs(=44100)は微小な時間であり、y2がそのまま続いた場合とのずれはかなり小さいと思われます。実際、音を聞く限りでもほとんどノイズが確認できないため、本手法が有効に使えるケースは多いと考えられます。

また、図2と図4を比較すると、t=4.00秒地点において、赤い波の終点の位置が異なっています。これは、y2の初期位相を変更したため、その影響によりt=2.00以降の波形が後退したのだと考えられます。

（再掲）図1: y1からy2へ切り替わったとき

図3. y1からy2へ切り替わったとき

（再掲）図2: y2からy1へ切り替わったとき

図4. y1からy2へ切り替わったとき

まとめ

本記事では、時間的に連続した波を再生する際に発生するノイズに対し、初期位相を考慮することでノイズを低減できるかを確認しました。動作確認の結果、音を切り替えたときに生じていたノイズをかなり抑えることができました。また、図示によって、単純に音を切り替えるよりも滑らかに波形が切り替わっていることが確認できました。

今後の検討

本記事は音を切り替えたときのノイズに焦点を当てており、音の再生開始時や終了時のノイズには対処していません。そのため将来的にはこれらの状況にも対処できるようにしたいと考えています。

所感

実装を通じて、初期位相の役割やsin波の周期性について、より深く理解できたと思います。ただ、様々なサイトを参考にさせて頂きつつ、ChatGPTとのやり取りも行ないながら作ってみましたが、真に理解しきれているかは不安があります。今後はできればこういった物理的な知識も増やしたいため、物理とプログラミングをうまく取り込んだ文献などがあれば読みたいと思います。

もし私が認識を間違えている部分や、物理に関する良い参考書などをご存じの方は、ぜひコメントやご指摘を頂けますと幸いです。

参考