ディジタル信号処理｜周波数解析(FFT)をPythonで実装

これまでに整理した周波数解析や周波数応答の内容をもとに、サンプルプログラムを作って実際の使い方についてまとめます🐜

今回は周波数解析（FFT）を扱います。なお、対象とする信号は1次元の時系列（波形）データとします。

処理の大まかな流れ
開発環境
サンプルプログラム
解説
- 合成波作成から抽出まで
- FFT
  - sample_fft.pyのfft関数
  - sample_fft.pyのグラフプロットの結果
おわりに

処理の大まかな流れ

前処理を含めた処理の大まかな流れは以下の通りです。

測定した波形データを読み込む
読み込んだ波形データの一部を取り出す
FFTして周波数解析をする

ここで大事なのは、2.のデータ抽出です。

データ解析のしやすさから、適当な長さに抜き出してFFTにかけますが、ただ抜き出すだけでは不十分です。
なぜなら、FFTまたはDFTでは信号は周期的であることを想定しているため、ただ抽出してしまうと信号の両端が不連続となり、前提が崩れてしまいます。

そこで、例えばHann窓（ハン窓）やHamming窓（ハミング窓）などの窓関数を利用します。

抽出した波形データ全体に、上述の窓関数を畳み込むことで、強引に周期性をもつ連続信号に変換します。
この窓関数をかけた信号をFFTすることで、前提を崩すことなく解析することができます。

もちろん、窓関数をかけると大きな低周波が含まれることになりますので、それを頭に入れてFFTの結果を見ることになります。

開発環境

後述するサンプルプログラムを実行する開発環境は次の通りです。

項目	内容
OS	Windows 10
Python	3.7.6(Anaconda 4.8.1)
Numpy	1.17.4
Scipy	1.3.2
matplotlib	3.1.1

サンプルプログラム

Pythonを使ったFFTのサンプルプログラム(sample_fft.py)は以下の通りです。自由にコピペして、実際に動かしてみてください。

#!/usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack

class TestFFT():
    def __init__(self):
        self.msg = 'Sample script for fft'
        # 設定値
        self.fs = 250   # サンプリング周波数
        self.sec = 10   # 時間
        self.start_time = 2.25  # 抽出開始時間
        self.period = 3.5       # 抽出期間

        # 目的の信号
        self.A = 1  # 振幅
        self.f = 5  # 周波数

        # ノイズ｜sin波
        self.An = 0.1
        self.fn = 50

    def process(self):
        print(self.msg)
        # 目的の信号
        t = np.arange(0, self.sec, 1/self.fs)
        target_signal = self.A * np.sin(2.0 * np.pi * self.f * t)

        # ノイズ｜sin波
        noise_sin = self.An * np.sin(2.0 * np.pi * self.fn * t)

        # 合成波
        conbined_signal = target_signal + noise_sin

        # 波形の一部を抜き出す
        start = int(self.start_time * self.fs)  # インデックスは整数
        end = int(self.start_time + self.period) * self.fs
        y = conbined_signal[start:end]
        x = t[start:end]

        # FFT
        fft(y, self.fs)

        plt.show()


def fft(array, fs, ylabel_name='Data', xlabel_name='Time'):
    # 窓関数で波形を切り出す
    L = len(array)
    #window = np.hamming(L)  # ハミング窓
    window = np.hanning(L)  # ハン窓
    array_window = array * window

    # FFT計算
    # numpy.fftよりscipy.fftpackの方が速い
    NFFT = 2**nextpow2(L)  # 計算速度向上のため解析データ数に近い2の乗数を計算
    fft_amp = fftpack.fft(array_window, NFFT)  # 周波数領域のAmplitude
    fft_fq = fftpack.fftfreq(NFFT, d=1.0/fs)  # Amplitudeに対応する周波数
    # 正の領域のみ抽出
    fft_amp = fft_amp[0: int(len(fft_amp)/2)]
    fft_fq = fft_fq[0: int(len(fft_fq)/2)]
    fft_amp = db(abs(fft_amp))  # 複素数→デシベル変換

    # グラフ表示
    plt.figure(figsize=(8, 6*1.5))
    plt.subplots_adjust(hspace=0.4)  # x軸ラベル表示のため余白調整
    # 入力データプロット
    plt.subplot(3, 1, 1)
    plt.plot(array)
    plt.ylabel(ylabel_name)
    plt.grid()
    # 入力データ*窓関数プロット
    plt.subplot(3, 1, 2)
    plt.plot(array_window)
    plt.xlabel(xlabel_name)
    #plt.ylabel('Data*hamming')
    plt.ylabel('Data*hann')
    plt.grid()
    # FFTプロット
    plt.subplot(3, 1, 3)
    plt.plot(fft_fq, fft_amp)
    plt.xlim(0, fs/2)
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dB]')  # |X(omega)|
    plt.grid()

# MATLAB関数nextpow2
def nextpow2(n):
    m_f = np.log2(n)
    m_i = np.ceil(m_f)
    return int(np.log2(2**m_i))

# デシベル変換
def db(x, dBref=1):
    y = 20 * np.log10(x / dBref)
    return y


if __name__ == '__main__':
    proc = TestFFT()
    proc.process()

実行結果は次の通りです。

解説

合成波の条件とサンプルプログラムを解説します。

合成波作成から抽出まで

データ取得の条件は次の通りとします¹。

項目	内容
サンプリング周波数[Hz]	250
データ取得時間[s]	10

また、取得した波形は以下の2つからなる合成波としました。いわゆる、高周波ノイズが重畳した場合を考えます。

波形	波形情報
抽出したい目的の波形	振幅1, 周波数5[Hz]
除去したいノイズ波形	振幅0.1, 周波数50[Hz]

詳細な説明はしませんが、参考に連続時間のsin波と離散時間のsin波の式を以下に示します。

$\displaystyle x(t) = A \sin{ \omega t } \\ x[n] = A \sin{(2 \pi f \frac{n}{f_{s}} ) }$

$\omega$ は角周波数、 $f$ はサンプリング周波数、 $f_{s}$ は信号の周波数です。

抽出期間は、2.25秒後から3.5秒後の1.25秒間とします。中途半端な値を使っている理由は、実際にあるような両端が不連続である例を扱いたいだけで、それ以外に意味はありません。

FFT

基本的にはこれまでの下記ブログの内容を実装していけばよいので、軽くコメントする程度にします。

ディジタル信号処理｜周波数解析(FFT)とは - ari23の研究ノート

sample_fft.pyのfft関数

上から順に解説していきます。

def fft(array, fs, ylabel_name='Data', xlabel_name='Time'):
    # 窓関数で波形を切り出す
    L = len(array)
    #window = np.hamming(L)  # ハミング窓
    window = np.hanning(L)  # ハン窓
    array_window = array * window

今回はHann窓を使用し、Hamming窓はコメントアウトしてあります。

    # FFT計算
    # numpy.fftよりscipy.fftpackの方が速い
    NFFT = 2**nextpow2(L)  # 計算速度向上のため解析データ数に近い2の乗数を計算
    fft_amp = fftpack.fft(array_window, NFFT)  # 周波数領域のAmplitude
    fft_fq = fftpack.fftfreq(NFFT, d=1.0/fs)  # Amplitudeに対応する周波数

FFTの前にあるnextpow2関数は、引数の整数に一番近い2の累乗のべき指数（要は2の右肩）を計算します。
MATLABのnextpow2関数を模したものです。これを使ってFFTを計算量を最適化します。

    # 正の領域のみ抽出
    fft_amp = fft_amp[0: int(len(fft_amp)/2)]
    fft_fq = fft_fq[0: int(len(fft_fq)/2)]
    fft_amp = db(abs(fft_amp))  # 複素数→デシベル変換

FFTした結果は線対称な正負の値なので、正の値のみ使用します。
次に絶対値→デシベル値に変換して、最後にグラフプロットします。

sample_fft.pyのグラフプロットの結果

以下に再掲し、上から順に解説します。

まず、1つ目は切り取った波形をそのまま出しています。
両端が不連続になるよう切り取られていることが確認できます。

2つ目は、切り取った波形に窓関数をかけた（畳み込んだ）波形を出しています。
今回はHann窓をかけたので、波形の両端が連続になったことがわかります。

3つ目は、窓関数をかけた波形をFFTした結果です。
窓関数によって全体が右肩下がりの波形で、5Hzと50Hzにピークがあることが確認できます。

これでFFTの周波数解析のサンプルプログラムを作ることができました。

おわりに

こちらを参考に、サンプルプログラムを作成しました。

中身を知らずに使用すると間違って実装することが多々あるので、できる限り学問的な背景を確認するとよいと思います。

参考になれば幸いです(^^)

日本語だと秒を[sec]とよく表記されるが、英語論文では[s]が一般的↩