こんにちは、うさじん（@Rabi_Jin_）です。
今回はデータ分析シリーズとなります。
タイトル通り、PythonのNumpyを使って周波数分析をやってみました。
高速フーリエ変換（FFT）は波形が持つ周波数成分を分析するアルゴリズムです。
コードをつくる
最終的な出力は以下となります。
1つ目のグラフに正弦波、2つ目のグラフにFFTを使ったパワースペクトルを表示します。
Pythonのmatplolibを使ったグラフの表示については以下の記事でご紹介しています。
こんにちは、うさじんです。 本記事では前回、正弦波をpandas.DataFrame.plotでグラフを作りましたが、 それと同じことをmatplotlibのsubplotsで行ってみました。 前回の記[…]
こんにちは、うさじんです。 前回、Pythonのmatplotlibでグラフシートを作りました。 matplotlibはDataFrameとことなり、色々な設定項目があるのが良いですね。 本記事では、matplotlibの[…]
コードの作成のは以下のステップで勧めます。
- 正弦波を作る
- FFT（高速フーリエ変換）で分析する
- 各種分析値を抽出する
- グラフをつくる
- コード全体
正弦波をつくる
FFTを行うための模擬信号を作ります。
そのため、まずはじめに各種ライブラリをインポートします。
from matplotlib import pyplot as plt
import numpy as np
import japanize_matplotlib次に正弦波を決める条件を設定し、時刻列データ及び数値データを作成します。
n= 1024 #データ数
Fs = 1024 #サンプリング周波数
dt = 1/Fs #サンプリング周期
Frange = Fs/2.56 #分析周波数
F = 50 #正弦波の周波数
gain = 1 #正弦波の振幅
offset = 0 #直流成分
#時間軸データ
#０～１秒までを表示
time = np.arange(0, 1, dt)
#正弦波データ
#振幅を１、オフセットを0とします。
sine = gain * np.sin(2*np.pi*F*time) + offset1行目 データ数が1024なのはFFTが２のべき乗のデータ数を解析するアルゴリズムであるからです。
FFT（高速フーリエ変換）で分析する
先程作成した振幅のデータをFFTで処理します。
#高速フーリエ変換
FFT_sine = np.fft.fft(sine)
FFT_sine = abs(FFT_sine * 2 / n) #交流成分の正規化（2倍してサンプル数で割り算して絶対値を取る）
FFT_sine[0] = abs(FFT_sine[0] / 2) #直流成分の正規化（サンプル数（2）で割り算して絶対値を取る）
#周波数軸
Ffreq = np.fft.fftfreq(n, d = dt) #データ数、サンプリング周期を指定する。注意点として、numpyのFFTアルゴリズムで返されるのは複素数（実数と虚数で表す値）となります。
そのままグラフで表示しても、振幅値がなんのこっちゃわからない数字が表示されます。（下記）
まず、正弦波の振幅値と合わせるため、サンプル数で除し（正規化）、2倍します（鏡像部の足し込み）。
その後、FFTの戻り値である複素数の絶対値をとります。
また、FFTを行う際はサンプリング周波数の半分のナイキスト周波数からサンプリング周波数までの数値は使用しません。
今回はサンプリング周波数のFs:1024Hzの分析周波数であるFrange:400Hzを横軸とします。
各種分析値を抽出する
波形の数値を確認するため、作成した正弦波およびFFTのデータから分析値を抽出します。
#正弦波データの最大最小値
maxmin = "max:" + str(max(sine)) + " " + "min:" + str(min(sine)) #最大値最小値の取得
#FFT最大値と周波数の取得
Fmax = max(FFT_sine) #FFTの最大値取得
maxHz = np.argmax(FFT_sine) #FFTの最大値周波数位置の取得
fftmax = "peak:" + str(round(Fmax, 1)) + " " + "@" + str(maxHz) + "Hz"今回はFFT時の周波数分解能⊿Fが1Hzのため、6行目で返すmaxHzがそのまま周波数値となります。
np.argmaxはリストの番号を返すもののため、サンプリング周波数が1024以外の場合、Ffreqのリストから抽出しましょう。
print(Ffreq[maxHz])取得した値をグラフに挿入します。
グラフをつくる
グラフを作成します。
#Figureの作成
fig = plt.figure(figsize=(6, 8),tight_layout=True) #figの設定、余分な余白削除
#グラフタイトルの作成
fig.suptitle("正弦波:" + str(F) + "Hz", weight=5, fontsize=16) #グラフタイトルの設定、太さ、サイズ
#データ及び軸項目の作成
#正弦波のプロット
plt.subplot(2,1,1) #サブプロットの作成、位置
plt.title("時刻歴グラフ "+maxmin, loc='left') #グラフタイトルの設定、位置
plt.plot(time,sine) #グラフデータの指定（x, y）
plt.xlim(0,1) #横軸の最小値、最大値
plt.xlabel("Time[sec]") #横軸ラベル
plt.ylabel("振幅[-]") #縦軸ラベル
#FFTのプロット
plt.subplot(2,1,2) #サブプロットの作成、位置
plt.title("パワースペクトルグラフ " + fftmax, loc='left') #グラフタイトルの設定、位置
plt.plot(Ffreq[0:int(Frange)], FFT_sine[0:int(Frange)]) #グラフデータの指定（x, y）
plt.xlim(0, Frange) #横軸の最小値、最大値
plt.xlabel("Frequency[Hz]") #横軸ラベル
plt.ylabel("振幅[-]") #縦軸ラベルFFTのプロット時、分析周波数以上をプロットしないよう、配列をスライスします。
また、先程抽出した各種分析値もグラフに入れ込みます。
コード全体
コード全体と実行結果です。
from matplotlib import pyplot as plt
import numpy as np
import japanize_matplotlib
n= 1024 #データ数
Fs = 1024 #サンプリング周波数
dt = 1/Fs #サンプリング周期
Frange = Fs/2.56 #分析周波数
F = 50 #正弦波の周波数
gain = 1 #正弦波の振幅
offset = 0 #直流成分
#時間軸データ #０～１秒までを表示
time = np.arange(0, 1, dt)
#正弦波データ
#振幅を１、オフセットを0とします。
sine = gain * np.sin(2*np.pi*F*time) + offset
#高速フーリエ変換
FFT_sine = np.fft.fft(sine)
FFT_sine = abs(FFT_sine * 2 / n) #交流成分の正規化（2倍してサンプル数で割り算して絶対値を取る）
FFT_sine[0] = abs(FFT_sine[0] / 2) #直流成分の正規化（サンプル数（2）で割り算して絶対値を取る）
#周波数軸
Ffreq = np.fft.fftfreq(n, d = dt) #データ数、サンプリング周期を指定する。
#正弦波データの最大最小値
maxmin = "max:" + str(max(sine)) + " " + "min:" + str(min(sine))
#最大値最小値の取得 #FFT最大値と周波数の取得
Fmax = max(FFT_sine) #FFTの最大値取得
maxHz = np.argmax(FFT_sine) #FFTの最大値周波数位置の取得
fftmax = "peak:" + str(round(Fmax, 1)) + " " + "@" + str(maxHz) + "Hz"
#Figureの作成
fig = plt.figure(figsize=(6, 8),tight_layout=True) #figの設定、余分な余白削除
#グラフタイトルの作成
fig.suptitle("正弦波:" + str(F) + "Hz", weight=5, fontsize=16) #グラフタイトルの設定、太さ、サイズ
#データ及び軸項目の作成
#正弦波のプロット
plt.subplot(2,1,1) #サブプロットの作成、位置
plt.title("時刻歴グラフ "+maxmin, loc='left') #グラフタイトルの設定、位置
plt.plot(time,sine) #グラフデータの指定（x, y）
plt.xlim(0,1) #横軸の最小値、最大値
plt.xlabel("Time[sec]") #横軸ラベル
plt.ylabel("振幅[-]") #縦軸ラベル
#FFTのプロット
plt.subplot(2,1,2) #サブプロットの作成、位置
plt.title("パワースペクトルグラフ " + fftmax, loc='left') #グラフタイトルの設定、位置
plt.plot(Ffreq[0:int(Frange)], FFT_sine[0:int(Frange)]) #グラフデータの指定（x, y）
plt.xlim(0, Frange) #横軸の最小値、最大値
plt.xlabel("Frequency[Hz]") #横軸ラベル
plt.ylabel("振幅[-]") #縦軸ラベル正弦波の振幅は0-pで1です。
FFTのパワースペクトルでも50Hzに1を確認しました。
おわりに
以上、高速フーリエ変換（FFT）による周波数分析でした。
波形データのFFT分析はExcelでもできますが、
Pythonの方がライブラリも豊富でより詳細に分析できます。
次はランダム信号や複数の周波数を持つ信号に対しても上記のFFTを使って分析してみたいですね。
しっかりした分析ができるか確認したいです。
よろしくお願いします。
Pythonを使ったデータ処理にご興味の方は下記の本がとても参考になります。
|
［エンジニアのための］データ分析基盤入門 データ活用を促進する！ プラットフォーム＆データ品質の考え方 [ 斎藤 友樹 ] 価格:2,992円 |
|
Pythonによるあたらしいデータ分析の教科書 第2版【電子書籍】[ 寺田学 ] 価格:2,838円 |
|
Pythonによるデータ分析入門 第2版 NumPy、pandasを使ったデータ処理ー [ Wes McKinney ] 価格:4,180円 |