フーリエ解析
フーリエ解析
フーリエ解析(Fourier Analysis) は,複雑な関数を周波数成分に分解して解析する手法です. 音,電磁波,光などは波としての特徴を持つことから,フーリエ変換(Fourier Transform) によってsin関数やcos関数に分解することが可能であり,それらの関数の周波数成分の強さも求めることができます.
ここで,波を表すsin関数を次のように定義します. ωは角速度(rad/s)であり,1秒間に進むラジアン角の大きさを表します. 例えば,1秒間に1周期だけ進む場合はω=2πとなります. また,aは振幅の最大値を表します.
f(t)=a⋅sin(ωt)
例として,1秒間に10周期のsin関数が含まれる10Hzの波を考えます. このとき,ω=10⋅2πとなります. また,振幅の最大値はa=1とします.
f(t)=sin(10⋅2πt)=sin(20πt)
同様に,1秒間に30周期のsin関数が含まれる30Hzの波を考えます. このとき,ω=30⋅2πとなります. また,振幅の最大値はa=2とします.
f(t)=2⋅sin(30⋅2πt)=2⋅sin(60πt)
上記の10Hzと30Hzの波を合成すると次のようになります. このように,sin関数やcos関数を合成することで,任意の関数(波)を表現することが可能です.
f(t)=sin(20πt)+2⋅sin(60πt)
このf(t)から,各周波数の振幅を算出するのが フーリエ変換 です. フーリエ変換の結果は次のグラフで表されます. このグラフは パワースペクトル とも呼ばれます. 横軸は周波数であり,10Hzと30Hzの振幅が大きくなっていることがわかります. また,10Hzの振幅に対して,30Hzの振幅が2倍になっていることも確認できます.
先に述べたように音は空気の振動であり波の性質を持ちます. 音の高低(音階)は周波数で決まり,「ド」「レ」「ミ」の周波数は次の表となっています. フーリエ変換して,対応する周波数が取得できるかどうか,後程,確認してみましょう.
音階 | 周波数 |
---|---|
ド(C4) | 261.63Hz |
レ(E4) | 293.67Hz |
ミ(F4) | 329.63Hz |
データセット
ここでは,次のCSV形式のデータセットを利用します. 下記のURLからファイルをダウンロードしてください. 上記の例で用いた10Hzと30Hzのsin関数を合成したデータです. データ数はN=128であり,これは1秒間に128回のサンプリングをしたとこと意味します.
また,魔王魂で提供されているピアノの「ド」「レ」「ミ」の音声ファイルをフーリエ変換して周波数成分を取得してみます. 各音声ファイルをダウンロードしておいてください.
Excelでフーリエ解析
Excelの分析ツールを利用してフーリエ解析を試してみましょう. ダウンロードしたファイルをExcelで開いてください.
分析ツール
分析ツールから「フーリエ解析」を選びましょう. 入力範囲にB1:B129を入力します. このとき,「先頭行をラベルとして使用」にチェックを入れます.
次のような結果が算出されます(A1には「fft」と入力しておきます).
上記の算出された値は複素数です.
そこで,imabs
関数を利用して,B列に振幅を計算します.
データ数N=128であることから,0Hzから127Hzまでの周波数に対する振幅が求められることになります.
0Hzから63Hzまでを対象に折れ線グラフを描くと次のようになります(複素共役のためN/2を対象とする). 10Hzと30Hzの振幅が大きくなっていることが確認できます.
Pythonでフーリエ解析
Pythonを利用してフーリエ解析を試してみましょう.
Jupyter Labを起動して,Python 3のノートブックを開きます.
ノートブックの名前は chapter15.ipynb とします.
pandas
,matplotlib
,numpy
, SciPy
などのライブラリをインポートしておきましょう.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from scipy.io import wavfile
データセットの読込
read_csv
関数でデータセットを読み込みます.
df = pd.read_csv("dataset14.csv")
display(df)
折れ線グラフで可視化してみましょう. 上述の10Hzと30Hzのsin関数を合成したデータであることが確認できます.
df.plot.line(x="time", y="amplitude")
フーリエ変換
np.fft.fft
でフーリエ変換を適用します.
フーリエ変換された結果を折れ線グラフで表示します.
Excelの結果と同様に,10Hzと30Hzの振幅が大きくなっていることが確認できます.
fft_data = np.abs(np.fft.fft(df["amplitude"]))
plt.plot(fft_data)
plt.xlim(0, 63)
plt.xlabel("周波数")
plt.ylabel("振幅")
音声ファイルをフーリエ変換
最後にピアノの「ド」の音声ファイル(do.wav)をフーリエ変換してみましょう.
WAVファイルの読み込みにはwavfile.read
を利用します.
do = wavfile.read("do.wav") # 261.63Hz
このWAVファイルはステレオになっていますが,
左のチャネル(left_channel
)のみを対象とします.
データ数は86,272,サンプリングレートは44,100Hzであることから,
86272/44100≃1.96秒の音声データであることがわかります.
rate = do[0]
left_channel = do[1][:,0]
right_channel = do[1][:,1]
size = len(left_channel)
print(size) # データ数
print(rate) # サンプリングレート
print(left_channel[0:10]) # 先頭10個のデータを表示
86272
44100
[-327680 -327680 -196608 -262144 -393216 -196608 -131072 -131072 65536
196608]
音声データを折れ線グラフで可視化します. 振幅が徐々に減衰する波のデータであることがわかります.
second = np.linspace(0, (size / rate), size) # 横軸
plt.plot(second, left_channel)
plt.xlabel("時間(秒)")
plt.ylabel("振幅")
音声データをフーリエ変換します. この結果,261Hz付近に大きな振幅が確認できます. また,その2倍の522Hz,3倍の783Hzにも振幅が表れています. これは 倍音 と呼ばれる現象で,倍音の違いが楽器や歌声の音色を決めると言われます.
fft_data = np.abs(np.fft.fft(left_channel))
fft_freq = np.fft.fftfreq(size, d=1.0/rate)
plt.plot(fft_freq, fft_data)
plt.xlim(0, 1000)
plt.xlabel("周波数(Hz)")
plt.ylabel("振幅")