初めに
本記事では，時間周波数領域の解析などに使われるComplex Demodulation法（以下，CD法と表記）についてなんとなく解説します．CD法は，時系列データを対象に，ある周波数に対する振幅と位相の時間的変化を抽出する手法です．
注意
筆者は，CD法についてなんとなく理解した程度の人間です．
内容が間違っている可能性もありますので，ご了承ください．
また，間違っていたらコメント欄にてお知らせ頂けると幸いです．
数学的解釈
CD法をなんとなく理解するための数学的解釈です．
まず，簡単のために周波数$f$をもち振幅$A$と位相$\theta$が時間の関数で与えられる
次のような時系列データ（信号）$P(t)$を考えてみましょう．
P(t)=A(t)\sin\{2{\pi}ft+{\theta}(t)\}
この信号の，周波数$f$が持つ振幅と位相を見てましょう．
（もちろん，この信号には周波数$f$以外の信号は含まれていませんよね？）
$P(t)$にそれぞれ，$\sin{2{\pi}ft}$，$\cos{2{\pi}ft}$を乗算すると，
P(t)\sin{2{\pi}ft}=\frac{A}{2}\{\cos\theta-\cos(4{\pi}ft+\theta)\}
P(t)\cos{2{\pi}ft}=\frac{A}{2}\{\sin\theta-\sin(4{\pi}ft+\theta)\}
と書き表せます．
上式は，$\cos\theta$や$\sin\theta$という周波数$0$の波と
$\cos(4{\pi}ft+\theta)$や$\sin(4{\pi}ft+\theta)$という周波数$2f$の波に分解できたことを表します．
そして，この2つの式にlow-passフィルター処理$F$を施して，周波数$0$の波だけを取り出してみましょう．
x(t)=F\{p(t)\sin2{\pi}ft\}=\frac{A}{2}\cos\theta
y(t)=F\{p(t)\cos2{\pi}ft\}=\frac{A}{2}\sin\theta
となります．
この得られた2つの式から
{A(t)}^2=4(\frac{A^2}{4}{\cos}^2\theta+\frac{A^2}{4}{\sin}^2\theta)=4(x^2(t)+y^2(t))
\theta(t)={\tan}^{-1}\frac{y(t)}{x(t)}
と計算できます．これにて，振幅$A$と位相$\theta$を求めることができました．
まとめると，
対象の信号に対して所望の周波数$f$の$\sin{2{\pi}ft}$，$\cos{2{\pi}ft}$を掛けて，
low-passフィルター処理を行い，2乗和と$\tan^{-1}$を取ると，
周波数$f$の振幅$A(t)$と位相$\theta(t)$を得ることができます．
これがCD法の原理となります．
なぜComplex Demodulationと呼ばれるか？
ここまで読んできた読者の人の中には，こう思う人もいるでしょう．
なぜComplex Demodulation（直訳：複素変調）と呼ぶのか．
単純に$\sin$と$\cos$を掛けているだけではないか．
これに対する答えは，複素関数とオイラーの公式について考えるとすぐにわかります．
オイラーの公式は，
e^{ix}={\cos}x+i{\sin}x
という式です．
オイラーの公式を利用することで，今まで$\sin$と$\cos$をべつべつに乗算して計算していたところを，
$e^{ix}$を1回乗算するだけでこの計算を終わらせることができます．
乗算した後は，実部と虚部の部分を取り出して2乗和と$\tan^{-1}$を計算するだけです．
これが，Complex Demodulation（直訳：複素変調）と呼ばれる所以だと思います．
実際は，プログラムにするとき複素数を扱うとわかりずらくなる（ような感じ）がするので，
最初説明した方法で行います．
（計算負荷的な観点から考えると，複素数を使った方が1回で済むので早いかもしれませんが．）
Pythonでの実装
それでは，Pythonにてクラスの実装を行いましょう．
ライブラリには，numpyとscipyを使っています．
今回，low-passフィルターとしてscipyのFIRフィルターを使っています．
タップ数は255，カットオフ周波数は4Hzとしています．
import numpy as np
from scipy import signal
class time_frequency_analysis:
__fs = 0
__time = []
__data = []
def __init__(self, sampling_freq: int, time: np.ndarray, data: np.ndarray):
self.__fs = sampling_freq
self.__time = time
self.__data = data
def complex_demodulation(self, freq: int):
# amp,phase:numpy.ndarray
freq_time = 2 * np.pi * freq * np.array(self.__time)
cos_array = np.cos(freq_time)
sin_array = np.sin(freq_time)
temp_array_cos = cos_array * np.array(self.__data)
temp_array_sin = sin_array * np.array(self.__data)
# FIRフィルタ
fir_filter = signal.firwin(numtaps=255, fs=self.__fs, cutoff=4, pass_zero=True)
filtered_temp_array_cos = signal.lfilter(fir_filter, 1, temp_array_cos)
filtered_temp_array_sin = signal.lfilter(fir_filter, 1, temp_array_sin)
amp = np.sqrt(4*filtered_temp_array_cos**2 + 4*filtered_temp_array_sin**2)
phase = np.arctan2(filtered_temp_array_sin, filtered_temp_array_cos)
return amp, phase
注意
フィルターの設定によって，抽出された振幅の大きさが変化しますので，最適なフィルターを設定する必要があります．
今回のフィルターは不十分な可能性があります．
実際に動作させてみましょう．
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
class time_frequency_analysis:
__fs = 0
__time = []
__data = []
def __init__(self, sampling_freq: int, time: np.ndarray, data: np.ndarray):
self.__fs = sampling_freq
self.__time = time
self.__data = data
def complex_demodulation(self, freq: int):
# amp,phase:numpy.ndarray
freq_time = 2 * np.pi * freq * np.array(self.__time)
cos_array = np.cos(freq_time)
sin_array = np.sin(freq_time)
temp_array_cos = cos_array * np.array(self.__data)
temp_array_sin = sin_array * np.array(self.__data)
# FIRフィルタ
fir_filter = signal.firwin(numtaps=255, fs=self.__fs, cutoff=4, pass_zero=True)
filtered_temp_array_cos = signal.lfilter(fir_filter, 1, temp_array_cos)
filtered_temp_array_sin = signal.lfilter(fir_filter, 1, temp_array_sin)
amp = np.sqrt(
4 * filtered_temp_array_cos**2 + 4 * filtered_temp_array_sin**2
)
phase = np.arctan2(filtered_temp_array_sin, filtered_temp_array_cos)
return amp, phase
if __name__ == "__main__":
f_test1 = 10 # 波形の周波数
fs = 1000 # サンプリング周波数
N = 50000 # サンプリング点数
t = np.arange(N) / fs # 時間
y = []
for i in t:
y.append(math.sin(2 * (math.pi) * f_test1 * i))
# 波形の出力
plt.figure()
plt.plot(t, y)
plt.legend()
plt.xlabel("t[s]")
plt.ylabel("y[m]")
plt.xlim(0, 1)
plt.grid()
plt.show()
CD = time_frequency_analysis(fs, t, y)
amp, phase = CD.complex_demodulation(10)
plt.figure()
plt.plot(t, amp)
plt.legend()
plt.xlabel("t[s]")
plt.ylabel("amp")
plt.xlim(0, 1)
plt.grid()
plt.show()
実行すると，解析対象の信号$\sin20{\pi}t$は
こういう波形を描きます．
この波形の10Hzの振幅を表示すると，
よって，振幅が1に収束していくことがわかります．ただし，0.2秒までの部分が徐々に変化しています．
これはおそらくFIRフィルターによるものです．完璧に0秒から1にするのは難しいので，最初の方の値は無視するなどの対策が必要となりそうです．