サンプリング定理の証明

サンプリング定理の証明

 最近勉強しているとサンプリング定理(標本化定理)が出てきて,それについての証明をしてみたので未来の自分に向けてメモを残しておこうと思う.
 今回はこういうふうに考える.ある信号波S(t)に対して搬送波として矩形波p(t)を用いて変調を行い,PAM波f(t)=S(t)p(t)を得る(パルス振幅変調:PAM).このPAM波を復調したときに信号波を完全に復元できるためには搬送波である矩形波の周波数は,
\displaystyle{\omega_c\ge2\omega_{max}}
\omega_c矩形波の周波数.\omega_{max}:信号波に含まれる周波数の最大値)
を満たす必要がある.これを示してみる.

f:id:Ibetoge:20210503233513p:plain
信号波,矩形波,PAM波

証明

 まず,搬送波として振幅A,周期T=2\pi/\omega_c,一周期内での非ゼロの時間\tauである矩形波p(t)を考える.

f:id:Ibetoge:20210503215342p:plain
矩形波

この矩形波フーリエ級数展開すると,以下のようになる(矩形波フーリエ級数展開のやり方はこのページの後ろを参照).
\displaystyle{p(x)=\frac{A\omega_c\tau}{\pi}\left(1+\sum_{m=1}^{\infty}Sinc(m\pi\tau/2)\cos(m\omega_ct)\right)}
 さらに,信号波をフーリエ展開すると以下のようになるとする.
\displaystyle{S(t)=\sum_1^N\left\{\alpha_n\cos(n\omega_st)+\beta_n\sin(n\omega_st)\right\}}
このとき信号波に含まれる周波数の最大値は,\omega_{max}=N\omega_sとなる.
そうするとPAM波が得られ,
\displaystyle{
f(t)=S(t)p(t)=S(t)\frac{A\omega_c\tau}{\pi}\left(1+\sum_{m=1}^{\infty}Sinc(m\pi\tau/2)\cos(m\omega_ct)\right)}

\displaystyle{
=\frac{A\omega_c\tau}{\pi}S(t)+\frac{A\omega_c\tau}{\pi}S(t)\sum_{m=1}^{\infty}Sinc(m\pi\tau/2)\cos(m\omega_ct)}
となる.ここで,式を見やすくするために B=\frac{A\omega_c\tau}{\pi}, a_m=Sinc(m\pi\tau/2)と置いて,式の第二項を展開すると,

\displaystyle{
f(t)=BS(t)+B\sum_{m=1}^{\infty}\sum_{n=1}^N\left\{a_m\alpha_n\cos(m\omega_ct)\cos(n\omega_st)+a_m\beta_n\cos(m\omega_ct)\sin(n\omega_st)\right\}}
\displaystyle{
=BS(t)+}
\displaystyle{
B\sum_{m=1}^{\infty}\sum_{n=1}^N\left\{\frac{a_m\alpha_n}{2}(\cos(m\omega_ct+n\omega_st)+\cos(m\omega_ct-n\omega_st))+\frac{a_m\beta_n}{2}(\sin(m\omega_ct+n\omega_st)-\sin(m\omega_ct-n\omega_st))\right\}}
さて,ここで上式の第一項に関しては,Bの値が既知なので第二項さえ取り除くことができれば,PAM波をもとの信号波に完全に復元できたと言える.ではどのような手法で第二項を取り除くか?またその手法を用いることのできる条件は何か?
 ここで用いるのが低域通過フィルタである.これは,合成波のなかで,ある周波数以上の周波数成分を取り除く電子回路である(このページの後ろに簡単な説明を載せときます).ではここではいくつ以上の周波数成分を取り除けばよいかというと,もちろん\omega\gt\omega_{max}となる成分を取り除けば良い.なぜなら,\omega_{max}以下の周波数成分を取り除くと上式の第一項の成分も取り除かれて,PAM波から信号波を完全に取り除くことができなくなるからである.更にいうと,\omega\gt\omega_{max}を取り除く低域通過フィルタを作用させたときに第二項が完全に除去される必要がある.つまり,第二項の成分のうち,周波成の最小値が\omega\gt\omega_{max}を満たしていれば良い.では,上式第二項の最小周波数はいくつなのか.整数mの範囲は1\le m \le \inftyで,n1\le n \le Nなので,最小周波数はm\omega_c-n\omega_sm=1, n=Nを代入した\omega_c-N\omega_s=\omega_c-\omega_{max}となる.したがって,低域通過フィルタを用いてPAM波から信号波を完全に復調できる条件は,低域通過フィルタ\omega\gt\omega_{max}で上式第二項の最小周波数成分\omega_c-\omega_{max}を除去できることであり,数式にすると,
\displaystyle{\omega_c-\omega_{max}\gt\omega_{max}\Longleftrightarrow\omega_c\gt2\omega_{max}}
つまり,搬送波として用いる矩形周波数は,信号波に含まれる周波数の最大値の二倍より大きくなくてはならない.また,\omega_c\gt2\omega_{max}を満たすような周波数\omega_cナイキスト周波数とよんだりする.

エイリアシング

 写真などでよくわからない周期的な模様ができたりする.これをエイリアシングとか折返し雑音とか呼んだりする.これは上の式で矩形波の周波数がナイキスト周波数より小さかったとき,第二項が十分に取り除かれずに残ってしまったときに生じるものであるらしい.この場合上の式のt が写真内の距離xに対応して,周波数はピクセル数に対応するのだと思う.

矩形波フーリエ展開

 矩形波は周期Tとして,区間-T/2\le t\le T/2において,
\displaystyle{
p(t) = \left\{
\begin{array}{cc}
A,& -\tau/2\le t\le\tau/2  \\
0, &  otherwise
\end{array}
\right.
}
で表される.周期関数をフーリエ展開する場合,
\displaystyle{
f(t)=a_0+\sum_{n=1}^\infty(a_n\cos(2n\pi t/T)+b_n\sin(2n\pi t/T)),
}
\displaystyle{
a_n=\frac{2}{T}\int_{-T/2}^{T/2}f(t)\cos(2n\pi t/T)dt,
}
\displaystyle{
b_n=\frac{2}{T}\int_{-T/2}^{T/2}f(t)\sin(2n\pi t/T)dt,
}
である.ここで,上の式のf(x)矩形波p(x)を代入するが,p(x)は偶関数なのでb_n=0となる.次にa_n を求める.
\displaystyle{a_n=\frac{2}{T}\int_{-T/2}^{T/2}p(t)\cos(2n\pi t/T)dt}
\displaystyle{=\frac{2}{T}\int_{-\tau/2}^{\tau/2}A\cos(2n\pi t/T)dt}
\displaystyle{=\frac{A}{n\pi}\left\{\sin(n\pi\tau/T)-sin(-n\pi\tau/T)\right\}}
\displaystyle{=\frac{2A}{n\pi}\sin(n\pi\tau/T)}
\displaystyle{=\frac{2A}{n\pi}\sin(n\pi\tau/T)\cdot\frac{\tau/T}{\tau/T}}
\displaystyle{=\frac{2A\tau}{T}\frac{\sin(n\pi\tau/T)}{n\pi\tau/T}}
\displaystyle{=\frac{2A\tau}{T}Sinc(n\pi\tau/T)}

低域通過フィルター

低域通過フィルタはローパスフィルタ,LPSなどと呼ばれる電子回路である.WIkipediaによると,「ローパスフィルタ(英語: Low-pass filter: LPF、低域通過濾波器)とは、フィルタの一種で、なんらかの信号のうち、遮断周波数より低い周波数の成分はほとんど減衰させず、遮断周波数より高い周波数の成分を逓減させるフィルタである。ハイカットフィルタ等と呼ぶ場合もある。電気回路・電子回路では、フィルタ回路の一種である。」らしい.

Plotlyで線グラフをたくさん描く

Plotlyとは?

Pythonのグラフ描画ライブラリーで描画したグラフの拡大縮小,軸のスケールの変更などなどができる.
Plotlyのグラフ描画モジュールにはexpressとgraph_objectsがあるらしいが,基本はgraph_objectsでexpressモジュールはサクッとグラフを書くためのモジュールで内部ではgraph_objectsの関数を呼んでいるらしい.(公式ドキュメントPlotly Express | Python | Plotlyの雑和訳)

f:id:Ibetoge:20210427124756g:plain
Plotlyの動くグラフ

今回はgraph_objectsでいい感じのグラフを色々書いてみる

expressでの描写はこちら→
ibetoge.hatenablog.com

線グラフ

基本的なグラフを描写してみる.とりあえず使うモジュールをインポートする.

import numpy as np
import plotly.graph_objects as go

ひとまず0\le t\le2\pi\sin{4t},e^tを描く.

# データを用意
t = np.linspace(0, 2*np.pi, 10)
sint = np.sin(4*t)
ex = np.exp(t/np.pi)

# グラフのインスタンス作成
fig = go.Figure(data=[
    go.Scatter(x=t, y=sint, name=r'$\sin t$'),
    go.Scatter(x=t, y=ex, name=r'$e^t$')
])

# 動かせるグラフをブラウザーで表示
fig.show()

基本的にはこれでグラフが描ける.

f:id:Ibetoge:20210427170338p:plain
基本的な線グラフ

あとは,軸に名前をつけるのと,作成したグラフのインスタンスに新しいグラフを追加してみたい.

# グラフ追加
fig.add_trace(go.Scatter(x=t, y=t**2/4, name=r'$t^2/4$', mode='markers'))  # ドットのみ
fig.add_trace(go.Scatter(x=t, y=sint*ex, name=r'$e^t\sin t$', mode='lines'))  # 線のみ
fig.add_trace(go.Scatter(x=t, y=1/t, name=r'$1/t$', mode='lines+markers'))  # ドットと線

# 軸の名前を追加
fig.update_xaxes(title="t")
fig.update_yaxes(title=r"$f(t)$")

こんな感じで追加できる.このとき,modeで線グラフの種類を選べる.modeを指定しないとPlotlyがいい感じに決めてくれる.グラフを追加できるので,まず空のインスタンスを作ってから追加していってもOK.つまり,

# 空のインスタンス作成
fig = go.Figure()

# グラフ追加
fig.add_trace(go.Scatter(x=t, y=t**2/4, name=r'$t^2/4$')) 
fig.add_trace(go.Scatter(x=t, y=sint*ex, name=r'$e^t\sin t$'))
fig.add_trace(go.Scatter(x=t, y=1/t, name=r'$1/t$' )) 

fig.show()

みたいな感じで書いてもいい.

f:id:Ibetoge:20210427171237p:plain
いろいろな線グラフ
f:id:Ibetoge:20210427171830p:plain
いろいろな線グラフ(横軸のメッシュ細かめ)

Plotlyのグラフをファイルに保存する

fig.show()

を用いると,実行した際に動かせるグラフがブラウザーに表示されるが,グラフをhtmlファイルとして保存したい場合は,

#htmlで保存
fig.write_html("figure.html")

でOK.画像ファイルとして保存するには,kaleidoパッケージが必要なので,次のどちらかのコマンドを叩いてkaleidoをインストールする.(kaleidoの代わりにorcaパッケージでも良いらしいが非推奨)
pipの場合

 pip install -U kaleido

anacondaの場合

 conda install -c conda-forge python-kaleido

kaleidoをインストールしたら,

#pngで保存
fig.write_image("figure.png")
#jpegで保存
fig.write_image("figure.jpeg")

と書けばすれば良し.こうすれば保存したhtmlファイルをブラウザで開けばいつでも動かせるグラフが見られる.もちろん画像ファイルとして保存した場合は動かせない.

棒グラフとは円グラフとかも書こうと思ったけど疲れたからここまで

ソースコード

いろいろな線グラフ

import numpy as np
import plotly.graph_objects as go

# データを用意
t = np.linspace(0, 2*np.pi, 100)
sint = np.sin(4*t)
ex = np.exp(t/np.pi)

# グラフの空インスタンス作成
fig = go.Figure()

# グラフ追加
fig.add_trace(go.Scatter(
    x=t, y=t**2/4, name=r'$t^2/4$', mode='markers'))  # ドットのみ
fig.add_trace(go.Scatter(x=t, y=sint*ex,
                         name=r'$e^t\sin t$', mode='lines'))  # 線のみ
fig.add_trace(go.Scatter(x=t, y=1/t, name=r'$1/t$',
                         mode='lines+markers'))  # ドットと線
# 軸の名前を追加
fig.update_xaxes(title="t")
fig.update_yaxes(title=r"$f(t)$")

# 動かせるグラフをブラウザーで表示
fig.show()
#htmlで保存
fig.write_html("figure.html")
#pngで保存
fig.write_image("figure.png")

Plotlyで簡単なグラフを描いてみる(Python)

Plotlyとは?

Pythonのグラフ描画ライブラリーで描画したグラフの拡大縮小,軸のスケールの変更などなどができる.
Plotlyのグラフ描画モジュールにはexpressとgraph_objectsがあるらしいが,基本はgraph_objectsでexpressモジュールはサクッとグラフを書くためのモジュールで内部ではgraph_objectsの関数を呼んでいるらしい.(公式ドキュメントPlotly Express | Python | Plotlyの雑和訳)

f:id:Ibetoge:20210427124756g:plain
Plotlyの動くグラフ

今回はexpressでサインカーブを書いてみる.
graph_objectsでの描写はこちら→
ibetoge.hatenablog.com

plotly.expressでとりあえず描く

環境 macOSpython 3.7.6,Plotly 4.14.3
Plotlyがない人はターミナルで下のコマンドを叩いてインストールしてね.

pip install plotly

モジュールをインポートする.

import numpy as np
import plotly.express as px

まずはデータを用意する.今回は横軸 tの範囲を0\le t\le 4\piとしてサインカーブf(t)=\sin tを書くのでそれを用意する.

t = np.linspace(0, 4*np.pi, 100)
sint = np.sin(t)

次にplotly.expressモジュールで描画する.

#グラフのインスタンスを作成
fig = px.line(x=t, y=sint)
#グラフをブラウザーで表示
fig.show()

これで動かせるグラフが描画できる.
グラフの軸に名前を設定するには,インスタンス作成時に

fig = px.line(x=t, y=sint, labels={'x': r"$t$", 'y': r"$\sin t$"})

とすればOK.ここでr"$\sin t$"とかのr"$$"の部分ですが,Plotlyで文字列にLatex数式を使うときはr"$$"を使用します.
もっとしっかり目盛りなども設定したい人向け→

f:id:Ibetoge:20210427144425p:plain
Expressで書いたサインカーブ(これはPNGなので動きません)

Plotlyのグラフをファイルに保存する

fig.show()

を用いると,実行した際に動かせるグラフがブラウザーに表示されるが,グラフをhtmlファイルとして保存したい場合は,

#htmlで保存
fig.write_html("figure.html")

でOK.画像ファイルとして保存するには,kaleidoパッケージが必要なので,次のどちらかのコマンドを叩いてkaleidoをインストールする.(kaleidoの代わりにorcaパッケージでも良いらしいが非推奨)
pipの場合

 pip install -U kaleido

anacondaの場合

 conda install -c conda-forge python-kaleido

kaleidoをインストールしたら,

#pngで保存
fig.write_image("figure.png")
#jpegで保存
fig.write_image("figure.jpeg")

と書けばすれば良し.こうすれば保存したhtmlファイルをブラウザで開けばいつでも動かせるグラフが見られる.もちろん画像ファイルとして保存した場合は動かせない.

最後にソースコード

import plotly.express as px
import numpy as np

# 横軸は0<t<4pi,グラフはsin(t)
t = np.linspace(0, 4*np.pi, 100)
sint = np.sin(t)

# グラフのインスタンス作成(軸ラベル付き)
fig = px.line(x=t, y=sint, labels={'x': r"$t$", 'y': r"$\sin t$"})
fig.update_xaxes(title="t")
fig.update_yaxes(title=r"$\sin t$")

# グラフをブラウザーに表示
fig.show()
# htmlで保存
fig.write_html("figure.html")

極座標形式の2次元データをMatplotlibで描写する

極座標からデカルト座標に直して描写する.

f:id:Ibetoge:20201219144408p:plainf:id:Ibetoge:20201219144020p:plain
右の図を作ろう

数値解析を極座標で行った後,2次元の極座標形式データを描写したい時がありました.
例えば,
{\displaystyle 
\begin{eqnarray}
  \left\{
    \begin{array}{l}
      r=x^2+y^2 \\
     \theta=\arccos\left(\frac{x}{\sqrt{x^2+y^2}}\right)
    \end{array}
  \right.
\end{eqnarray}
}
のような座標でのデータを描写したい時があります.つまり

double r[imax];
double theta[jmax];
double T[jmax][imax];

のような配列データを格納したcsvファイルやらが手元にあります.それらをpythonを用いてmatplotlibで描写しましょう.
(r,~\theta方向にそれぞれ{\rm imax,~jmax}分割してるような空間を配列r,thetaによって表現していて,
そのグリッド点で物理量T[jmax][imax]を定義しています.)
上の配列はそれぞれr.csv,theta.csv,T.csvに入っていることにしておきます.
今回は
{\displaystyle 
\begin{equation}
T=\mathrm{e}^{-r^2}
\end{equation}
}
として,Tを描写します.

流れとしては,

  1. 極座標配列をcsvから読み込む
  2. x-y平面のデータに変換(スプライン補完)
  3. 描写

という感じです.
まずは必要なモジュールたちをインポートします.

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

scipyはスプライン補完に使い,mpl_toolkits.mplot3dと matplotlib.pyplotは描写,numpyは色々なことに使います.

次にnumpyでcsvファイルを読み込んで配列に格納しましょう.

r = np.loadtxt("r.csv", delimiter=',')
theta = np.loadtxt("theta.csv", delimiter=',')
T_polar = np.loadtxt("T.csv", delimiter=',')

描写したいx-y平面とx-y平面でのTを定義します.今回は配列r,thetaで表せる部分をすべて含めて,グリッド幅は配列rと同じにしておきます.

size_r = r.shape[0]
x = np.linspace(-r[-1], r[-1], 2*size_r)
y = np.linspace(-r[-1], r[-1], 2*size_r)
T_cartesian = np.empty((2*size_r, 2*size_r))

スプライン補完でx-y平面でのTを計算

spl_T = interpolate.interp2d(theta, r, T_polar, kind='cubic')
for ix in range(0, 2*size_r):
    for iy in range(0, 2*size_r):
        r_tmp = (x[ix]**2+y[iy]**2)**0.5
        costh = x[ix]/(x[ix]**2+y[iy]**2)**0.5
        th_tmp = np.arccos(costh)
        T_cartesian[ix, iy] = spl_T(th_tmp, r_tmp)

最後に描写

x, y = np.meshgrid(x, y, indexing='ij')
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.contourf(x, y, T_cartesian)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

これで完成.

f:id:Ibetoge:20201219144408p:plainf:id:Ibetoge:20201219144020p:plain
左図:r-theta平面,右図:x-y平面

割と一般化できるやり方なので,これで任意の座標系から任意の座標系へ変換して描写できるようになります.

ちなみに今回はカラーマップの等値線図を書きましたが

ax = fig.add_subplot(1, 1, 1)
ax.contourf(x, y, T_cartesian)

の部分を変えるといろいろな図にできます.

グラフの種類

等値線図
ax = fig.add_subplot(1, 1, 1)
ax.contour(x, y, T_cartesian)
f:id:Ibetoge:20201219150159p:plain
等値線図
3次元ワイヤーフレーム
ax = fig.add_subplot(1, 1, 1, projection="3d")
ax.plot_wireframe(x, y, T_cartesian)
f:id:Ibetoge:20201219150408p:plain
ワイヤー
3次元曲面
ax = fig.add_subplot(1, 1, 1, projection="3d")
ax.plot_surface(x, y, T_cartesian)
f:id:Ibetoge:20201219150501p:plain
三次元曲面

コード

最後にcsv作るように書いたコードとプロットするのに使ったコードを載せときます.

CSV作るよう
#include <cmath>
#include <fstream>
#include <iostream>

using namespace std;

int main() {
  int imax = 200;
  int jmax = 128;
  double r[imax], theta[jmax], T[jmax][imax];
  double pi;
  pi = acos(-1.0);

  for (int i = 0; i < imax; i++) {
    r[i] = (i + 0.50) / imax;
  }
  for (int j = 0; j < jmax; j++) {
    theta[j] = (j + 0.50) / jmax * 2.0 * pi;
  }

  for (int i = 0; i < imax; i++) {
    for (int j = 0; j < jmax; j++) {
      T[j][i] = exp(-1.0 * pow(r[i], 2));
    }
  }
  ofstream ofsr("r.csv");
  for (int i = 0; i < imax - 1; i++) {
    ofsr << r[i] << ",";
  }
  ofsr << r[imax - 1] << endl;

  ofstream ofst("theta.csv");
  for (int j = 0; j < jmax - 1; j++) {
    ofst << theta[j] << ",";
  }
  ofst << theta[jmax - 1] << endl;

  ofstream ofsv("T.csv");
  for (int i = 0; i < imax; i++) {
    for (int j = 0; j < jmax; j++) {
      if (j == jmax - 1) {
        ofsv << T[j][i] << endl;
      } else {
        ofsv << T[j][i] << ",";
      }
    }
  }
}
プロットする
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

# csvを読み込む
r = np.loadtxt("r.csv", delimiter=',')
theta = np.loadtxt("theta.csv", delimiter=',')
T_polar = np.loadtxt("T.csv", delimiter=',')

# x-y平面を定義してスプライン補完
size_r = r.shape[0]
x = np.linspace(-r[-1], r[-1], 2*size_r)
y = np.linspace(-r[-1], r[-1], 2*size_r)
T_cartesian = np.empty((2*size_r, 2*size_r))

spl_T = interpolate.interp2d(theta, r, T_polar, kind='cubic')
for ix in range(0, 2*size_r):
    for iy in range(0, 2*size_r):
        r_tmp = (x[ix]**2+y[iy]**2)**0.5
        costh = x[ix]/(x[ix]**2+y[iy]**2)**0.5
        th_tmp = np.arccos(costh)
        T_cartesian[ix, iy] = spl_T(th_tmp, r_tmp)

# matplotlibで描写する
x, y = np.meshgrid(x, y, indexing='ij')
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.contourf(x, y, T_cartesian)
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('cartesian.png')

# 極座標での描写
theta, r = np.meshgrid(theta, r, indexing='ij')
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.contourf(theta, r, T_polar)
plt.xlabel(r'$\theta$')
plt.ylabel('r')
plt.savefig('polar.png')