Python : Matplotlib の psd, csd 関数と周波数応答

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

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

関数の種類

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

Table 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()

15062001.png

各データを単純に 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()

15062002.png

さて、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()

15062003.png

コメント

Comments powered by Disqus
書籍更新情報
2016-10-21
Pythonによる科学技術計算 基礎編
PDF版の販売を開始しました。
販売ページはこちら

2016-09-09
Pythonによる科学技術計算 基礎編
1.2版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
Kindle ストア、Leanpubで販売中です
Pythonによる科学技術計算 基礎編
PDF版の販売はこちら
同人誌
技術書典(2016/6/25)
Emacs/org-modeのPDF作成術
電子版をBOOTHで販売中です
Emacs/org-modeのPDF作成術
Share