Pythonで信号ファイルを高速フーリエ変換して振幅スペクトルをもとめる

はじめに

大学の研究で信号ファイルをフーリエ変換して振幅スペクトルをみる必要があったので、その過程で学んだことや、詰まった点を初学者目線でまとめようと思います。
自分と同じように初めてFFTのプログラムをつくろうと思ってる人のお役に立てたら幸いです。

高速フーリエ変換(FFT)

今回はCSVファイルで記述される離散的なデータを扱っていくので、離散的フーリエ変換を行います。
離散的フーリエ変換の定義式は

X(m) = \sum_{n=0}^{N-1} x(n) \exp \left(-j \frac{2\pi}{N} mn \right) ~~~(m = 0, 1, 2, ..., N-1)

となり、特定の$n$の値に対するフーリエ変換$X(n)$の計算に$N$回、$m = 0$から$m = N-1$までのフーリエ変換を求めるのに$N^2$回の乗算が必要になります。例えばサンプリング点数$N$が$N = 8,192$の場合、計算回数は67,108,864回になります。

上のように大量になってしまう計算回数を改善するために考えられたのが高速フーリエ(fast Fourier transform, FFT)です。原理について簡単にまとめると、$N$個のデータについて前半と後半に分け、それぞれについてフーリエ変換していくという方法です。最終的に計算回数は$N$回の加算と$N/2$回の乗算に減らせます。$N = 8,192$の場合、FFTでは53,243回の計算ですみます。

FFTを行うデータの長さは計算測度の観点から2の冪乗が望ましいですが、それ以外の長さのデータでも計算をすることはできます。

開発環境

Python	Python3.8.2
numpy	1.22.4
pandas	1.4.2
matplotlib	3.5.2

実際にプログラムを書いてみる

自分が詰まった点

先に自分が苦戦した点についてまとめておきます。特に2つ目は、原因が特定するのに時間がかかってしまいました。

• CSVファイルから自分が抽出したい値だけをフーリエ変換する方法
• numpyのfft関数を使っても出力される値が変わらない

CSVファイルから自分が抽出したい値だけをフーリエ変換する方法

自分が用いたデータは以下のようなCSVファイルです。（冒頭の数点のみ）
最初の数行にサンプリング周波数やデータの単位などの情報があり第５行から測定値が記されています。
また、列方向でも第１列はサンプリング点数の番号、第２第３列は測定値ではなく、第４から第８列が計算に用いたい数値になっています。

Sampling Rate(Hz)	1000
Type	AD	AD	1/1	1/1	1/1	1/1	1/1
Unit			uV	uV	uV	uV	uV
POINT	EVENT		EEG	EEG	EEG	EEG	EEG
00000000	0	529	0.00	0.00	0.00	0.00	0.00
00000001	0	530	0.00	0.00	0.00	0.11	0.00
00000002	0	530	-0.16	-0.21	0.16	0.75	0.16
00000003	0	526	-0.38	-0.43	0.64	2.63	0.80
00000004	0	520	-0.05	-0.21	1.66	6.38	2.68
00000005	0	517	1.13	0.97	3.38	12.55	6.60
00000006	0	517	3.49	3.43	5.63	20.49	12.39
00000007	0	517	6.49	6.65	7.89	29.08	19.10
00000008	0	517	9.66	10.14	9.76	37.23	25.97
00000009	0	517	12.71	13.68	11.16	44.47	32.46
00000010	0	517	15.77	17.22	12.55	50.86	38.52
00000011	0	517	18.88	20.71	14.06	56.33	44.05
00000012	0	516	21.78	23.87	15.56	60.78	48.98
00000013	0	511	24.30	26.61	16.79	64.48	53.22
00000014	0	503	26.82	29.29	18.08	68.29	57.67
00000015	0	502	29.88	32.51	19.74	72.91	62.98

このようなCSVファイルから自分が欲しい数値だけを読み込む方法は以下のようになります。
ファイル名はsample.csvであるとします。

import pandas as pd

df = pd.read_csv(sample.csv, encoding='shift-jis', skiprows=3, usecols=[3])

pd.read_csvでCSVファイルを選択して、skiprowsに最初の何行を飛ばしたいか、usecolsに読み込みたい列を入力します。この数字は0番目から始まることに注意してください。encodingは文字コードの選択です。

最初の行から計算に用いる数値が並んでいる場合は、

df = pd.read_csv(sample.csv, encoding='shift-jis', usecols=[3], header=None)

のように header = None と指定すればよいです。

numpyのfft関数を使っても出力される値が変わらない

FFTの実装にはnumpyのnumpy.fft.fft関数を使いました。失敗した時のコードは以下です。

import numpy as np
import pandas as pd

df = pd.read_csv(sample.csv, encoding='shift-jis', skiprows=3, usecols=[3])
F = np.fft.fft(df)

このコードでは、fft関数にデータを入力しても変わらない値が出力されてしまいました。（ちなみにネットで調べると多くがこのコードとおなじものでした。）

調べていくと、このコードがダメな理由はデータの型にありました。上のコードのdfという配列データは(5000, 1)の２次元配列でした。この２次元配列を行方向すなわち要素数１の方向にフーリエ変換していたというわけです...。
Numpyのマニュアルでも、fft関数の引数にフーリエ変換する軸方向を指定する axis という値が明記されていました。したがって、正しいコードは以下のようになります。

df = pd.read_csv(sample.csv, encoding='shift-jis', skiprows=3, usecols=[3])
F = np.fft.fft(df, axis=0)

振幅スペクトルの計算

振幅スペクトルとは、フーリエ変換した$X(m)$の値の絶対値$\left|X(m)\right|$です。したがって、入力値の絶対値を出力してくれるnp.abs関数によって求めれます。

import numpy as np
import pandas as pd

df = pd.read_csv(sample.csv, encoding='shift-jis', skiprows=3, usecols=[3])
F = np.fft.fft(df)

# 振幅スペクトルAmp
Amp = np.abs(F)

自分が書いたコード

自分の研究では脳波信号を用いてFFTを実装しており、その過程で得られた結果の加算平均などをとっていますが、平均をとる必要がない場合はデータ配列を一つにすれば良いです。
また、データのCSVファイルが複数あったのでsys.argvで引数を読み込んでいます。

fft.py

"""
python3 fft.py (csvファイル 複数可)
サンプリング周波数: 1000Hz
電極は P3, P4, O1, O2, Pz の５つ
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys


# 脳波データ取得
def import_data(i):
	df_P3 = pd.read_csv(sys.argv[i], encoding='shift-jis', skiprows=3, usecols=[3])
	df_P4 = pd.read_csv(sys.argv[i], encoding='shift-jis', skiprows=3, usecols=[4])
	df_O1 = pd.read_csv(sys.argv[i], encoding='shift-jis', skiprows=3, usecols=[5])
	df_O2 = pd.read_csv(sys.argv[i], encoding='shift-jis', skiprows=3, usecols=[6])
	df_Pz = pd.read_csv(sys.argv[i], encoding='shift-jis', skiprows=3, usecols=[7])

	df = [df_P3, df_P4, df_O1, df_O2, df_Pz]
	return df


# fft計算
def fft(index):
	df = import_data(index)
	data_length = len(df[:][0])
	Amp = np.zeros((data_length, 1))

	F_P3 = np.fft.fft(df[0], axis=0)
	F_P4 = np.fft.fft(df[1], axis=0)
	F_O1 = np.fft.fft(df[2], axis=0)
	F_O2 = np.fft.fft(df[3], axis=0)
	F_Pz = np.fft.fft(df[4], axis=0)

	Amp += np.abs(F_P3)
	Amp += np.abs(F_P4)
	Amp += np.abs(F_O1)
	Amp += np.abs(F_O2)
	Amp += np.abs(F_Pz)

	Amp /= 5 * (data_length / 2)

	return Amp


def main():
	data_num = len(sys.argv) - 1
	df_sample = import_data(1)
	data_length = len(df_sample[:][0])
	freq = np.linspace(0, 1000, data_length)

	index = 1
	Amp = np.zeros((data_length, 1))
	while index <= data_num:
		Amp += fft(index)
		index += 1

	Amp /= data_num
	
	# グラフ表示
	plt.figure()
	plt.rcParams['figure.subplot.bottom'] = 0.2
	plt.plot(freq, Amp)
	plt.xticks([0, 5, 10, 12, 15, 20, 30, 50])
	plt.xlabel('Frequency [Hz]', fontsize=20)
	plt.ylabel('Amplitude [μV]', fontsize=20)
	plt.xlim(5, 30)
	plt.ylim(0, 2.5)
	plt.grid()
	plt.show()


if __name__ == "__main__":
	main()

このコードに実験のデータファイルを入力してみると以下の図が得られました。
15Hzの周波数成分が多いか確認したかったため、目的が達成できました。

まとめ

この記事ではCSVファイルのデータを高速フーリエ変換する方法について記しました。
今回解決策が見つかるのに時間がかかってしまった原因としては、自分でマニュアルを読むことをせず、ネットにあるFFTのやり方を参考にプログラムを作成し続けたからだと反省しました。以後は、まず公式のマニュアルに目を通し、自力でプログラムを書けるようにトレーニングしていきたいと思います。
最後まで読んでいただきありがとうございました。