はじめに
信号処理において、相互スペクトル密度(Cross Spectral Density, CSD)は、周波数領域における 2 つの信号間の相関を測定する指標です。これは、2 つの信号が周波数成分に関してどれだけ関連しているかを判断するために使用されます。この実験では、Python の Matplotlib ライブラリを使用して 2 つの信号の CSD を計算する方法を学びます。
VM のヒント
VM の起動が完了したら、左上隅をクリックしてNotebookタブに切り替え、Jupyter Notebook を開いて練習を行ってください。
時々、Jupyter Notebook の読み込みが完了するまで数秒待つ必要がある場合があります。Jupyter Notebook の制限により、操作の検証を自動化することはできません。
学習中に問題が発生した場合は、いつでも Labby に質問してください。セッション終了後にフィードバックを提供していただければ、迅速に問題を解決します。
必要なライブラリをインポートする
以下のライブラリをインポートする必要があります。numpy と matplotlib.pyplot です。
import matplotlib.pyplot as plt
import numpy as np
信号を生成する
2 つの信号を生成する必要があります。これらの信号は、コヒーレントな部分とランダムな部分を含んでいます。両方の信号のコヒーレントな部分の周波数は 10 Hz です。信号のランダムな部分は、白色雑音をローパスフィルタに通して有色雑音を生成することで作成されます。
dt = 0.01
t = np.arange(0, 30, dt)
## Fixing random state for reproducibility
np.random.seed(19680801)
nse1 = np.random.randn(len(t)) ## white noise 1
nse2 = np.random.randn(len(t)) ## white noise 2
r = np.exp(-t / 0.05)
cnse1 = np.convolve(nse1, r, mode='same') * dt ## colored noise 1
cnse2 = np.convolve(nse2, r, mode='same') * dt ## colored noise 2
## two signals with a coherent part and a random part
s1 = 0.01 * np.sin(2 * np.pi * 10 * t) + cnse1
s2 = 0.01 * np.sin(2 * np.pi * 10 * t) + cnse2
信号をプロットする
生成した 2 つの信号を、Matplotlib の plot 関数を使ってプロットすることができます。
fig, ax = plt.subplots()
ax.plot(t, s1, label='s1')
ax.plot(t, s2, label='s2')
ax.set_xlabel('Time')
ax.set_ylabel('Amplitude')
ax.legend()
plt.show()
CSD を計算する
2 つの信号の相互スペクトル密度(Cross Spectral Density, CSD)を計算するには、Matplotlib の csd 関数を使用する必要があります。この関数は、2 つの信号、FFT(高速フーリエ変換)のポイント数、およびサンプリング周波数を入力として受け取ります。
fig, ax = plt.subplots()
cxy, f = ax.csd(s1, s2, 256, 1. / dt)
ax.set_ylabel('CSD (dB)')
plt.show()
結果を解釈する
得られたプロットは、2 つの信号の相互スペクトル密度(CSD)を示しています。x 軸は周波数を表し、y 軸はその周波数における 2 つの信号間の相関の強さを表しています。プロットは 10 Hz にピークを示しており、これは信号のコヒーレントな部分の周波数です。これは、2 つの信号がこの周波数で強く相関していることを示しています。
まとめ
この実験では、Python の Matplotlib ライブラリを使用して 2 つの信号の相互スペクトル密度(Cross Spectral Density, CSD)を計算する方法を学びました。コヒーレントな部分とランダムな部分を持つ 2 つの信号を生成し、信号をプロットし、CSD を計算し、結果を解釈しました。CSD は、周波数領域における 2 つの信号間の相関を判断するための有用なツールです。