伝達関数・ラプラス変換・周波数特性

関連する定理

スケールした確率変数に対する確率密度関数

連続値をとる確率変数 Xがありその確率密度関数が P_{X}(x)であったとする。このとき Xをスケールして定義される新たな確率変数 Y\equiv aXの確率密度関数 P_{Y}(y)は

\begin{eqnarray} P_Y(y) = \frac{1}{a}P_{X}(x)\ \ \ (ただしy=ax) \end{eqnarray}

である。

証明 確率変数 Xに対して区間 [x, x+dx]に収まる確率は

\begin{eqnarray} P_{X}(x)dx \end{eqnarray}

である。 Y=aXと定義される確率変数 Yに対して、区間 [ax, ax+adx]に収まる確率は

\begin{eqnarray} P_{Y}(ax)adx \end{eqnarray}

であり、この両者が等しくないといけないことから

\begin{eqnarray} P_{X}(x)dx &=& P_{Y}(ax)adx \\ P_{Y}(ax) &=& \frac{1}{a}P_{X}(x) \end{eqnarray}

が成立する。

スケールしたガウシアンの確率密度に対する定理

G(x; \mu, \sigma)を期待値が \muで標準偏差が \sigmaのガウシアンに従う分布であるとすると

\begin{eqnarray} \frac{1}{a}G\left(\frac{x}{a}; \mu, \sigma \right) = G(x; a\mu, a\sigma) \end{eqnarray}

が成立する。

証明

\begin{eqnarray} \frac{1}{a}G\left(\frac{x}{a}; \mu, \sigma \right) &=& \frac{1}{a}\frac{1}{\sqrt{2\pi \sigma^{2}}} \exp \left(-\frac{(\frac{x}{a}-\mu)^{2}}{2\sigma^{2}} \right) \\ &=& \frac{1}{\sqrt{2\pi (a\sigma)^{2}}} \exp \left( -\frac{(x-a\mu)^{2}}{2(a\sigma)^{2}} \right) \\ &=& G(x; a\mu, a\sigma) \end{eqnarray}

たたみこみ積分のフーリエ変換に関する定理

2つの関数をたたみこみしたもののフーリエ変換は、それぞれフーリエ変換したものの積となる。
例えば f(t), g(t), (f\otimes g)(t)をフーリエ変換したものをそれぞれ F(k), G(k), FG(k)とすると

\begin{eqnarray} FG(k) = F(k) \times G(k) \end{eqnarray}

が成立する。

証明

\begin{eqnarray} FG(k) &=& \int dx \left[ \int ds f(s) g(x-s) \right] e^{-ikx} \\ &=& \int dx \left[ \int ds f(s) g(x-s) \right] e^{-iks} e^{-ik(x-s)} \\ &=& \int dx \left[ \int ds f(s) g(t) \right] e^{-iks} e^{-ikt} \ \ \ \ \ (t \equiv x-s)\\ &=& \int ds f(s) e^{-iks} \int dt g(t) e^{-ikt} \\ &=& F(k) \times G(k) \end{eqnarray}

が成立する。

ガウシアンのフーリエ変換をするとガウシアンになる

期待値が0、標準偏差が \sigmaのガウシアンに従う関数をフーリエ変換すると期待値が0で標準偏差が \frac{1}{\sigma}のガウシアンになる。

証明 まずは期待値が0のガウシアンを考え、 f(x) = \frac{1}{2\pi\sigma^{2}}e^{-\frac{x^{2}}{2\sigma^{2}}}に対してフーリエ変換 F(\omega)を求める。

\begin{eqnarray} F(\omega) &=& \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi \sigma^{2}}} e^{-\frac{x^{2}}{2\sigma^{2}}} e^{-i\omega x} dx \\ &=& \frac{1}{\sqrt{\pi \eta}} \int_{-\infty}^{\infty} e^{-\eta x^{2}-i\omega x} dx \ \ \ \ \ (\eta \equiv \frac{1}{2\sigma^{2}})\\ &=& \frac{1}{\sqrt{\pi \eta}} \int_{-\infty}^{\infty} e^{-\eta (x+i\lambda)^{2} - \frac{\omega^{2}}{4\eta}} dx \ \ \ \ \ (\lambda \equiv \frac{\omega}{2\eta})\\ &=& \frac{1}{\sqrt{\pi \eta}} \int_{-is \infty}^{\infty} e^{-\eta (x+i\lambda)^{2} - \frac{\omega^{2}}{4\eta}} dx \ \ \ \ \ (\lambda \equiv \frac{\omega}{2\eta})\\ &=& \frac{1}{\sqrt{\pi \eta}} e^{- \frac{\omega^{2}}{4\eta}} \int_{-\infty}^{\infty} e^{-\eta (x+i\lambda)^{2}} dx \end{eqnarray}

積分の部分は g(z) \equiv e^{-\eta z^{2}}の複素積分から求める。 \lambda \gt 0すなわち \omega \gt 0のケースを考えると（ \omega \lt 0の場合も同様である）、 g(z)は図の閉領域内に極値を持たないので、コーシーの積分定理から

\begin{eqnarray} 0 &=& \int_{L1} f(z) dz + \int_{L2} f(z) dz + \int_{L3} f(z) dz + \int_{L4} f(z) dz \\ \end{eqnarray}

となるが、

\begin{eqnarray} \int_{L2} f(z) dz &=& \lim_{N \rightarrow \infty} \int_{y=\lambda}^{0} f(N+iy) dy \\ &=& \lim_{N \rightarrow \infty} \int_{\lambda}^{0} e^{-i\eta (N+iy)^{2}} dy \\ &=& \lim_{N \rightarrow \infty} \int_{\lambda}^{0} e^{-i\eta (N^{2}-y^{2}+iyN)} dy \\ &\rightarrow& 0 \end{eqnarray}

L3についても同様に0になる。したがって、

\begin{eqnarray} F(\omega) &=& \frac{1}{\sqrt{\pi \eta}} e^{- \frac{\omega^{2}}{4\eta}} \int_{L1} f(z) dz \\ &=& \frac{1}{\sqrt{\pi \eta}} e^{- \frac{\omega^{2}}{4\eta}} \int_{-L3} f(z) dz \\ &=& \frac{1}{\sqrt{\pi \eta}} e^{- \frac{\omega^{2}}{4\eta}} \int_{\infty}^{\infty} e^{-\eta x^{2}} dx \\ &=& \frac{1}{\sqrt{\pi \eta}} e^{- \frac{\omega^{2}}{4\eta}} \sqrt{\frac{\pi}{\eta}} \\ &\propto& e^{- \frac{\omega^{2} \sigma^{2}}{2}} \end{eqnarray}

となり、期待値0で標準偏差 \frac{1}{\sigma}のガウシアンに従うことが示せた。
一般に期待値が0ではなく \muの場合を考えると

\begin{eqnarray} F(\omega) &=& \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi \sigma^{2}}} e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}} e^{-i\omega x} dx \\ &=& \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi \sigma^{2}}} e^{-\frac{t^{2}}{2\sigma^{2}}} e^{-i\omega (t+\mu)} dt \\ &\propto & e^{-i\omega \mu} e^{-frac{\omega^{2}\sigma^{2}}{2}} \\ &=& e^{-\frac{\sigma^{2}}{2}(\omega + i\frac{\mu}{\sigma^{2}})} \\ \end{eqnarray}

となり、期待値が -i\frac{\mu}{\sigma^{2}}で標準偏差は変わらず \frac{1}{\sigma^{2}}である。

ガウシアンに従う変数の二乗の期待値

X \sim N(\mu, \sigma^{2})であるとする。このとき、 \lt X^{2} \gt = \sigma^{2}である。

証明 カイ二乗分布の定義より

\begin{eqnarray} \frac{X^{2}}{\sigma^{2}} &\sim& \chi(1) \\ <\frac{X^{2}}{\sigma^{2}}> &=& 1 \\ \lt X^{2} \gt &=& \sigma^{2} \end{eqnarray}

細かい定理

\begin{eqnarray} \sum_{n=0}^{N-1}e^{-i\frac{2\pi nk}{N}} = N\delta_{k, 0} \\ \end{eqnarray}

証明 k=0のとき

\begin{eqnarray} \sum_{n=0}^{N-1} 1 = N \end{eqnarray}

となる。 k \ne 0のとき

\begin{eqnarray} \sum_{n=0}^{N-1}e^{-i\frac{2\pi nk}{N}} &=& \frac{1-e^{-2\pi i k}}{1-e^{-i\frac{2\pi k}{N}}} \\ &=& \frac{1-(e^{-2\pi i})^{k}}{1-e^{-i\frac{2\pi k}{N}}} \\ &=& \frac{1-1}{1-e^{-i\frac{2\pi k}{N}}} \\ &=& 0 \end{eqnarray}

ガウシアン分布に従うノイズの周波数分布

x(0), x(1), \dots, x(N-1)を期待値0で標準偏差 \sigmaのガウシアンに従うノイズとする。このとき、 x(n)をフーリエ変換した振幅の周波数分布は N\sigma^{2}で周波数に対して一定となる。

証明 x(n)のそれぞれがガウシアンに従うので、

\begin{eqnarray} x(n) &\in& N(\mu, \sigma^{2}) \\ P_{X}(x) &=& G(x; 0, \sigma) \end{eqnarray}

このとき、 x(n)のフーリエ変換は

\begin{eqnarray} X(k) &=& \sum_{n=0}^{N-1} x(n) e^{-i\frac{2\pi nk}{N}} \\ \mathrm{Re}\left[ X(k) \right] &=& \sum_{n=0}^{N-1}x(n) \cos\left( \frac{2\pi nk}{N} \right) \\ &\equiv& \sum_{n=0}^{N-1}y(n) \end{eqnarray}

「スケールした確率変数に対する確率密度関数」より、

\begin{eqnarray} P_{Y_{n}} (y) &=& \frac{1}{\cos \left( \frac{2\pi nk}{N} \right)}P_{X}(x) \\ &=& \frac{1}{\cos \left( \frac{2\pi nk}{N} \right)} P_{X}\left(\frac{y}{\cos \left( \frac{2\pi nk}{N} \right)}\right) \\ &=& \frac{1}{\cos \left( \frac{2\pi nk}{N} \right)} G\left(\frac{y}{\cos \left( \frac{2\pi nk}{N} \right)}; 0, \sigma \right) \\ &=& G\left(y; 0, \sigma \cos \left( \frac{2\pi nk}{N} \right) \right) \ \ \ \ \ (スケールしたガウシアンに関する定理) \\ \end{eqnarray}

したがって、

\begin{eqnarray} F\left[ P_{Re[X(k)]} \right](x) &=& \prod_{n=0}^{N-1} F[P_{Y_{n}}](x) \\ &=& \prod_{n=0}^{N-1} F\left [ G\left(y; 0, \sigma \cos \left( \frac{2\pi nk}{N} \right) \right) \right ] (x) \\ &=& \prod_{n=0}^{N-1} G\left(y; 0, \frac{1}{\sigma \cos \left( \frac{2\pi nk}{N} \right)} \right) (x) \\ &=& \prod_{n=0}^{N-1} \sqrt{\frac{\sigma^{2} \cos^{2} \left( \frac{2\pi nk}{N} \right)}{2\pi}} \exp \left( -\frac{x^{2} \sigma^{2} \cos^{2} \left( \frac{2\pi nk}{N} \right) }{2} \right) \\ &=& \left( \frac{\sigma^{2}}{2\pi} \right)^{\frac{N}{2}} \prod_{n=0}^{N-1} \cos \left( \frac{2\pi nk}{N} \right) \exp \left( -\frac{x^{2} \sigma^{2} \cos^{2} \left( \frac{2\pi nk}{N} \right) }{2} \right) \\ &=& \left( \frac{\sigma^{2}}{2\pi} \right)^{\frac{N}{2}} \prod_{n=0}^{N-1} \cos \left( \frac{2\pi nk}{N} \right) \exp \left( -\frac{x^{2} \sigma^{2}}{2} \sum_{n=0}^{N-1} \cos^{2} \left( \frac{2\pi nk}{N} \right) \right) \\ &=& \left( \frac{\sigma^{2}}{2\pi} \right)^{\frac{N}{2}} \prod_{n=0}^{N-1} \cos \left( \frac{2\pi nk}{N} \right) \exp \left( -\frac{x^{2} \sigma^{2}}{2} \sum_{n=0}^{N-1} \left( \frac{1}{2} + \frac{1}{2} \cos \left( \frac{4\pi nk}{N} \right) \right) \right) \\ \\ &=& \left( \frac{\sigma^{2}}{2\pi} \right)^{\frac{N}{2}} \prod_{n=0}^{N-1} \cos \left( \frac{2\pi nk}{N} \right) \exp \left( -\frac{x^{2} \sigma^{2}}{2} \left[ \frac{N}{2} + \frac{N}{2}\delta_{k, 0} \right] \right) \\ \end{eqnarray}

\expの部分のみ着目すると、 k=0のとき

\begin{eqnarray} F\left[ P_{Re[X(0)]} \right](x) &\sim& \exp \left( - \frac{N\sigma^{2}x}{2}\right) \\ &\sim& G(x; 0, \frac{1}{\sqrt{N}\sigma} ) \\ P_{Re[X(0)]} &\sim & G(x; 0, \sqrt{N}\sigma ) \ \ \ (ガウシアンのフーリエ変換もガウシアンであることを利用) \\ \end{eqnarray}

となり、 k>0のときは

\begin{eqnarray} F\left[ P_{Re[X(k)]} \right](x) &\sim& \exp \left( - \frac{N\sigma^{2}x}{4}\right) \\ &\sim& G(x; 0, \sqrt{\frac{2}{N}}\frac{1}{\sigma} ) \\ P_{Re[X(k)]} &\sim & G(x; 0, \sqrt{\frac{N}{2}}\sigma ) \ \ \ (ガウシアンのフーリエ変換もガウシアンであることを利用) \\ \end{eqnarray}

となる。またこの議論は虚部に対してもほぼ同様で

\begin{eqnarray} P_{Im[X(0)]} (x) &=& 0 \\ P_{Im[X(k)]} (x) &\sim& G(x; 0, \sqrt{\frac{N}{2}}\sigma ) \ \ \ (k>0) \\ \end{eqnarray}

以上の議論から

\begin{eqnarray} \lt |X(0)|^{2} \gt &=& \lt Re[ X(0)]^{2} + Im [ X(0) ]^{2} \gt \\ &=& \lt Re[ X(0)]^{2} \gt + \lt Im[ X(0)]^{2} \gt \\ &=& N\sigma^{2} \\ \lt |X(k)|^{2} \gt & = & \lt Re[ X(k)]^{2} + Im [ X(k) ]^{2} \gt \\ &=& \lt Re[ X(k)]^{2} \gt + \lt Im[ X(k)]^{2} \gt \\ &=& \frac{N}{2}\sigma^{2} + \frac{N}{2}\sigma^{2} \\ &=& N\sigma^{2} \end{eqnarray}

となり、振幅2乗の期待値は一定値 N\sigma^{2}になることが示せた。

もし従うガウシアンの期待値が \mu (\ne 0)であった場合、すなわち x^{'}(n) \sim G(x; \mu, \sigma)だった場合は、そのフーリエ変換 X^{'}(k)は k=0のみピークを持つ形になる。

\begin{eqnarray} X^{'}(k) &=& \sum_{n=0}^{N-1}x^{'}(n)e^{-i\frac{2\pi nk}{N}} \\ &=& \sum_{n=0}^{N-1}(x(n)+\mu) e^{-i\frac{2\pi nk}{N}} \\ &=& X(k) + \mu \sum_{n=0}^{N-1} e^{-i\frac{2\pi nk}{N}} \\ &=& X(k) + N\mu \delta_{k, 0} \end{eqnarray}

実験

実験に使用したコード

import matplotlib.pyplot as plt
import numpy as np


N = 900
mu = 0.3
sigma = 1
dt=0.1

array = np.random.normal(mu, sigma, N)

# Fourier transform                                                                        
f = np.fft.fft(array)
amp = np.abs(f)

fig = plt.figure()
ax1 = fig.add_subplot(2, 1, 1)
ax1.plot([i for i in range(N)], array)
ax1.set_title(f'Gaussian (mean={mu}, sigma={sigma}, N={N}, dt={dt})')
ax1.set_xlabel('Step')
ax1.set_ylabel('Value')
ax2 = fig.add_subplot(2, 1, 2)
freq = np.fft.fftfreq(len(amp), d=dt)
ax2.plot(freq, amp)
ax2.set_xlabel('Frequency')
ax2.set_ylabel('Amplitude')
plt.show()

f:id:salpik:20210807115138p:plain f:id:salpik:20210807115145p:plain f:id:salpik:20210807115153p:plain f:id:salpik:20210807115200p:plain

関係図

f:id:salpik:20210719231305p:plain
関係図

線形時不変システム

線形時不変システムに正弦波を入力すると、同じ周波数の正弦波が出力される。

系の安定性

伝達関数の極の実部の符号が系の安定性を決める。負であれば安定、正であれば発散する。

周波数特性

。伝達関数が G(s)であるとし、その極が p_{1}, p_{2}, \dots, p_{n}であるとする。安定な形であるとして Re(p_{1}) \lt 0, Re(p_{2}) \lt 0, \dots, Re(p_{n}) \lt 0が成立するとする。この場合 G(s)は

\begin{eqnarray} G(s) = \sum_{i=1}^{n}\frac{K_{i}}{s-p_{i}} \end{eqnarray}

と変形できる。この系に正弦波を入力したときの応答が周波数特性を考えてみると

\begin{eqnarray} y(t) &=& \mathcal{L}^{-1}[G(s)\mathcal{L}[u(t)]] \\ &=& \mathcal{L}^{-1}[G(s)\mathcal{L}[e^{i\omega t}]] \\ &=& \mathcal{L}^{-1}[G(s) \frac{1}{s-i\omega}] \ \ \ \ \ (ラプラス変換を実施)\\ &=& \mathcal{L}^{-1}\left[\frac{K_{0}}{s-i\omega} + \sum_{i=1}^{n}\frac{K_{i}}{s-p_{i}} \right] \\ &=& K_{0}e^{i\omega t} + \sum_{i=1}^{n}K_{i}e^{p_{i}t} \end{eqnarray}

となるが、 p_{i}の実部は負であるから t \rightarrow \inftyにおいて e^{p_{i}t} \rightarrow 0となるので無視すると

\begin{eqnarray} y(t) &=& K_{0}e^{i\omega t} \\ \end{eqnarray}

である。 K_{0}は

\begin{eqnarray} K_{0} &=& \lim_{s\rightarrow i\omega} [ (s-i\omega)G(s) \frac{1}{s-i\omega} ] \\ &=& G(i\omega) \end{eqnarray}

として求まる。以上をまとめると

\begin{eqnarray} y(t) = G(i\omega)e^{i\omega t} \end{eqnarray}

となるので G(i\omega)が周波数特性となる。

参考

https://www.ishikawa-nct.ac.jp/lab/E/y_kawai/www/data/course/CE2/21CE2/handouts/21CE2_s_lect01/21CE2_s_lect01_slide.pdf

一次遅れの系

入出力が以下のような微分方程式で表される系を一次遅れと呼ぶ。

\begin{eqnarray} \tau \frac{dy(t)}{dt} + y(t) = Kx(t) \end{eqnarray}

時間ステップが \Delta tであるような場合、この微分方程式を置き換えると

\begin{eqnarray} \tau \frac{y(t+\Delta t) - y(t)}{\Delta t} + y(t) &=& Kx(t) \\ y(t+\Delta t) &=& \left( 1 - \frac{\Delta t}{\tau} \right) y(t) + K\frac{\Delta t}{\tau} x(t) \end{eqnarray}

という漸化式となる。
この一次遅れの系の伝達関数 G(s)を求めてみると

\begin{eqnarray} \tau (sY(s)-Y(0)) + Y(s) &=& KX(s) \\ (\tau s + 1) Y(s) &=& KX(s) \ \ \ \ \ (Y(0)=0) \\ G(s) &=& \frac{Y(s)}{X(s)} \\ &=& \frac{K}{1+\tau s} \end{eqnarray}

したがって、この系の周波数特性は

\begin{eqnarray} G(i\omega) &=& \frac{K}{1+i\tau \omega} \\ &=& \frac{K}{\sqrt{1+(\omega \tau)^{2}}e^{i\phi}} \\ &=& \frac{K}{\sqrt{1+(\omega \tau)^{2}}} e^{-i\phi} \end{eqnarray}

となる。

参考

dsp.stackexchange.com