Python: Matplotlib の cohere 関数でコヒーレンスを推定

以前 に Matplotlib の psd, csd 関数で簡単に周波数応答が求められることを紹介しました。 Matplotlib にはコヒーレンスを推定する cohere 関数もあるので、簡単にコヒーレンスも求められます。

実際に周波数応答を見る場合にはコヒーレンスも求め、信号にノイズ等の影響が生じていないかを判定します。 詳しくは 小野測器さん のページ等で説明されています。

cohere(x, y)でコヒーレンス \(C_{xy}\) を推定します。関数の内部では Welch 法で各信号の PSD, CSD を計算しています。

import numpy as np
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)

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

Cyx, f = mlab.cohere(y, x, NFFT=1024, Fs=fs,
                  window=mlab.window_hanning, noverlap=1024/2)

plt.figure()
# plt.plot(f, 10*np.log10(np.abs(FRF)))
plt.plot(f, Cyx)
plt.xlim(0, fs/2)
plt.xlabel("Frequency[Hz]")
15071001.png

コメント

Comments powered by Disqus
書籍更新情報
2017-02-18
Pythonによる科学技術計算 基礎編
1.4版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
電子書籍
線形代数(1): Pythonによる科学技術計算 実践編