魔王魂「シャイニングスター」を数式化してみた
この記事はFUN Advent Calendar 2023 の12月12日の枠に投稿したものです。他の方が投稿された記事も是非ご覧ください。
他の方の記事は以下のリンクからアクセスできます。
0. 目次
1. はじめに
1.1. 自己紹介と概要
皆様こんにちは。はじめましての方ははじめまして。
複雑系コース3年のかすまるです。
プロジェクト学習の成果発表会が終わったかと思ったら、グループ報告書、生協の仕事、院試に向けた勉強など重いタスクがまだまだ残っていました。自由の身となれるのはいつになるのでしょうか。
去年は積分と行列を使って、きのこの山よりもたけのこの里の方が優れていることを証明しました。
今年はもう少し高度な数学系ネタに挑戦してみようと思います。
去年の記事↓
皆様は「シャイニングスター」という楽曲をご存知でしょうか。
「シャイニングスター」とは、森田交一様が作詞作曲されたフリーBGMサイト「魔王魂」に掲載されているフリーBGMの一つです。
今では、YouTubeなど様々な場面で使用されている楽曲となっています。
そんな「シャイニングスター」などの楽曲は空気の振動、つまり波によって構成されています。
そして、これらの波は\( \sin \)や\( \cos \)を使って数式で表すことができます。
そこで今回は、「シャイニングスター」を数式にすることで、普段聞いている楽曲を別の形で見てみようと思います。
1.2. お詫び
今回のアドカレのテーマを決めるにあたって、X（旧Twitter）にてアンケートを行いました。
その結果、「YOASOBI『アイドル』を数式化してみた」というテーマが選ばれました。
しかし、数式化するにあたって楽曲の音源をダウンロードする必要がありましたが、著作権やYouTubeの利用規約を侵さずにダウンロードする方法を見つけることができなかったため、楽曲をフリーBGMである「シャイニングスター」に勝手ながら変更しました。
「アイドル」の数式化を楽しみにしてくださった皆様、申し訳ありませんでした。
2. 数式化までの流れ
※わかりやすさを重視して説明するため、数学的に厳密ではない表現が含まれているかもしれません。ご了承ください。
※スマートフォンでは数式が最後まで表示されない場合があります。PCなどの大きな画面で閲覧することを推奨します。
2.1. フーリエ変換と数式化
データから数式を推測する手法はいくつかありますが、今回は音楽などの信号処理の分野で一般的に使われているフーリエ変換を使います。
まず、フーリエ変換の基礎となるフーリエ級数展開について紹介します。
フーリエ級数展開は、未知の数式や扱うことが難しい形の周期的な数式を無理やり\( \sin \)と\( \cos \)の和で表す手法です。
例えば、周期が\( T \)である未知の数式\( f(t) \)にフーリエ級数展開を行うと、
\begin{align} a_k = \frac{2}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) \cos \left( \frac{2 \pi k t}{T} \right) dt \end{align}
\begin{align} b_k = \frac{2}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) \sin \left( \frac{2 \pi k t}{T} \right) dt \end{align}
と定義される係数\( a_k, b_k \)を用いて
\begin{align} f(t) = \frac{a_0}{2} + \sum_{k=1}^{\infty} a_k \cos \left( \frac{2 \pi k t}{T} \right) + b_k \sin \left( \frac{2 \pi k t}{T} \right) \end{align}
となります。
これらの式(1)～(3)がフーリエ級数展開の公式となります。
この式では\( a_k \)と\( b_k \)という2つの係数について計算することになり面倒なので、以下のようなオイラーの公式を用いて係数を1つにします。
\begin{align} e^{i \theta} = \cos \theta + i \sin \theta \end{align}
実際にオイラーの公式(4)を用いると、フーリエ級数展開の公式は
\begin{align} c_k = \frac{1}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) e^{-\frac{2 \pi i k t}{T}} dt \end{align}
と定義される係数\( c_k \)を用いて
\begin{align} f(t) = \sum_{k = -\infty}^{\infty} c_k e^{\frac{2 \pi i k t}{T}} \end{align}
と変形できます。
このように、フーリエ級数展開を使うことで、数式を\( \sin \)と\( \cos \)で表すことができます。
しかし、フーリエ級数展開は周期的な数式しか扱うことができません。
そこで、フーリエ級数展開の公式を周期的でない数式でも扱えるように変形していきます。
まず、周期的でない数式を周期が\( \infty \)の数式として扱います。
また、計算を簡単にするために\( T = 2 l \pi \)とおいて、フーリエ級数展開の公式を変形すると
\begin{align} c_k = \int_{-l \pi}^{l \pi} f(t) e^{-\frac{i k t}{l}} dt \end{align}
\begin{align} f(t) = \frac{1}{2 l \pi} \sum_{k = -\infty}^{\infty} c_k e^{\frac{i k t}{l}} \end{align}
となります。
式(7)について、\( \omega = n / l \)とおいて、\( l \rightarrow \infty \)の極限をとると、\( c_k \)は
\begin{align} F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i \omega t} dt \end{align}
という関数\( F(\omega) \)に変形できます。
このような変形をフーリエ変換と呼びます。
ここで、\( f(t) \)を音楽の波を表す式と考えると、
\begin{align} \omega = \frac{n}{l} = \frac{2 \pi n}{T} \end{align}
と表すことができます。
よって、\( \omega \)は\( n \)の値ごとに決定される\( f(t) \)の周波数と考えることができます。
このことから、音楽を分析するときにおいては、フーリエ変換は音楽の波を周波数ごとの成分に分解して並べたもの、つまり周波数のスペクトルに変換するという操作であると言えます。
また、区分求積法（総和の式\( \sum \)の極限をとると積分になる）を用いて式(8)を変形すると
\begin{align} f(t) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) e^{i \omega t} d \omega \end{align}
となります。
このような関数\( F(\omega) \)から関数\( f(t) \)への変形を逆フーリエ変換と呼びます。
また、音楽を分析するときにおいては、逆フーリエ変換は周波数のスペクトルを元の音楽の波に戻す操作であると言えます。
以上のことを踏まえると、「シャイニングスター」の数式\( f_{\mathrm{out}}(t) \)は、元の「シャイニングスター」の波を\( f_{\mathrm{in}}(t) \)と置くと
\begin{align} F(\omega) = \int_{-\infty}^{\infty} f_{\mathrm{in}}(t) e^{-i \omega t} dt \end{align}
\begin{align} f_{\mathrm{out}}(t) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) e^{i \omega t} d \omega \end{align}
という２つの数式で表すことができます。
2.2. 離散フーリエ変換と数式化
前の節では、音楽の波\( f_{\mathrm{in}}(t) \)が連続した値を取る（\( t \)がどんな値であっても\( f_{\mathrm{in}}(t) \)の値が存在する）という前提で話を進めました。
しかし、実際に得られる波のデータは飛び飛びの値、つまり離散的な値を取ります。
このような離散的な値でもフーリエ変換を行えるように、式を変形していきます。
ここでは、有限な音楽の波のデータを取ったとき、\( n \)番目のデータの値を\( f(n) \)とおきます。
式(5)の連続的な変数\(t, T \)を離散的な変数\(n, N \)と置き換えます。ただし、\( N \)は波のデータの総数です。すると、
\begin{align} c_k = \frac{1}{N} \sum_{n=0}^{N-1} f(n) e^{-\frac{2 \pi i k n}{N}} \end{align}
と変形できます。ここで、係数\( c_k \)を関数\( F(k) \)に置き換えると
\begin{align} F(k) = \frac{1}{N} \sum_{n=0}^{N-1} f(n) e^{-\frac{2 \pi i k n}{N}} \end{align}
となります。この\( F(k) \)は連続的な値のフーリエ変換における周波数のスペクトル\( F(\omega) \)に相当します。
このような離散的な波のデータをスペクトルに変換する操作を離散フーリエ変換と呼びます。
また、式(6)の連続的な変数\(t, T \)を離散的な変数\(n, N \)に、係数\( c_k \)を関数\( F(k) \)と置き換えると、
\begin{align} f(t) = \sum_{k = 0}^{N-1} F(k) e^{\frac{2 \pi i k n}{N}} \end{align}
となります。
このような離散的なスペクトルを波のデータに変換する操作を逆離散フーリエ変換と呼びます。
式(15)と式(16)を変形すると、離散フーリエ変換と逆離散フーリエ変換の式は以下のように表すことができます。
\begin{align} F(k) = \sum_{n=0}^{N-1} f(n) e^{-\frac{2 \pi i k n}{N}} \end{align}
\begin{align} f(t) = \frac{1}{N} \sum_{k = 0}^{N-1} F(k) e^{\frac{2 \pi i k n}{N}} \end{align}
以上のことから、「シャイニングスター」の数式\( f_{\mathrm{out}}(n) \)は、元の「シャイニングスター」の波のデータを\( f_{\mathrm{in}}(n) \)と置くと
\begin{align} F(k) = \sum_{n=0}^{N-1} f_{\mathrm{in}}(n) e^{-\frac{2 \pi i k n}{N}} \end{align}
\begin{align} f_{\mathrm{out}}(n) = \frac{1}{N} \sum_{k = 0}^{N-1} F(k) e^{\frac{2 \pi i k n}{N}} \end{align}
という2つの数式で表すことができます。
3. プログラムによる検証
3.1. 音声データの入手
それでは、実際にこれらの数式に「シャイニングスター」のデータを代入して、数式\( f_{\mathrm{out}}(n) \)を作成していきます。
今回は1番のサビの部分（56.5秒～1分27.5秒）のみを数式化します。
最初に、数式化に必要なデータを入手します。
まず、公式サイトから「シャイニングスター」の音源をダウンロードします。
次に、ダウンロードした音源をAudacityという音楽編集ソフトに取り込みます。
Audacityの公式サイト↓
取り込んだ音源を左耳から聞こえる音と右耳から聞こえる音に分けて、それぞれの時刻ごとの音量をtxtファイルに出力します。
このとき、サンプリング定数は44100Hzとします。
これは、音量を1秒間に44100回txtファイルに出力することを意味します。
この操作によって入手できるデータは大きすぎるため、Pythonを用いてデータを0.01秒単位のデータに圧縮します。
以下のプログラムは、今回使用するモジュールのインポートとデータの取り込み、圧縮を行うプログラムです。
# %%
#ライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt
import os
#各自音楽データを入れたファイルのパスを入力
os.chdir("***")
# %%
#データパスの読み込み
#1番サビ部分：00:56:500~01:27:500
path1l = 'source1l.txt' #Audacityによって出力されたtxtファイル
path1r = 'source1r.txt'
path2l = 'source2l.txt'
path2r = 'source2r.txt'
# %%
#圧縮した音声データを入れるリストの作成
left = []
right = []
# %%
#音声データを圧縮する関数
RANGE = 0.01 # 0.01秒単位の平均を求める
def compress(f1, f2, data):
fprev = 0
for f in [f1, f2]:
list = []
prev = 0
for line in f.read().split("\n"):
line = line.strip()
if line == "":
break
time, dB = line.split("\t")
time = float(time)
try:
dB = float(dB)
except ValueError:
dB = 0 # -ooなので無音
if time > prev + RANGE:
data.append(np.mean(list))
list = []
prev += RANGE
fprev += RANGE
else:
list.append(dB)
#print(line)
# %%
#音声データの圧縮
compress(open(path1l), open(path2l), left)
compress(open(path1r), open(path2r), right)
sizel = len(left)
sizer = len(right)
# %%
#圧縮したデータのプロット
t = np.arange(0, sizel*RANGE, RANGE)
plt.plot(t, left, label="left")
plt.plot(t, right, label="right")
plt.xlabel("time[s]")
plt.ylabel("wave")
plt.legend()
圧縮されたデータの波形は以下の図1のようになります。
ただし、「left」が左耳から聞こえる音の波形、「right」が右耳から聞こえる音の波形を表しています。
このようにして、「シャイニングスター」の音源から波のデータを入手することができました。
3.2. データのフーリエ変換
入手したデータをもとに、実際に離散フーリエ変換を実行していきます。
2章で説明した式(19)の計算です。
Pythonでは、numpyというモジュールに含まれるfft関数によって離散フーリエ変換を行うことができます。
（正確には、fft関数は離散フーリエ変換をより高速に行うための計算方法である高速フーリエ変換を行う関数です。）
以下のプログラムは、fft関数を用いてデータに離散フーリエ変換を行い、その結果得られた周波数のスペクトルを表示するプログラムです。
# %%
#離散フーリエ変換の実施
spl = np.fft.fft(left)
spr = np.fft.fft(right)
freq = np.fft.fftfreq(sizel, d=RANGE)
ampl = abs(spl/(sizel/2))
ampr = abs(spr/(sizer/2))
# %%
#スペクトルのプロット
plt.plot(freq[1 : int(sizel/2)], ampl[1 : int(sizel/2)], label="left")
plt.plot(freq[1 : int(sizer/2)], ampr[1 : int(sizer/2)], label="right")
plt.xlabel("frequency[Hz]")
plt.ylabel("Amplitude")
plt.xscale("log")
plt.legend()
\( \mathrm{freq} \)は、離散フーリエ変換によって生まれた関数\( F(k) \)に紐づいた周波数の配列です。
\( \mathrm{ampl}, \mathrm{ampr} \)は、左耳、右耳から聞こえる音に関する周波数ごとの音量を格納した配列です。
具体的な計算方法は以下のサイトをご覧ください。
そして、以下の図2は、実際に表示された「シャイニングスター」の周波数ごとのスペクトルです。
横軸は周波数の対数軸、縦軸は周波数ごとの音量です。
これによって、数式化に必要な式(19)の計算をして、関数\( F(k) \)を求めることができました。
3.3. スペクトルの逆フーリエ変換
得られたスペクトルに逆離散フーリエ変換を行います。
2章で説明した式(20)の計算です。
Pythonでは、numpyのifft関数によって逆離散フーリエ変換を行います。
以下のプログラムは、実際にスペクトルに対して逆離散フーリエ変換を行って、その結果をグラフに表示させるプログラムです。
# %%
#逆離散フーリエ変換の実施
xl = np.fft.ifft(spl)
xr = np.fft.ifft(spr)
# %%
#変換後のデータのプロット
plt.plot(t, xl.real, label="left")
plt.plot(t, xr.real, label="right")
plt.xlabel("time[s]")
plt.ylabel("wave")
plt.legend()
このプログラムによって、式(20)の\( f_{\mathrm{out}} (n) \)、つまり「シャイニングスター」の数式を得ることができます。
そして、図3は、\( f_{\mathrm{out}} (n) \)をグラフに表した図です。
元のデータの値をグラフ化した図1と数式化した関数\( f_{\mathrm{out}} (n) \)の値をグラフ化した図3を比較すると、ほとんど同じグラフを作成できていることがわかります。
よって、作成した数式が概ね正しいことを確かめることができました。
4. まとめ
4.1. 記事全体のまとめ
今回は、魔王魂さんの「シャイニングスター」を数式化する手法について書いてみました。
難しいところもあったかと思いますが、楽しんでもらえたら幸いです。
来年は時間があるかわかりませんが、もし機会があればまたアドカレの記事作成に挑戦してみたいと思います。
それでは、最後までご覧いただきありがとうございました。
4.2. 12日前後の記事
11日Part1 こたさん 「私と、アートと、それからAI」
（去年は僕の記事への感想ありがとうございました。こたさんはきのこ派でしたか。「先」で待ってますね。）
11日Part2 とみすけさん 「はこだて未来大学に合格した話(他、入りたいサークルの話など)」
（未来大への合格おめでとうございます。キーボード配列自作は強すぎますね。入学後が楽しみです。）
12日Part1 当記事
12日Part2 kazu8823さん 「神ゲー紹介(ガチな方)」
（最近Steamに全く触れていませんでしたが、神ゲーを探しに漁ってみるのも面白そうですね。）
13日Part1 パン屋さん。さん 「高専プロコンで行った手探りチーム開発のお話」
（すでにチーム開発の経験が...。アドカレに参加している新入生のレベル高すぎませんか？）
13日Part2 ろくさん 「留学準備についてか一問一答、part1と統合可能性あり」
（留学も一度は行ってみたいですね。一番の不安要素は海外の食事が口に合うかどうかです。）
