畳み込み積分の確認

畳み込み積分がいまいちイメージしづらかったので実際にプログラムを書いて確認してみた。

畳み込み積分の式

伝達関数

1自由度振動系にt=0でインパルス入力した応答を伝達関数として使う。
インパルス入力は周波数領域で全周波数で大きさ1のスペクトルを持つので、インパルス応答の周波数スペクトルがそのまま伝達関数となる性質を使う。
減衰しながら振動が収束する系。

入力

t<1までは0で1<=tで1となるステップ入力。インパルス入力がずっと続く入力とみなせる。

出力(畳み込み積分)

畳み込み積分によって最初の図のような応答を示す系に入力を入れた際の時系列の応答が求まった。ステップ入力を入れた瞬間オーバーシュートして、振動しながら一定値に収束するのでイメージ通り。

例えばt=2.0sec時点の出力はt=2.0での入力と、t=1.0からステップ入力を入れた出力の名残りの影響を受けていて、それはt=2.0secまでの畳み込み積分で求まる。
イメージとしては入力波形の時間軸を左右ひっくり返して、応答と入力の内積をとる感じが個人的にしっくりくる。

使ったコードはこちら

dt = 0.05
time = np.arange(0, 5, dt)
m = 400  # [kg]
omega = 2 * np.pi * 1.5
zeta = 0.2
sigma = omega * zeta

# -----　粘性1自由度振動系 インパルス入力の応答
h = np.exp(-sigma * time) * np.sin(omega * time) / (m * omega)

# ----- 0<=t<1まで0で 1<=t<=5まで1が続くステップ入力(インパルスが連続している状態)
f = np.hstack([np.zeros(int(1 /dt)), np.ones(int((5 -1 + dt) / dt))])

# ----- ステップ入力と粘性1自由度振動系の畳み込み積分
y = np.zeros(time.shape[0])
for t in range(time.shape[0]):
    for tau in range(t+1):
        y[t] = y[t] + f[tau] * h[t - tau]

# ----- グラフ化
fig, ax = plt.subplots(3, 1, gridspec_kw=dict(hspace=0.4), figsize=(5, 6))
ax[0].plot(time, h)
ax[0].set_xlabel("Time[sec]")
ax[1].plot(time, f)
ax[1].set_xlabel("Time[sec]")
ax[2].plot(time, y)
ax[2].set_xlabel("Time[sec]")