はじめに
フーリエ級数展開とはある周期をもった波形をsinとcosの足し合わせで表そうというものです
これによってある信号において、周波数が異なるsin波とcos波がそれぞれどれくらいの成分で含まれているか知ることができます。
実用例としては、ある特定の周波数の音をカットしたり、画像の圧縮などに使われたりします。
今回はこの技術を使って、ピコピコ音と称されるノコギリ波の音波を生成してみましょう
フーリエ級数展開
フーリエ級数展開ではある周期信号 f(t)をsinとcosで表します。角周波数をω, sinとcosの振幅をそれぞれa, b、またオフセットをcとすると下記の用に表現できます
f(x) = \sum_{n = 1}^\infty (a_n \sin n\omega t + b_n \cos n\omega t) + c
フーリエ級数展開をするには、各角周波数における振幅とオフセットを求めれば良いことになります
オフセットを求める
sin, cosともに角周波数がどんな値でも周期Tで積分をすると0になります
ゆえにf(t)を一周期分積分した値は幅T, 高さcの長方形の面積となります
以上よりcは次の用に求められます
\int_{0}^T f(t) dt = T \times c \Leftrightarrow c = \frac{1}{T} \int_0 ^T f(t) dt
角周波数nωtにおけるsinの振幅を求める
三角関数の直交性より、信号f(t)に対して振幅が1で取り出したい角周波数を掛けたものを積分すると取り出したい波以外の波の面積が0になります
もし詳細に知りたい方は、図を交えてわかりやすく書かれているので、こちらの書籍を参照ください。
また取り出された波の成分の面積は、振幅 x T/2となるので、振幅の値を求めることができます
\int_{0}^T f(t) \sin n\omega t dt = \frac{a_nT}{2} \Leftrightarrow a_n = \frac{2}{T} \int_0 ^T f(t) \sin n\omega t dt
角周波数nωtにおけるcosの振幅を求める
cosもsinと同様に計算ができるので以下のようになります
b_n = \frac{2}{T} \int_0 ^T f(t) \cos n\omega t dt
ノコギリ波をフーリエ級数展開する
ノコギリ波は $f(t) = \frac{A}{T}t \, (0 \leq t < T)$ で表されるので、この関数に対してフーリエ級数展開を行ってみましょう
各係数について、計算すると下記のようになります
\begin{align}
c &=\frac{1}{T} \int_0 ^T f(t) dt =\frac{1}{T} \int_0 ^T \frac{A}{T} t \,dt \\
&=\frac{A}{T^2} \Big[t^2\Big]^T_0 \\
&=\frac{A}{2}
\end{align}
\begin{align}
a_n &= \frac{2}{T} \int_0 ^T f(t) \sin n\omega t dt = \frac{2}{T} \int_0 ^T \frac{A}{T}t \sin n\omega t dt \\
&= \frac{2A}{T^2}\int_0^T t\sin n \omega t dt = \frac{2A}{T^2}\Big[ - \frac{t \cos n \omega t}{n \omega} + \frac{\sin n \omega t}{(n \omega)^2} \Big]^T_0 \\
&= - \frac{A}{n \pi} \, (\because \omega = \frac{2 \pi}{T})
\end{align}
\begin{align}
b_n &= \frac{2}{T} \int_0 ^T f(t) \cos n\omega t dt = \frac{2}{T} \int_0 ^T \frac{A}{T}t \cos n\omega t dt \\
&= \frac{2A}{T^2}\int_0^T t\cos n \omega tdt = \frac{2A}{T^2}\Big[ \frac{t \sin n \omega t}{n \omega} + \frac{\cos n \omega t}{(n \omega)^2} \Big]^T_0 \\
&= 0 \, (\because \omega = \frac{2 \pi}{T})
\end{align}
したがって、このノコギリ波は以下のように表現することができます
f(t) = -\frac{A}{\pi}\sum_{n=1}^\infty \frac{1}{n} \sin n\omega t + \frac{A}{2}
グラフでノコギリ波を表示してみる
上記で求められた結果がノコギリ波になっているか描画して見ましょう
コンピュータで無限は扱うことができないので、有限のnの数を設定してノコギリ波として近似します。
pythonでグラフを描画するコードを書くと以下のような感じ
def sawtooth_wave(N):
dt = 0.01
T = 2 * np.pi / 4.
A = 1
omega = 2 * np.pi / T
t = np.arange(0, 2*np.pi, 0.01)
y = np.zeros(len(t))
for n in range(1,N):
y += - A / (np.pi * n) * np.sin( n * 2 * np.pi / T * t)
y += A / 2.
plt.plot(t, y, label="n = {}".format(N))
plt.yticks([0.5], ["A/2"])
plt.xticks([T, 2*T, 3*T, 4*T], ["T", "2T", "3T", "4T"])
plt.legend(loc="upper left")
下記にn=0, 1, 5, 100とした場合のノコギリ波を描画してみます
図を見て分かるようにnを増やしていくと、ノコギリ波に近づいていくことがわかります
ノコギリ波を実際に聞いてみる
ノコギリ波がどんな音が聞こえるかローカルに保存してみましょう
import numpy as np
# 定数
SR = 44100 # サンプリング周波数
f = 440. # 周波数
A = 1 # 振幅
T = 5 # 音声の時間
N = 5000 # sin波の重ねる数
# 各時間
t = np.linspace(0, T, SR * T)[:-1]
# ノコギリ波の計算
y = np.zeros(len(t))
for n in range(1,N):
y += - A / (np.pi * n) * np.sin( n * 2 * np.pi * f * t)
y += A / 2.
# 16bitで量子化
y2 = (y * np.iinfo(np.int16).max).astype(np.int16)
# ステレオに変換
y3 = np.array(np.c_[y2,y2])
# 保存
import scipy.io.wavfile as siw
siw.write("sawtooth.wav", SR, y3)
Qiitaにオーディオファイルが貼ることができないので、はてぶの方か自分で上記のプログラムを実行してみて確認してみてください
pythonで音を扱う~フーリエ級数展開で作るノコギリ波~ - すからすっからすっからかん
波形がギザギザしてると、電子音っぽく聞こえますね
まとめ
フーリエ級数展開を紹介し、ノコギリ波を実際にフーリエ級数展開してみた
以上！