UENEのノート

最近やったことを書いたりするやつ。Qiitaも書いてる。

PythonによるFFTを用いたパワースペクトル推定

目次

  1. 概観
  2. 準備
  3. 実装
  4. 結果
  5. 解説
  6. 補足
  7. 参考

概観

はじめに

パワースペクトル(もしくはパワースペクトル密度)の推定方法は様々です。大きく分けて二つの手法が存在します。

論文やライブラリではパワースペクトル密度を推定することがほとんどだと思います。パワースペクトルパワースペクトル密度の違いについては別記事で解説する予定です。ちなみにパワースペクトル密度からパワースペクトルへの変換は比較的簡単に行うことが出来ます。

  1. パラメトリック手法
  2. ノンパラメトリック手法

セミパラメトリックという手法も存在するようですが、私には理解できませんでした。

私は統計学をあまり得意ではないので、それぞれの用語の正しい意味は詳しい人に聞いてみてください。今回はパワースペクトルパラメトリック手法の中のピリオドグラム法を用います。

ピリオドグラム法

ピリオドグラム法とはFFTを用いてパワースペクトルを算出する方法です。パーシバルの定理から時間信号のパワー(二乗平均値)と、振幅スペクトルの二乗値の合計値が等しいことを考えてパワースペクトルとしています。

複素フーリエ係数とパワースペクトルの関係(パーシバルの定理)

一般に信号の強さ(パワー)は時間信号の二乗平均値を表します。ここで、有限長かつ離散信号に関して時間信号は複素フーリエ係数$c_k$を利用して以下のように定義することができます。

\begin{equation} x(t) = \sum_{k=-\infty}^{\infty} c_k e^{j 2\pi k t /T} \end{equation}

二乗平均値$P$は

\begin{equation} P = \frac{1}{T} \int_{0}^{T} | x(t) |^2 dt \end{equation}

と表すことができるので、

\begin{align} P & = \frac{1}{T} \int_{0}^{T} | \sum_{k=-\infty}^{\infty} c_k e^{j 2\pi k t /T} |^2 dt \\ & = \frac{1}{T} \sum_{k=-\infty}^{\infty} | c_k |^2 \int_{0}^{T} { e^{j 2\pi k t /T} }^2 dt \\ & = \frac{1}{T} \sum_{k=-\infty}^{\infty} | c_k |^2 \int_{0}^{T} { \cos{\frac{2\pi k t}{T}} + j\sin{\frac{2\pi k t}{T}} }^2 dt \\ & = \frac{1}{T} \sum_{k=-\infty}^{\infty} | c_k |^2 \int_{0}^{T} 1 + 2j\cos{\frac{2\pi k t}{T}}\sin{\frac{2\pi k t}{T}} dt \\ & = \frac{1}{T} \sum_{k=-\infty}^{\infty} | c_k |^2 \left[ t + j\sin{\frac{4\pi k t}{N}} \right]_{0}^{T} \\ & = \sum_{k=-\infty}^{\infty} | c_k |^2 \\ & = \sum_{k=-\infty}^{-1} | c_k |^2 + \sum_{k=1}^{\infty} | c_k |^2 + | c_0 |^2 \\ & = 2\sum_{k=1}^{\infty} | c_k |^2 + | c_0 |^2 \\ & = \sum_{k=0}^{\infty} \mathrm{PowerSpectra}(k) \end{align}

ここで、 \mathrm{PowerSpectra}(k)は以下のように定義しました。

\begin{equation} \mathrm{PowerSpectra}(k) = \begin{cases} 2| c_k |^2 & (k>0)\\ | c_0 |^2 & (k=0) \end{cases} \end{equation}

複素フーリエ係数はDFTによって求めることができるので、DFTからパワースペクトルを算出することができます。ちなみに上の等式をパーシバルの定理といいます。scipy - Fourier Transformsscipyで実装されているDFTの詳細が記されています。それらに注意して、ソースコード中では以下のようにパワースペクトルを算出しています。ちなみに要素数は基本的に偶数個です。なぜならサンプリング周波数は基本的にkHzオーダーだからです。

For N even, the elements $y[1] \dots y[N/2-1]$ contain the positive-frequency terms, and the elements $y[N/2] \dots y[N-1]$ contain the negative-frequency terms, in order of decreasingly negative frequency.

power = np.hstack( \
        ( \
             np.abs( fft_signal[0] /(T*fs) )**2.0, \
        2.0* np.abs( fft_signal[1:int(T*fs/2-1.0)] /(T*fs) )**2.0 \
        ) \
        )

準備

pipを利用してインストールします。仮想環境を作っておいた方がなにかと良いですよね?!

  1. numpy

    python3(or py) -m pip install numpy

  2. scipy

    python3(or py) -m pip install scipy

  3. matplotlib

    python3(or py) -m pip install matplotlib

インストールする方法については、pypa - Installing Packagesを参考にしてください!

実装

import numpy as np
from scipy.fft import fft, fftfreq
import matplotlib.pyplot as plt

# Defining Constants
T = 100.0
fs = 192.0e3

amp1=1.0
freq1 = 5.0e3
pha1=-2*np.pi/3

amp2=np.sqrt(2.0)
freq2 = 16.0e3
pha2=2*np.pi/3

# Generating Signal
t = np.linspace(0.0, T, int(T*fs))
x = amp1* np.cos(2.0* np.pi* freq1* t+ pha1) \
        + amp2* np.cos(2.0* np.pi* freq2* t+ pha2)

# Doing FFT
fft_signal = fft(x)

# Generating Frequency Vector
f = fftfreq(int(T*fs), 1.0/fs)[0:int(T*fs/2-1.0)]/1.0e3

# Generating Amplitude & Phase
power = np.hstack( \
        ( \
             np.abs( fft_signal[0] /(T*fs) )**2.0, \
        2.0* np.abs( fft_signal[1:int(T*fs/2-1.0)] /(T*fs) )**2.0 \
        ) \
        )
power_dB = 10.0*np.log10(power/(20.0e-6)**2.0)

# Plotting
size_f=20

fig = plt.figure(figsize=[10,7], dpi=300, linewidth=0.0, tight_layout=T)

ax0 = fig.subplots(1, 1)

ax0.plot(f, power_dB)
ax0.set_xlabel('frequency[kHz]', fontsize=size_f)
ax0.set_ylabel('power spectra[dB SPL]', fontsize=size_f)
ax0.set_xlim(f[0], f[-1])
ax0.set_ylim(-70, 120)

fig.align_ylabels()
fig.suptitle('Power Spectra', fontsize=size_f)
fig.savefig('power_spectra.png')

結果

power_spectra.png

解説

dB SPLについて

株式会社 小野測器 - 騒音計とはでは

「音圧レベルはSound Pressure Level(SPL)なので、特に音圧レベルの単位である dB(デシベル)を明示的に表現するために、以前は“dB SPL” と表記する場合がありました。 しかし、最近の JIS 規格や計量法などの法律なども、単に dB と表記するのが正しい単位表記となっています。」

とあり、dB SPLという表記は正しい表記ではないそうですが、基準音圧を  20 \mathrm{\mu Pa} としたことを強調するために本記事ではdB SPLと表記することにします。

結果について

 1 Pa \simeq 93.98 dB SPL ですが、 16 kHz の部分で 93.88 dB SPL であったので、ある程度正しく算出できていると考えています。

参考