Processing math: 100%
メニュー

フーリエ解析

フーリエ解析

フーリエ解析(Fourier Analysis) は,複雑な関数を周波数成分に分解して解析する手法です. 音,電磁波,光などは波としての特徴を持つことから,フーリエ変換(Fourier Transform) によってsin関数やcos関数に分解することが可能であり,それらの関数の周波数成分の強さも求めることができます.

ここで,波を表すsin関数を次のように定義します. ωは角速度(rad/s)であり,1秒間に進むラジアン角の大きさを表します. 例えば,1秒間に1周期だけ進む場合はω=2πとなります. また,aは振幅の最大値を表します.

f(t)=asin(ωt)

例として,1秒間に10周期のsin関数が含まれる10Hzの波を考えます. このとき,ω=102πとなります. また,振幅の最大値はa=1とします.

f(t)=sin(102πt)=sin(20πt)

Image from Gyazo

同様に,1秒間に30周期のsin関数が含まれる30Hzの波を考えます. このとき,ω=302πとなります. また,振幅の最大値はa=2とします.

f(t)=2sin(302πt)=2sin(60πt)

Image from Gyazo

上記の10Hzと30Hzの波を合成すると次のようになります. このように,sin関数やcos関数を合成することで,任意の関数(波)を表現することが可能です.

f(t)=sin(20πt)+2sin(60πt)

Image from Gyazo

このf(t)から,各周波数の振幅を算出するのが フーリエ変換 です. フーリエ変換の結果は次のグラフで表されます. このグラフは パワースペクトル とも呼ばれます. 横軸は周波数であり,10Hzと30Hzの振幅が大きくなっていることがわかります. また,10Hzの振幅に対して,30Hzの振幅が2倍になっていることも確認できます.

Image from Gyazo

先に述べたように音は空気の振動であり波の性質を持ちます. 音の高低(音階)は周波数で決まり,「ド」「レ」「ミ」の周波数は次の表となっています. フーリエ変換して,対応する周波数が取得できるかどうか,後程,確認してみましょう.

音階 周波数
ド(C4) 261.63Hz
レ(E4) 293.67Hz
ミ(F4) 329.63Hz

データセット

ここでは,次のCSV形式のデータセットを利用します. 下記のURLからファイルをダウンロードしてください. 上記の例で用いた10Hzと30Hzのsin関数を合成したデータです. データ数はN=128であり,これは1秒間に128回のサンプリングをしたとこと意味します.

dataset14.csv

また,魔王魂で提供されているピアノの「ド」「レ」「ミ」の音声ファイルをフーリエ変換して周波数成分を取得してみます. 各音声ファイルをダウンロードしておいてください.

Excelでフーリエ解析

Excelの分析ツールを利用してフーリエ解析を試してみましょう. ダウンロードしたファイルをExcelで開いてください.

分析ツール

分析ツールから「フーリエ解析」を選びましょう. 入力範囲にB1:B129を入力します. このとき,「先頭行をラベルとして使用」にチェックを入れます.

Image from Gyazo

次のような結果が算出されます(A1には「fft」と入力しておきます).

Image from Gyazo

上記の算出された値は複素数です. そこで,imabs関数を利用して,B列に振幅を計算します. データ数N=128であることから,0Hzから127Hzまでの周波数に対する振幅が求められることになります.

Image from Gyazo

0Hzから63Hzまでを対象に折れ線グラフを描くと次のようになります(複素共役のためN/2を対象とする). 10Hzと30Hzの振幅が大きくなっていることが確認できます.

Image from Gyazo

Pythonでフーリエ解析

Pythonを利用してフーリエ解析を試してみましょう. Jupyter Labを起動して,Python 3のノートブックを開きます. ノートブックの名前は chapter15.ipynb とします. pandasmatplotlibnumpy, 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)

Image from Gyazo

折れ線グラフで可視化してみましょう. 上述の10Hzと30Hzのsin関数を合成したデータであることが確認できます.

df.plot.line(x="time", y="amplitude")

Image from Gyazo

フーリエ変換

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("振幅")

Image from Gyazo

音声ファイルをフーリエ変換

最後にピアノの「ド」の音声ファイル(do.wav)をフーリエ変換してみましょう. WAVファイルの読み込みにはwavfile.readを利用します.

do = wavfile.read("do.wav") # 261.63Hz

このWAVファイルはステレオになっていますが, 左のチャネル(left_channel)のみを対象とします. データ数は86,272,サンプリングレートは44,100Hzであることから, 86272/441001.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("振幅")

Image from Gyazo

音声データをフーリエ変換します. この結果,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("振幅")

Image from Gyazo

愛知県名古屋市にある椙山女学園大学 文化情報学部 向研究室の公式サイトです. 専門は情報科学であり,人工知能やデータベースなどの技術要素を指導しています. この公式サイトでは,授業で使用している教材を公開すると共に, ベールに包まれた女子大教員のミステリアスな日常を4コマ漫画でお伝えしていきます. サイトに関するご意見やご質問はFacebookまたはTwitterでお問い合わせください.