6
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

matplotlibとscipyのフーリエ変換(FFT)を用いてパーシバルの定理を確認する方法

Last updated at Posted at 2020-03-31

python でパワースペクトルを計算する方法

matplotlib と scipy で FFT を計算する方法

matplotlib と scipy のどちらもパワースペクトルを計算する方法があるが、両者でデフォルトで持っている関数と定義が違う。大概の使い分けは次の通り。

  • 単純にパワースペクトルを計算する場合、mlab.psd が一番簡単。
  • scipy.fftpack を使うと実部と虚部が得られる。位相の計算など細かいことが可能だか、規格化は自分で行う必要あり。

どちらでやっても、パーセバルの定理から全パワーは同じになるが、両者の規格化の違いを考えないとパワースペクトルの縦軸が完璧に同じにならない。

ここでは単純な三角関数を両方の方法でフーリエ変換して、完璧にmatplotlibとscipyのFFTで結果が同じになる例を示す。おまけで、パワースペクトルから三角波の振幅を求めて、それがもとの三角関数の振幅と同じになること、フィルターを入れるとどういう影響があるかを紹介する。

問題設定

フーリエ変換する前の関数

縦軸は特に意味はないが、電圧V(t)という信号で単一の三角関数を考える。最も簡単な表式として$$y=V(t)=\sin(t)$$を用いる。周期 $f$ の三角関数は、$$\sin(2\pi ft)$$と書けるので、周期 $$f=1/2\pi$$ という設定になる。振幅は1.0である。

これを mlab.psdscipy.fftpackの両方でフーリエ変換してパワーを計算する。それぞれ、hanning filter をかけた場合も計算し、違いを見てみる。

時間分解能 dt は0.1秒とする。ナイキスト周波数 f_n は、その逆数の半分で、5Hzとなる。

サンプルコード

qiita_fft_matplotlib_scipy_comp.py
#!/usr/bin/env python

"""
#2013-03-12 ; ver 1.0; First version, using Sawada-kun's code
#2013-08-28 ; ver 1.1; refactoring
#2020-03-31 ; ver 1.2; updated for python3
"""
__version__=  '1.2'

import numpy as np
import scipy as sp
import scipy.fftpack as sf
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import optparse

def scipy_fft(inputarray,dt,filtname=None, pltfrag=False):

    bin = len(inputarray)
    if filtname == None:
        filt = np.ones(len(inputarray))
    elif filtname == 'hanning':
#        filt = sp.hanning(len(inputarray))
        filt = np.hanning(len(inputarray))                
    else:
        print('No support for filtname=%s' % filtname)
        exit()

    freq = sf.fftfreq(bin,dt)[0:int(bin/2)]
    df = freq[1]-freq[0]
    fft = sf.fft(inputarray*filt,bin)[0:int(bin/2)]

    real = fft.real
    imag = fft.imag
    psd = np.abs(fft)/np.sqrt(df)/(bin/np.sqrt(2.))

    if pltfrag: 
        binarray = range(0,bin,1)
        F = plt.figure(figsize=(10,8.))
        plt.subplots_adjust(wspace = 0.1, hspace = 0.3, top=0.9, bottom = 0.08, right=0.92, left = 0.1)
        ax = plt.subplot(1,1,1)
        plt.title("check fft of scipy")
        plt.xscale('linear')
        plt.grid(True)
        plt.xlabel(r'number of bins')
        plt.ylabel('FFT input')
        plt.errorbar(binarray, inputarray, fmt='r', label="raw data")
        plt.errorbar(binarray, filt, fmt='b--', label="filter")
        plt.errorbar(binarray, inputarray * filt, fmt='g', label="filtered raw data")
        plt.legend(numpoints=1, frameon=False, loc=1)
        plt.savefig("scipy_rawdata_beforefft.png")
        plt.show()

    return(freq,real,imag,psd)


usage = u'%prog [-t 100] [-d TRUE]'
version = __version__
parser = optparse.OptionParser(usage=usage, version=version)
parser.add_option('-f', '--outputfilename', action='store',type='string',help='output name',metavar='OUTPUTFILENAME', default='normcheck')
parser.add_option('-d', '--debug', action='store_true', help='The flag to show detailed information', metavar='DEBUG', default=False)
parser.add_option('-t', '--timelength', action='store', type='float', help='Time length of input data', metavar='TIMELENGTH', default=20.)
options, args = parser.parse_args()
argc = len(args)

outputfilename =  options.outputfilename
debug =  options.debug
timelength =  options.timelength

timebin = 0.1 # timing resolution
# Define Sine curve
inputf = 1/ (2 * np.pi) # f = 1/2pi -> sin (t)
inputt = 1 / inputf # T = 1/f 
t = np.arange(0.0, timelength*np.pi, timebin)
c = np.sin(t)
fNum = len(c)
Nyquist = 0.5 / timebin
Fbin = Nyquist / (0.5 * fNum)
fftnorm =  1. / (np.sqrt(2) * np.sqrt(Fbin))

print("................. Time domain ................................")   
print("timebin = " + str(timebin) + " (sec) ")
print("frequency = " + str(inputf) + " (Hz) ")
print("period = " + str(inputt) + " (s) ")
print("Nyquist (Hz) = " + str(Nyquist) + " (Hz) ")
print("Frequency bin size (Hz) = " + str(Fbin) + " (Hz) ")
print("FFT norm  = " + str(fftnorm) + "\n")

# Do FFT.
# Hanning Filtered
psd2, freqlist = mlab.psd(c,
                         fNum,
                         1./timebin,
                         window=mlab.window_hanning,
                         sides='onesided',
                         scale_by_freq=True
                         )
psd = np.sqrt(psd2)    

spfreq,spreal,spimag,sppower = scipy_fft(c,timebin,'hanning',True)

# No Filter
psd2_nofil, freqlist_nofil = mlab.psd(c,
                                     fNum,
                                     1./timebin,
                                     window=mlab.window_none,
                                     sides='onesided',
                                     scale_by_freq=True
                                     )
psd_nofil = np.sqrt(psd2_nofil)    
spfreq_nf,spreal_nf,spimag_nf,sppower_nf = scipy_fft(c,timebin,None)

# Get input norm from FFT intergral
print("................. FFT results ................................")       
amp = np.sqrt(2) * np.sqrt(np.sum(psd * psd) * Fbin)
print("Amp  = " + str(amp))
peakval = amp / (np.sqrt(2) * np.sqrt(Fbin) )
print("Peakval  = " + str(peakval))

# (1) Plot Time vs. Raw value
F = plt.figure(figsize=(10,8.))
plt.subplots_adjust(wspace = 0.1, hspace = 0.3, top=0.9, bottom = 0.08, right=0.92, left = 0.1)

ax = plt.subplot(3,1,1)
plt.title("FFT Norm Check")
plt.xscale('linear')
plt.xlabel(r'Time (s)')
plt.ylabel('FFT input (V)')
plt.errorbar(t, c, fmt='r', label="raw data")
plt.errorbar(t, amp * np.ones(len(c)), fmt='b--', label="Amplitude from FFT", alpha=0.5, lw=1)
plt.ylim(-2,2)
plt.legend(numpoints=1, frameon=True, loc="best", fontsize=8)

# (2a) Plot Freq vs. Power (log-log)
ax = plt.subplot(3,1,2)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'Frequency (Hz)')
plt.ylabel(r'FFT Output (V/$\sqrt{Hz}$)')
plt.errorbar(freqlist_nofil , psd_nofil, fmt='r-', label="Window None, matplotlib", alpha=0.8)    
plt.errorbar(freqlist, psd, fmt='b-', label="Window Hanning, matplotlib", alpha=0.8)
plt.errorbar(spfreq_nf, sppower_nf, fmt='mo', label="Window None, scipy", ms = 3, alpha=0.8)    
plt.errorbar(spfreq, sppower, fmt='co', label="Window Hanning, scipy", ms = 3, alpha=0.8)
plt.errorbar(freqlist, fftnorm * np.ones(len(freqlist)), fmt='r--', label=r"1./($\sqrt{2} \sqrt{Fbin}$)", alpha=0.5, lw=1)
plt.legend(numpoints=1, frameon=True, loc="best", fontsize=8)

# (2b) Plot Freq vs. Power (lin-lin)
ax = plt.subplot(3,1,3)
plt.xscale('linear')
plt.xlim(0.1,0.2)
plt.yscale('linear')
plt.xlabel(r'Frequency (Hz)')
plt.ylabel(r'FFT Output (V/$\sqrt{Hz}$)')    
plt.errorbar(freqlist_nofil , psd_nofil, fmt='r-', label="Window None, matplotlib", alpha=0.8)    
plt.errorbar(freqlist, psd, fmt='b-', label="Window Hanning, matplotlib", alpha=0.8)
plt.errorbar(spfreq_nf, sppower_nf, fmt='mo', label="Window None, scipy", ms = 3, alpha=0.8)    
plt.errorbar(spfreq, sppower, fmt='co', label="Window Hanning, scipy", ms = 3, alpha=0.8)
plt.errorbar(freqlist, fftnorm * np.ones(len(freqlist)), fmt='r--', label=r"1./($\sqrt{2} \sqrt{Fbin}$)", alpha=0.5, lw=1)
plt.legend(numpoints=1, frameon=True, loc="best", fontsize=8)

plt.savefig("fft_norm_check.png")
plt.show()

実行方法は、qiita_fft_matplotlib_scipy_comp.py をダウンロードし、macであればターミナル、windows であれば anaconda prompt で、下記のように、単に何も指定しないで実行するか、、-t で時間幅を変えるなどして実行する。

$ python qiita_fft_matplotlib_scipy_comp.py -h
Usage: qiita_fft_matplotlib_scipy_comp.py [-t 100] [-d TRUE]

Options:
  --version             show program's version number and exit
  -h, --help            show this help message and exit
  -f OUTPUTFILENAME, --outputfilename=OUTPUTFILENAME
                        output name
  -d, --debug           The flag to show detailed information
  -t TIMELENGTH, --timelength=TIMELENGTH
                        Time length of input data

フーリエ変換する前の時系列データ

何もオプションを指定しないで実行すると、次のような時系列データのフーリエ変換を実行する。

scipy_rawdata_beforefft.png

hanningフィルタがかかる場合、元の時系列データ(赤)に対して、hanning フィルタ(青)がかかって、緑のデータをフーリエ変換することになる。フィルターを使わない場合は、赤を周期境界条件でフーリエ変換することになる。フィルターをかける場合は、オフセット(周波数=0の成分)が乗っている場合は、フィルターによる人為的な低周波ノイズがのる場合があるので、フィルターを使う場合はDC成分は引いてから使うほうがよい。

パワースペクトルの計算結果

fft_norm_check.png

(上段) フーリエ変換する時系列データ。パワースペクトルから計算した時系列データの振幅を重ねて表示した。(中段)赤がフィルターなしの場合、青がhanning filter ありの場合。実線がmatplotlibで計算した場合で、ドットはscipyで計算したもの。赤い薄い線は時系列データから推定したパワーの値。(下段)中段の縦軸リニアバージョン。

パーセバルの定理の確認

パーセバルの定理より、周波数空間と時空間の積分値のパワーは一致する。三角関数の場合は、変動のパワー(RMS)と振幅は$\sqrt{2}$で結ばれているので、振幅とRMSは相互に変換できる。ここでは、周波数空間を全積分したパワーから時系列データの振幅(この例では1)を計算し、それが一致することを上の絵の青点線で示した。

周波数空間で積分している箇所は、コードの、

amp = np.sqrt(2) * np.sqrt(np.sum(psd * psd) * Fbin)
print("Amp  = " + str(amp))

に該当する。

amplitude = \sqrt{2} \sqrt{ \sum_i psd(f_i)^2 \Delta f }

で計算している。この値がこの例では1になっていて、

plt.errorbar(t, amp * np.ones(len(c)), fmt='b--', label="Amplitude from FFT", alpha=0.5, lw=1)

で上段に青い点線でプロットしている。

hanning フィルターの効果

フィルタがなければ、mlab.psdscipy.fftpack は完璧に一致する。hanning filterを使うとパワーはある割合で下がるのでその分を補正する必要があるのだが、matplotlibとscipyで定義が違うのか(深追いしてない)、デフォルトのまま使うと両者の結果は一致はしない。

コードの補足

mlab.psdの使い方

# Do FFT.
# Hanning Filtered
psd2, freqlist = mlab.psd(c,
                         fNum,
                         1./timebin,
                         window=mlab.window_hanning,
                         sides='onesided',
                         scale_by_freq=True
                         )
psd = np.sqrt(psd2)    

ここで、c はFFTしたい同時間感覚幅で取得されたデータ。FFTはサンプルが全く同じ時間間隔でされたことを前提としている。(もし、時間間隔が変わる場合にはダイレクトに離散フーリエ変換をするが、計算時間はかかる。精密な実験などではそうする場合もある。) fNumはサンプル数、1/timebin で時間幅の逆数を指定する、window は None とすれば何もなし、mlabのオプションにいろんなフィルターがあるが、ここでは hanningの例を示した。sidesは onesided で片面(周波数が正の値)で積分したときにトータルのパワーが時系列のパワーと一致するように規格化するという意味。scale_by_freq で、周波数あたりに直す。これにより、返り値はpower^2/Hzの単位となる。ここでは最後に、psd = np.sqrt(psd2) をすることで、power/sqrt(Hz)の単位にしている。

scipy.fftpack の使い方

    freq = sf.fftfreq(bin,dt)[0:int(bin/2)]
    df = freq[1]-freq[0]
    fft = sf.fft(inputarray*filt,bin)[0:int(bin/2)]

fftfreq という関数を使うことで、周波数を計算してくれる。これは実はとっても大事な機能である。なぜかというと、ここscipyのfftでは複素フーリエ変換を実行するのであるが、返り値は実部と虚部の係数であり、それをどういう周波数の順番でユーザーに返すのかは自明ではないので、対応のとれた周波数を返してくれる関数が必要になる。

    real = fft.real
    imag = fft.imag
    psd = np.abs(fft)/np.sqrt(df)/(bin/np.sqrt(2.))

で、実部を虚部を取り出し、片側(周波数が正)で規格化して帳尻があるように規格化を行う。

6
8
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
6
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?