前回は SciPy のパワースペクトル密度(PSD)を推定する関数についてまとめましたが、Matplotlib にも同様の関数があるうえに、クロススペクトル密度(CSD)を求める関数まであったので紹介します。

[2015-07-10 金] Matlplotlib の cohere 関数でコヒーレンスも簡単に推定できます。

関数の種類

Matplotlib の表 1 の関数で PSD 推定、CSD 推定ができます。

表1: Matplotlib の psd, csd 関数
関数	数値解法
psd	パワースペクトル密度を推定(Welch 法)
csd	クロススペクトル密度を推定(Welch 法)

psd は SciPy の welch に相当しますが、若干オプション引数の指定が異なります。 csd は csd(x, y)で \(P_{xy}\) を推定します。

使い方の例

例として適当なデータ x とそれを FIR フィルタで処理したデータ y 間の周波数応答を計算してみます。

まず、SciPy の freqz で FIR フィルタの周波数応答を計算すると以下のようになります。

# coding: utf-8
import numpy as np
from scipy import fftpack
from scipy import signal
import matplotlib.pyplot as plt
from matplotlib import mlab

n = 1024*4
fs = 512
dt = 1/fs

b = signal.firwin(30, 0.2, window="boxcar")
x = np.random.randn(n)
y = signal.filtfilt(b, 1, x)

w, h = signal.freqz(b, 1, fs)

plt.plot(w*fs/(2*np.pi), 20*np.log10(np.abs(h)), "b")
plt.xlim(0, fs/2)
plt.ylim(-120, 0)
plt.xlabel("Frequency[Hz]")
plt.ylabel("Amplitude[dB]")
# plt.show()

各データを単純に FFT 処理し、周波数応答を求めてみると以下のようになります。ここでは FFT 処理に平均処理(オーバラップ処理)をしていないため、結果を見るとピークが鋭く立っていますが全体的に微妙に滑らかでないのがわかります。

w = signal.hanning(n)
yfIn = fftpack.fft(x*w)
yfOut = fftpack.fft(y*w)
freq = fftpack.fftfreq(n, dt)

FRF = yfOut/yfIn
# FRF = (yfOut*np.conj(yfIn))/(yfIn*np.conj(yfIn))

plt.plot(freq[1:n/2], 10*np.log10(np.abs(FRF[1:n/2])), "r")
plt.xlim(0, fs/2)
plt.ylim(-120, 0)
plt.xlabel("Frequency[Hz]")
plt.ylabel("Amplitude[dB]")
# plt.show()

さて、psd, csd を使って周波数応答を求めると以下のようになります。ここではデータ長の 1/4 をセグメント長とし、50% のオーバーラップで平均処理しています。そのため、若干ピークが鈍りますが、全体的には平滑化し、誤差が小さくなっています。どの程度平均化処理するかは実際に計測データから周波数応答を計算してみて決めればいいと思います。

Pyx, f = mlab.csd(y, x, NFFT=1024, Fs=fs,
                  window=mlab.window_hanning, noverlap=1024/2)
Pxx, f = mlab.psd(x, NFFT=1024, Fs=fs,
                  window=mlab.window_hanning, noverlap=1024/2)
FRF = Pyx/Pxx

plt.plot(f, 10*np.log10(np.abs(FRF)), "g")
plt.xlim(0, fs/2)
plt.ylim(-120, 0)
plt.xlabel("Frequency[Hz]")
plt.ylabel("Amplitude[dB]")
# plt.show()