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

はじめに

今回は、下のような信号のパワースペクトルを求めたい場合どのようにすれば求められるかについて記事にしました。

概観

パワースペクトル～数学上～

パワースペクトル～プログラム上～

複素フーリエ係数はDFTによって求めることができるので、DFTからパワースペクトルを算出することができます。scipy - Fourier Transformsにscipyで実装されているDFTの詳細が記されています。それらに注意して、ソースコード中では以下のようにパワースペクトルを算出しています。

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 \
        ) \
        )

実装

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

# Defining Constants
T = 10.0
fs = 192.0e3

amp1=1.0
freq1 = 5.0e3
pha1=-2.0*np.pi/3.0

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

# 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 Power
power = np.hstack( \
        ( \
             np.abs( fft_signal[0] /(T*fs) )**2.0, \
        2.0* np.abs( fft_signal[1:int(T*fs/2.0-1.0)] /(T*fs) )**2.0 \
        ) \
        )
power_min = (20.0e-6)**2.0

power_dB = 10.0*np.log10(power/power_min)

# Generating Frequency Vector
f = np.linspace(0.0, fs/2.0, len(power)) /1.0e3

# Plotting
size_f = 20
fig = plt.figure(figsize=[10,7], dpi=300, linewidth=0.0, tight_layout=T)
axs = fig.subplots(2, 1)

ax0 = axs[0]
ax1 = axs[1]

ax0.plot(t, x)
ax0.set_xlabel('time[sec]', fontsize=size_f)
ax0.set_ylabel('amplitude[Pa]', fontsize=size_f)
ax0.tick_params(axis='x', labelsize=size_f)
ax0.tick_params(axis='y', labelsize=size_f)
ax0.set_xlim(t[0], t[-1])
ax0.set_ylim(-3, 3)

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

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

結果

解説

dB SPLについて

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

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

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

結果について

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

参考

• SPECTRAL AUDIO SIGNAL PROCESSING - The Periodogram