Matplotlib Python データ解析モモノキ＆ナノネと学習高速フーリエ変換FFT

Pythonで高速フーリエ変換（FFT）の練習-5 周波数ピークを自動で検出

Pythonで高速フリーエ変換（FFT）を行う方法をモモノキ＆ナノネと一緒に学習していきます。

モモノキ＆ナノネと一緒にPythonでFFTの使い方を覚えよう（5）

極大値と極小値の取得方法を覚えてピークの自動検出に挑戦しよう

高速フーリエ変換（FFT）練習の続きを始めるよ。今回はFFT分析結果からピークの周波数を自動で検出してみよう。

グラフを見れば、周波数はだいだいわかるんじゃない？

でも、FFT分析結果のグラフにピークがマーキングできたりしたらチョットうれしくない？

そうかなぁ？

簡易的な自動ピーク検出だけど、先に結果を出してみるね。たとえばこんな感じでできるよ。

In [1]:

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

N = 512 # サンプル数
dt = 0.01 # サンプリング間隔
fq1, fq2 = 3, 8 # 周波数
amp1, amp2 = 1.5, 1 # 振幅 

t = np.arange(0, N*dt, dt) # 時間軸
# 時間信号作成
f = amp1*np.sin(2*np.pi*fq1*t)
f += amp2*np.sin(2*np.pi*fq2*t)
f += np.random.randn(N)*0.5 # ノイズ
freq = np.linspace(0, 1.0/dt, N) # 周波数軸

# 高速フーリエ変換
F = np.fft.fft(f)

F_abs = np.abs(F) # 複素数 ->絶対値に変換
# 振幅を元の信号のスケールに揃える
F_abs = F_abs / (N/2) # 交流成分
F_abs[0] = F_abs[0] / 2 # 直流成分
    
# グラフ表示（時間軸）
plt.figure(figsize=(8,6))

plt.subplot(211)
plt.plot(t, f)
plt.xlabel('Time(s)')
plt.ylabel('Signal')

# FFTデータからピークを自動検出
maximal_idx = signal.argrelmax(F_abs, order=1)[0] # ピーク（極大値）のインデックス取得
# ピーク検出感度調整もどき、後半側（ナイキスト超）と閾値より小さい振幅ピークを除外
peak_cut = 0.3 # ピーク閾値
maximal_idx = maximal_idx[(F_abs[maximal_idx] > peak_cut) & (maximal_idx <= N/2)]

# グラフ表示（周波数軸）
plt.subplot(212)
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

plt.axis([0,1.0/dt/2,0,max(F_abs)*1.5])
plt.plot(freq, F_abs)
plt.plot(freq[maximal_idx], F_abs[maximal_idx],'ro')

# グラフにピークの周波数をテキストで表示
for i in range(len(maximal_idx)):
    plt.annotate('{0:.0f}(Hz)'.format(np.round(freq[maximal_idx[i]])),
                 xy=(freq[maximal_idx[i]], F_abs[maximal_idx[i]]),
                 xytext=(10, 20),
                 textcoords='offset points',
                 arrowprops=dict(arrowstyle="->",connectionstyle="arc3,rad=.2")
                )

plt.subplots_adjust(hspace=0.4)
plt.show()
print('peak', freq[maximal_idx])

peak [ 2.93542074  8.02348337]

同じコードを使って信号をもう少し複雑すると、こんな感じ。

In [2]:

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

N = 512 # サンプル数
dt = 0.01 # サンプリング間隔
fq1, fq2 , fq3 = 12, 24, 33 # 周波数
amp1, amp2, amp3 = 1.5, 1, 2 # 振幅 

t = np.arange(0, N*dt, dt) # 時間軸
# 時間信号作成
f = amp1*np.sin(2*np.pi*fq1*t)
f += amp2*np.sin(2*np.pi*fq2*t)
f += amp3*np.sin(2*np.pi*fq3*t)
f += np.random.randn(N)*0.3 # ノイズ
freq = np.linspace(0, 1.0/dt, N) # 周波数軸

# 高速フーリエ変換
F = np.fft.fft(f)

F_abs = np.abs(F) # 複素数 ->絶対値に変換
# 振幅を元の信号のスケールに揃える
F_abs = F_abs / (N/2) # 交流成分
F_abs[0] = F_abs[0] / 2 # 直流成分
    
# グラフ表示（時間軸）
plt.figure(figsize=(8,6))

plt.subplot(211)
plt.plot(t, f)
plt.xlabel('Time(s)')
plt.ylabel('Signal')

# FFTデータからピークを自動検出
maximal_idx = signal.argrelmax(F_abs, order=1)[0] # ピーク（極大値）のインデックス取得
# ピーク検出感度調整もどき、後半側（ナイキスト超）と閾値より小さい振幅ピークを除外
peak_cut = 0.3 # ピーク閾値
maximal_idx = maximal_idx[(F_abs[maximal_idx] > peak_cut) & (maximal_idx <= N/2)]

# グラフ表示（周波数軸）
plt.subplot(212)
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

plt.axis([0,1.0/dt/2,0,max(F_abs)*1.5])
plt.plot(freq, F_abs)
plt.plot(freq[maximal_idx], F_abs[maximal_idx],'ro')

# グラフにピークの周波数をテキストで表示
for i in range(len(maximal_idx)):
    plt.annotate('{0:.0f}(Hz)'.format(np.round(freq[maximal_idx[i]])),
                 xy=(freq[maximal_idx[i]], F_abs[maximal_idx[i]]),
                 xytext=(10, 20),
                 textcoords='offset points',
                 arrowprops=dict(arrowstyle="->",connectionstyle="arc3,rad=.2")
                )

plt.subplots_adjust(hspace=0.4)
plt.show()
print('peak', freq[maximal_idx])

peak [ 11.93737769  24.0704501   33.07240705]

ちょっとイイかも。

どうすればピークを検出できるの？

FFT解析データからピークを検出したい場合、グラフの極大値が分かれば目的を達成できそうだよね。

モモノキ、極大値って何箇所もあるけど簡単に求まるの？

scipyを使うと極大値や極小値を自動で求めてくれるよ。

自動で簡単ならやってみる。

FFT結果からピークを検出する例は、上のコード例を参考に自分で試してみて。今回は極大値と極小値を取得するキホン練習だけやってみよう。

OK、FFTのピーク検出はあとで試してみる。極大値と極小値の取得を練習する。

まずはデータ準備。ナノネ、乱数でノイズを入れた正弦波を一つ作ってグラフ表示してみて。

In [3]:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

N = 64 # データ数
x = np.linspace(0, 2*np.pi, N) # x軸
y = np.sin(x) # 正弦波
y_add_noise = y + np.random.randn(N)*0.3 # 正弦波＋ノイズ（ピーク検出テスト用）

plt.plot(x, y,'r',label='sin(x)', c='gray', alpha=0.5) # 正弦波
plt.plot(x, y_add_noise,label='sin(x)+noise', c='g') # 正弦波＋ノイズ

plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()

Out[3]:

<matplotlib.legend.Legend at 0xab9fdf6c>

ノイズ入り正弦波を作った（ノイズなし正弦波も、おまけでプロット）。

続いて、scipyのsignalをインポートしてね。極大値はsignal.argrelmaxで、極小値はsignal.argrelminを使って、それぞれ求めることができるよ。

In [4]:

from scipy import signal

signalのインポートが完了した。使い方は？

最初は極大値の取得方法。signal.argrelmax(numpy配列データ, order=1)とすれば、極大値のインデックスが配列で返ってくるよ。

ナノネ、さっき作ったノイズ入り正弦波から極大値を取得を試してみて。インデックスの値を利用して極大値にところにマーカーのプロットもできる？

In [5]:

# ピーク検出（極大値、極小値）
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

N = 64 # データ数
x = np.linspace(0, 2*np.pi, N) # x軸
y = np.sin(x) # 正弦波
y_add_noise = y + np.random.randn(N)*0.3 # 正弦波＋ノイズ（ピーク検出テスト用）

# ピーク検出
maximal_idx = signal.argrelmax(y_add_noise, order=1) # 極大値インデックス取得

plt.plot(x, y,'r',label='sin(x)', c='gray', alpha=0.5) # 正弦波
plt.plot(x, y_add_noise,label='sin(x)+noise', c='g') # 正弦波＋ノイズ

plt.plot(x[maximal_idx],y_add_noise[maximal_idx],'ro',label='peak_maximal') # 極大点プロット

plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()

Out[5]:

<matplotlib.legend.Legend at 0xabb63a2c>

signal.argrelmaxで配列データを渡すだけなら簡単だね。グラフの描き方は少し慣れてきたから大丈夫。

orderというパラメータがあったけど、1から値を変えるとどうなるの？

それならorderの値を少し増やして、結果がどうなるか実際試してみて（orderは1以上の整数を指定してね)。

In [6]:

# ピーク検出（極大値、オーダー値変更）
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

N = 128 # データ数
x = np.linspace(0, 2*np.pi, N) # x軸
y = np.sin(x) # 正弦波
y_add_noise = y + np.random.randn(N)*0.3 # 正弦波＋ノイズ（ピーク検出テスト用）

order_list = [1, 2, 4, 8]
axes_count = 4 # グラフの個数

fig = plt.figure(figsize=(12,8))

for i in range(axes_count):

    # ピーク検出
    maximal_idx = signal.argrelmax(y_add_noise, order=order_list[i]) # 極大値インデックス取得

    fig.add_subplot(2, 2, i+1)
    plt.title('order={}'.format(order_list[i]), fontsize=18)
    plt.plot(x, y,'r',label='sin(x)', c='gray', alpha=0.5) # 正弦波
    plt.plot(x, y_add_noise,label='sin(x)+noise', c='gray') # 正弦波＋ノイズ

    plt.plot(x[maximal_idx],y_add_noise[maximal_idx],'ro',label='peak_maximal') # 極大点プロット

plt.tight_layout()

orderの値を増やすと検出ピークの数が減ってる。

order値は、増減チェックを行うデータ数の幅（グラフ上のｘ軸）に相当するよ。1だと前後各1点、2だと前後各2点を対象にするみたいなイメージ。

次は極小値の取得方法。signal.argrelmin(numpy配列データ, order=1)とするればOK。使い方はargrelmaxと同じだよ。

In [7]:

# ピーク検出（極小値）
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

N = 64 # データ数
x = np.linspace(0, 2*np.pi, N) # x軸
y = np.sin(x) # 正弦波
y_add_noise = y + np.random.randn(N)*0.3 # 正弦波＋ノイズ（ピーク検出テスト用）

# ピーク検出
minimal_idx = signal.argrelmin(y_add_noise, order=1) # 極小値インデックス取得

plt.plot(x, y,'r',label='sin(x)', c='gray', alpha=0.5) # 正弦波
plt.plot(x, y_add_noise,label='sin(x)+noise', c='g') # 正弦波＋ノイズ

plt.plot(x[minimal_idx],y_add_noise[minimal_idx],'bo',label='peak_minimal') # 極小点プロット

plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()

Out[7]:

<matplotlib.legend.Legend at 0xab92882c>

極小値も簡単に取得できた。

極大値と極小値をいっぺんにプロットしてみる。

In [8]:

# ピーク検出（極大値、極小値）
from scipy import signal
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

N = 64 # データ数
x = np.linspace(0, 2*np.pi, N) # x軸
y = np.sin(x) # 正弦波
y_add_noise = y + np.random.randn(N)*0.3 # 正弦波＋ノイズ（ピーク検出テスト用）

#ピーク検出
maximal_idx = signal.argrelmax(y_add_noise, order=1) # 極大値インデックス取得
minimal_idx = signal.argrelmin(y_add_noise, order=1) # 極小値インデックス取得

plt.plot(x, y,'r',label='sin(x)', c='gray', alpha=0.5) # 正弦波
plt.plot(x, y_add_noise,label='sin(x)+noise', c='g') # 正弦波＋ノイズ

plt.plot(x[maximal_idx],y_add_noise[maximal_idx],'ro',label='peak_maximal') # 極大点プロット
plt.plot(x[minimal_idx],y_add_noise[minimal_idx],'bo',label='peak_minimal') # 極小点プロット
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()

Out[8]:

<matplotlib.legend.Legend at 0xab895c8c>

イイね。データ内容は下の例のように確認できるよ。

In [12]:

# 極大値情報の表示例
print('idx_length:', len(maximal_idx)) # peakの検出数
print('idx_value:', maximal_idx) # peakのindex
print('x_value:', x[maximal_idx]) # peakのx値
print('y_value:', y_add_noise[maximal_idx]) #peakのy値

idx_length: 1
idx_value: (array([ 3,  6,  9, 11, 14, 16, 18, 22, 24, 26, 28, 36, 43, 47, 49, 51, 55,
       57, 61], dtype=int32),)
x_value: [ 0.2991993   0.5983986   0.8975979   1.0970641   1.3962634   1.5957296
  1.7951958   2.1941282   2.3935944   2.5930606   2.7925268   3.5903916
  4.2885233   4.68745571  4.88692191  5.08638811  5.48532051  5.68478671
  6.08371911]
y_value: [ 0.96702313  0.84834658  0.90501103  1.32615386  1.11796672  1.4479129
  1.06884822  1.07126887  1.36109912  0.53516276  0.80185391 -0.06479653
 -0.32627323 -0.76644208 -1.04862984 -0.814803   -0.58736729 -0.47257849
 -0.30601209]

（最初の方に載せた、）FFT解析結果のピーク検出例は、極大値を一旦全部取得。そのあとで一定の振幅値を超えた極大値をピークとして採用しているよ。

あとは、グラフ上のピーク付近に周波数の値を表示すれば、簡易の自動ピーク検出は完成。引き出し線とテキスト挿入はMatplotlibのannotate()を利用しているよ。

検出条件も変更して練習してみる。

FFT簡易ピーク検出の練習はこれでおしまい。
お疲れ様でした。

お疲れ様でした。

またね！

〜　モモノキ＆ナノネと学習シリーズ（高速フーリエ変換）　〜

Pythonで高速フーリエ変換（FFT）の練習-5 周波数ピークを自動で検出

モモノキ＆ナノネと一緒にPythonでFFTの使い方を覚えよう（5）

極大値と極小値の取得方法を覚えてピークの自動検出に挑戦しよう

1 件のコメント :

人気記事

カテゴリー

このブログを検索

ブログアーカイブ

Pythonで高速フーリエ変換（FFT）の練習-5 周波数ピークを自動で検出

モモノキ＆ナノネと一緒にPythonでFFTの使い方を覚えよう（5）

極大値と極小値の取得方法を覚えてピークの自動検出に挑戦しよう

この記事をシェアする

1 件のコメント :

人気記事

カテゴリー

このブログを検索

ブログアーカイブ