2021-08-07

原始根

定義

ある素数 $p$ に対して、 $r, r^{2}, r^{3}, \dots, r^{p-1}\ (\mathrm{mod}\ p)$ が互いに異なるような $r$ のことを $p$ に対する原始根という。例えば、 $p=7$ に対しては

$\begin{eqnarray} 3 &\equiv& 3\ (\mathrm{mod}\ 7) \\ 3^{2} &\equiv& 2\ (\mathrm{mod}\ 7) \\ 3^{3} &\equiv& 6\ (\mathrm{mod}\ 7) \\ 3^{4} &\equiv& 4\ (\mathrm{mod}\ 7) \\ 3^{5} &\equiv& 5\ (\mathrm{mod}\ 7) \\ 3^{6} &\equiv& 1\ (\mathrm{mod}\ 7) \\ \end{eqnarray}$

となるので3は原始根であるが、

$\begin{eqnarray} 2 &\equiv& 2\ (\mathrm{mod}\ 7) \\ 2^{2} &\equiv& 4\ (\mathrm{mod}\ 7) \\ 2^{3} &\equiv& 1\ (\mathrm{mod}\ 7) \\ \end{eqnarray}$

となるので2は原子根ではない。

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

関連する定理

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

連続値をとる確率変数 $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

関係図

線形時不変システム

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

系の安定性

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

周波数特性

。伝達関数が $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

2021-06-20

ABC206 E Divide Both

問題

atcoder.jp

方針

互いに素でないペアの数を求める.
(5, 5), (5, 10), (10, 15) など
互いに素でないペアの中から、どちらかの値＝最大公約数となるペアの数を求める.
(5, 5), (5, 10) など（(10, 15)は列挙されない）
1の数から2の数を引いたものが答えになる

互いに素でないペアの数を求める

2つのやり方がある

メビウス関数を使う方法

メビウス関数 $\mu(n)$ は以下のように定義される。
- ある素数 $p$ が存在して $n$ が $p^2$ で割り切れる時 $\mu(n)=0$ .
- それ以外の時、 $n$ の素因数の個数を $k$ 個とすると $\mu(n)=(-1)^k$

メビウス関数には以下の性質がある.
任意の自然数 $n$ に対して

$\begin{eqnarray} \sum_{d|n}\mu(d) = 0 \end{eqnarray}$

ここで $d|n$ は $n$ に対する全ての約数を表す。これを変形すると

$\begin{eqnarray} -\sum_{d|n かつ d>2}\mu(d) = \mu(1) = 1 \end{eqnarray}$

となることがわかる。この性質を利用して、 $R$ 以下の任意の2つの数 $i$ 、 $j$ に対する全てのペア $(i, j)$ の中から互いに素でないものの個数を数える。

$\begin{eqnarray} \sum_{g=2}^{R} (gを公約数に持つペアの数) \times (-\mu(g)) \end{eqnarray}$

というものを考えると、この中において同じ項が複数のgで扱われることがある。例えば(6, 12)というペアは(1を除いた)公約数は2, 3, 6なので、g=2, 3, 6の時に出現する。この係数の和を考えるとこれは $-\sum_{d|n かつ d>2}\mu(d)$ となるので1となる。したがって、上で示した式において、「互いに素ではないペア」の係数が全て1となることが示されたので、この式が互いに素でないペアの数を表していることが示された。

漸化式的に求める方法

$F(g)$ を公約数（最大公約数ではない）がgであるペアの個数とし、 $G(g)$ を最大公約数がgであるペアの個数とする。 $F(g)$ は簡単に求まり、 $G(g)$ の方は

$\begin{eqnarray} G(g) = F(g) - G(2g) - G(3g) - \dots \end{eqnarray}$

となるので、自分より大きい $g$ に対して $G(g)$ が既に求まっているのであればこの漸化式で求められる。したがってgの大きい方から $G(g)$ を求めることができ、ここから互いに素でないペアの数が求まる。

参考

manabitimes.jp

2021-05-29

Dynamic Window Approach

元論文

www.researchgate.net

運動学の確認

本論文で扱われている、進行速度と角速度がコントロールできる平面上のロボットの運動学を整理する。時刻 $t$ におけるロボットの位置を $x(t), y(t)$ とすると

$\begin{eqnarray} x(t_{n}) &=& x(t_{0}) + \int_{t_{0}}^{t_{n}} v(t^{'}) \cos \theta(t^{'}) dt^{'} \\ &=& x(t_{0}) + \sum_{i=0}^{n-1} \int_{t_{i}}^{t_{i+1}} v(t^{'}) \cos \theta(t^{'}) dt^{'} \\ &=& x(t_{0}) + \sum_{i=0}^{n-1} \int_{t_{i}}^{t_{i+1}} \left[v(t_{i}) + \dot{v}({t_{i}}) \right] \cos \left[ \theta(t_{i}) + \omega(t_{i})(t^{'} - t_{i}) + \frac{1}{2} \dot{\omega}(t_{i})(t^{'} - t_{i})^{2} \right] dt^{'} \end{eqnarray}$

ここで、進行速度と角速度の変化率は $t_{i} \le t \le t_{i+1}$ の間で0であるとすると、 $\dot{v}(t_{i}) = 0, \dot{w}(t_{i}) = 0$ となり、

$\begin{eqnarray} x(t_{n}) &\sim& x(t_{0}) + \sum_{i=0}^{n-1} \int_{t_{i}}^{t_{i+1}} v(t_{i}) \cos \left[ \theta(t_{i}) + \omega(t_{i})(t^{'} - t_{i}) \right] dt^{'} \\ &\equiv& x(t_{0}) + \sum_{i=0}^{n-1} F_{x}^{i}(t_{i+1}) \end{eqnarray}$

$F_{x}^{i}(t_{i+1})$ は積分可能であり、 $w(t_{i}) \ne 0$ のとき、

$\begin{eqnarray} F_{x}^{i}(t_{i+1}) &=& \int_{t_{i}}^{t_{i+1}} v(t_{i}) \cos \left[ \theta(t_{i}) + \omega(t_{i})(t^{'} - t_{i}) \right] dt^{'} \\ &=& \frac{v(t_{i})}{\omega(t_{i})} \left[ \sin(\theta(t_{i}) + \omega(t_{i})(t_{i+1} - t_{i})) - \sin \theta(t_{i})\right] \end{eqnarray}$

$w(t_{i}) = 0$ のとき、

$\begin{eqnarray} F_{x}^{i}(t_{i+1}) &=& \int_{t_{i}}^{t_{i+1}} v(t_{i}) \cos \theta(t_{i}) dt^{'} \\ &=& v(t_{i}) (t_{i+1} - t_{i}) \cos \theta(t_{i}) \end{eqnarray}$

同様の計算をすれば、

$\begin{eqnarray} x(t_{n}) &\sim& y(t_{0}) + \sum_{i=0}^{n-1} F_{y}^{i}(t_{i+1}) \\ F_{y}^{i}(t_{i+1}) &=& \int_{t_{i}}^{t_{i+1}} v(t_{i}) \sin \left[ \theta(t_{i}) + \omega(t_{i})(t^{'} - t_{i}) \right] dt^{'} \\ &=& -\frac{v(t_{i})}{\omega(t_{i})} \left[ \cos(\theta(t_{i}) + \omega(t_{i})(t_{i+1} - t_{i})) - \cos \theta(t_{i})\right]\ \ \ \ \ (\omega(t_{i}) \ne 0) \\ &=& v(t_{i}) (t_{i+1} - t_{i}) \sin \theta(t_{i})\ \ \ \ \ (\omega(t_{i}) = 0) \end{eqnarray}$

2021-05-29

一目で思い出すUML

UML記法が一目で思い出せるように整理した。

クラス図

シーケンス図

2021-05-16

GPS解析

楕円におけるパラメータの定義

f:id:salpik:20210516230705p:plain — 楕円のパラメータ定義

f:id:salpik:20210516231732p:plain — 起動面の定義に関わるパラメータ（図はWikipediaより）

楕円軌道の軌道面を表すパラメータ

昇降点赤経 $\Omega$
起動傾斜角度 $i$
近地点引数 $\omega$ 楕円の長軸が軌道面のうちどの方向にあるかを規定する

楕円の形状を表すパラメータ

長軸
離心率

楕円軌道における位置を表すパラメータ

真近点角

解くべき方程式

$(x_{i}, y_{i}, z_{i})$ にある衛星 $i (i=1, 2, \dots N)$ から対象物までの距離の測定結果 $r_{i}$ が得られたとする。対象物の座標 $(x, y, z)$ を求めたい。方程式は

$\begin{eqnarray} r_{1}^{2} &=& (x-x_{1})^{2} + (y-y_{1})^{2} + (z-z_{1})^{2} \\ r_{2}^{2} &=& (x-x_{2})^{2} + (y-y_{2})^{2} + (z-z_{2})^{2} \\ r_{3}^{2} &=& (x-x_{3})^{2} + (y-y_{3})^{2} + (z-z_{3})^{2} \\ \cdots \end{eqnarray}$

である。非線形方程式なので解析的に解くのは難しいので、初期値から収束させる方法で解を探索する。初期解を $(x_{0}, y_{0}, z_{0})$ し、この位置において得られるはずの測定距離を $r_{i}^{'}$ とすると

$\begin{eqnarray} r_{1}^{'2} &=& (x_{0}-x_{1})^{2} + (y_{0}-y_{1})^{2} + (z_{0}-z_{1})^{2} \\ r_{2}^{'2} &=& (x_{0}-x_{2})^{2} + (y_{0}-y_{2})^{2} + (z_{0}-z_{2})^{2} \\ r_{3}^{'2} &=& (x_{0}-x_{3})^{2} + (y_{0}-y_{3})^{2} + (z_{0}-z_{3})^{2} \\ \cdots \end{eqnarray}$

$r_{i}^{'}$ は実際の測定結果 $r_{i}$ とは差分があるはずなので

$\begin{eqnarray} \Delta r_{1} \equiv r_{1}^{'2} - r_{1}^{'2} \\ \Delta r_{2} \equiv r_{2}^{'2} - r_{2}^{'2} \\ \cdots \end{eqnarray}$

が求まる。一方

$\begin{eqnarray} \Delta r_{1} = \frac{\partial r_{1}}{\partial x}\Delta x + \frac{\partial r_{1}}{\partial y}\Delta y + \frac{\partial r_{1}}{\partial z}\Delta z \\ \Delta r_{2} = \frac{\partial r_{2}}{\partial x}\Delta x + \frac{\partial r_{2}}{\partial y}\Delta y + \frac{\partial r_{2}}{\partial z}\Delta z \\ \cdots \end{eqnarray}$

この線形方程式を解いて $\Delta x, \Delta y, \Delta z$ が求まるので、これで初期解を更新できる。このステップを繰り返すことで厳密解に近い解を得ることができる。

エポック時刻からのずれの補正

受信したメッセージの時刻は、エポック時刻 $t^{epoc}$ を基準とした二次の補正量が必要となる。すなわち、真の時刻 $t^{true}$ と観測したメッセージの時刻[t^{measure}]の間には

$\begin{eqnarray} t^{measure} = t^{true} + a_{0} + a_{1}(t^{true}-t^{epoc}) + a_{2}(t^{true}-t^{epoc})^{2} \end{eqnarray}$

という関係式がある。補正係数は受信メッセージに含まれる。

メモ

初期化時に衛星探索をしやすいように、各衛星のGPSメッセージに全ての衛星の位置情報が含まれる。
電波が電離層を通過する時に電子密度に比例して遅延（電離層遅延）が発生する。遅延量はアップルトンハートリーの公式で計算され、メッセージには補正係数が含まれる。

参考

GPSのための実用プログラミング | 坂井丈泰 |本 | 通販 | Amazon

2021-03-11

中国剰余定理

定理の内容

$a_{1}, a_{2}, \dots , a_{n}$ および $m_{1}, m_{2}, \dots , m_{n}$ があるとき、

$\begin{eqnarray} x &\equiv& a_{1} (\mathrm{mod}\ m_{1}) \\ x &\equiv& a_{2} (\mathrm{mod}\ m_{2}) \\ &\vdots& \\ x &\equiv& a_{n} (\mathrm{mod}\ m_{n}) \\ \end{eqnarray}$

を全て満たす整数 $x$ が存在する必要十分条件は

$\begin{eqnarray} a_{1} \equiv a_{2} \equiv \dots \equiv a_{n} (\mathrm{mod}\ \mathrm{LCM}(m_{1}, m_{2}, \cdots , m_{n})) \end{eqnarray}$

であり、このとき $0 \le x \lt \mathrm{LCM}(m_{1}, m_{2}, \dots , m_{n})$ の範囲にただ一つの解 $x$ が存在する。

互いに素の2つのペアのケース

定理の内容を言い換えると、 $m_{1}, m_{2}$ が互いに素のとき、

$\begin{eqnarray} x &\equiv& a_{1} (\mathrm{mod}\ m_{1}) \\ x &\equiv& a_{2} (\mathrm{mod}\ m_{2}) \\ \end{eqnarray}$

を満たす解 $x$ が $0 \le x \lt m_{1} \times m_{2}$ の範囲にただ一つ存在する。

証明

解の存在性

拡張ユークリッドの互除法を用いることで、

$\begin{eqnarray} km_{1} + hm_{2} = 1 \end{eqnarray}$

となる整数 $k, h$ が存在する。これを変形すると

$\begin{eqnarray} km_{1} = 1 - hm_{2} \equiv 1\ (\mathrm{mod}\ m_{2}) \end{eqnarray}$

となるので

$\begin{eqnarray} a_{2}km_{1} \equiv a_{2}\ (\mathrm{mod}\ m_{2}) \end{eqnarray}$

となる。同様に

$\begin{eqnarray} a_{1}hm_{2} \equiv a_{1}\ (\mathrm{mod}\ m_{1}) \end{eqnarray}$

なので、

$\begin{eqnarray} x = a_{2}km_{1} + a_{1}hm_{2} \end{eqnarray}$

とすれば

$\begin{eqnarray} x &\equiv& a_{1} (\mathrm{mod}\ m_{1}) \\ x &\equiv& a_{2} (\mathrm{mod}\ m_{2}) \\ \end{eqnarray}$

を満たすことがわかる。

解の一意性

$0\le x_{1} \lt x_{2} \lt m_{1}m_{2}$ という2つの解 $x_{1}, x_{2}$ が存在したとする。このとき、

$\begin{eqnarray} x_{1} &\equiv& x_{2}\ (\mathrm{mod}\ m_{1}) \\ x_{1} &\equiv& x_{2}\ (\mathrm{mod}\ m_{2}) \\ \end{eqnarray}$

であるから

$\begin{eqnarray} x_{2} - x_{1} &=& rm_{1} \\ x_{2} - x_{1} &=& sm_{2} \\ \therefore rm_{1} &=& sm_{2} \end{eqnarray}$

となり、 $m_{1}$ と $m_{2}$ は互いに素であるから $r$ は $m_{2}$ の倍数である必要がある。一方解の範囲の条件より $x_{2} - x_{1} = rm_{1} \lt m_{1}m_{2}$ なので、必然的に $r=0$ となり、 $x_{1} = x_{2}$ である。

互いに素ではない2つのペアのケース

$\mathrm{GCD}(m_{1},\ m_{2})=g \gt 1$ のとき、

$\begin{eqnarray} x &\equiv& a_{1} (\mathrm{mod}\ m_{1}) \\ x &\equiv& a_{2} (\mathrm{mod}\ m_{2}) \\ \end{eqnarray}$

の解が存在する必要十分条件は、

$\begin{eqnarray} a_{1} &\equiv& a_{2} (\mathrm{mod} g) \end{eqnarray}$

であり、解は $0 \le x \lt LCM(m_{1},\ m_{2})$ の範囲に存在する。

解の存在性

拡張ユークリッドの互除法により、

$\begin{eqnarray} km_{1} + hm_{2} = g \end{eqnarray}$

となる整数 $k, h$ が存在する。また、 $a_{1} \equiv a_{2}\ (\mathrm{mod}\ g)$ なので、 $a_{2}-a_{1}=ng$ と表されるから、

$\begin{eqnarray} nkm_{1} + nhm_{2} = a_{2} - a_{1} \end{eqnarray}$

ここで $x=nkm_{1}+a_{1}$ と定義すると

$\begin{eqnarray} x=nkm_{1}+a_{1} = -nhm_{2} + a_{2} \end{eqnarray}$

となることから、

$\begin{eqnarray} x &\equiv& a_{1}\ (\mathrm{mod}\ m_{1}) \\ x &\equiv& a_{2}\ (\mathrm{mod}\ m_{2}) \\ \end{eqnarray}$

となることが示せる。また、もし $x$ が $\mathrm{LCM}(m_{1},\ m_{2})$ 以上であれば、 $x%LCM(m_{1},\ m_{2})$ を行っても $\mathrm{mod}\ m_{1}$ と $\mathrm{mod}\ m_{2}$ の値は変わらないので、解が $0 \le x \lt LCM(m_{1},\ m_{2})$ の範囲に存在することもわかる。