Complex Demodulationとは？

# 初めに
本記事では，時間周波数領域の解析などに使われるComplex Demodulation法（以下，CD法と表記）についてなんとなく解説します．CD法は，**時系列データを対象に，ある周波数に対する振幅と位相の時間的変化を抽出する**手法です．

:::note warn
注意
筆者は，CD法についてなんとなく理解した程度の人間です．
内容が間違っている可能性もありますので，ご了承ください．
また，間違っていたらコメント欄にてお知らせ頂けると幸いです．
:::


# 数学的解釈
CD法をなんとなく理解するための数学的解釈です．
まず，簡単のために周波数$f$をもち振幅$A$と位相$\theta$が時間の関数で与えられる
次のような時系列データ（信号）$P(t)$を考えてみましょう．
``` math
P(t)=A(t)\sin\{2{\pi}ft+{\theta}(t)\}
```
この信号の，周波数$f$が持つ振幅と位相を見てましょう．
（もちろん，この信号には周波数$f$以外の信号は含まれていませんよね？）

$P(t)$にそれぞれ，$\sin{2{\pi}ft}$，$\cos{2{\pi}ft}$を乗算すると，
``` math
P(t)\sin{2{\pi}ft}=\frac{A}{2}\{\cos\theta-\cos(4{\pi}ft+\theta)\} 
```
``` math
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$の波**だけを取り出してみましょう．
``` math
x(t)=F\{p(t)\sin2{\pi}ft\}=\frac{A}{2}\cos\theta
```
``` math
y(t)=F\{p(t)\cos2{\pi}ft\}=\frac{A}{2}\sin\theta
```
となります．
この得られた2つの式から
``` math
{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))
```
``` math
\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$を掛けているだけではないか．

これに対する答えは，複素関数と[オイラーの公式](https://ja.wikipedia.org/wiki/%E3%82%AA%E3%82%A4%E3%83%A9%E3%83%BC%E3%81%AE%E5%85%AC%E5%BC%8F)について考えるとすぐにわかります．
オイラーの公式は，
``` math
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としています．

``` python
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
    
```

:::note warn
注意

フィルターの設定によって，抽出された振幅の大きさが変化しますので，最適なフィルターを設定する必要があります．
今回のフィルターは不十分な可能性があります．
:::

実際に動作させてみましょう．
``` python
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$は
![Figure_1.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/709083/bd335ef7-23fa-633b-bca0-4306f34ee70f.jpeg)
こういう波形を描きます．

この波形の10Hzの振幅を表示すると，
![Figure_2.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/709083/bd48c914-fdee-8003-1b8b-c71bb698435b.jpeg)
よって，振幅が1に収束していくことがわかります．ただし，0.2秒までの部分が徐々に変化しています．
これはおそらくFIRフィルターによるものです．完璧に0秒から1にするのは難しいので，最初の方の値は無視するなどの対策が必要となりそうです．