sinc関数

sinc関数

正弦関数 $\sin x$ を $x$ で割った関数
 $$\begin{array}{}\text{(1)}& \mathrm{sinc}(x)=\frac{\sin x}{x}\end{array}$$
を sinc関数 とよびます。このままでは、$x=0$ で定義されない不連続な関数ですが、有名な極限値の公式
 $$\begin{array}{}\text{(2)}& \underset{x\to 0}{lim}\frac{\sin x}{x}=1\end{array}$$
を考慮して、
 $$\begin{array}{}\text{(3)}& \mathrm{sinc}(x)=\{\begin{array}{ll}\frac{\sin x}{x}& (x\ne 0)\\ 1& (x=0)\end{array}\end{array}$$
のように連続関数として定義するのが普通です。工学では全区間での積分が $1$ になるように正規化された
 $$\begin{array}{}\text{(4)}& \mathrm{sinc}(x)=\{\begin{array}{ll}\frac{\sin (\pi x)}{\pi x}& (x\ne 0)\\ 1& (x=0)\end{array}\end{array}$$
を使うのが一般的です。sinc関数は $x=0$ から正負の方向に減衰する振動関数です。
 
正規化された sinc関数は次の式で級数展開できます。
 $$\begin{array}{}\text{(5)}& \mathrm{sinc}(x)=\sum _{n=0}^{\mathrm{\infty }}\frac{(-1{)}^n{\pi }^{2\pi }}{(2n+1)!}{x}^{2n}\end{array}$$
後述する SymPy の sinc() は、正規化されていない sinc関数の任意の項までの級数展開を取得するメソッドを備えています。

NumPy には正規化された sinc 関数の値を返す numpy.sinc() が用意されています。$x$ が整数値をとるとき、
 $$\begin{array}{}\text{(6)}& \mathrm{sinc}(0)=1,\mathrm{sinc}(\pm 1)=\mathrm{sinc}(\pm 2)=\cdots =0\end{array}$$
となることを確認してみましょう。

# PYTHON_SINC

# In[1]

import numpy as np

# x=[-2, -1, 0, 1, 2]
x = np.arange(-2, 3)

# xが整数値をとるときのsinc(x)の値を確認
print(np.sinc(x))

# [-3.89817183e-17  3.89817183e-17  1.00000000e+00  3.89817183e-17  -3.89817183e-17]

実行結果を見ると、$x=0$ 以外の整数点では、”ほぼ” 0 になっています。数値計算の性質上、若干の誤差はあります。NumPy には正規化されていない sinc 関数は実装されていませんが、正規化された sinc 関数の変数 $x$ を $x/\pi$ に置き換えれば、正規化されていない sinc 関数が得られます。これらの関数が別々に定義されているのは煩わしいので、正規化・非正規化をオプションで選択できるような関数を定義しておくと便利です。

# In[2]

def sinc(x, norm=True):
    if norm == True:
        return np.sinc(x)
    else:
        return np.sinc(x/np.pi)

キーワード引数 norm を True に指定する、または省略すると、正規化された sinc 関数の値が返ります。norm にそれ以外の値を渡すと正規化されていない sinc 関数が返ります。Matplotlib をインポートして、それぞれの関数のグラフを同時に可視化してみます。

# In[3]

import matplotlib.pyplot as plt

# 区間[-15,15]を512分割した値を変数xに格納
x = np.linspace(-15, 15, 513)

# FigureとAxes
fig = plt.figure(figsize=(6.5, 5))
ax = fig.add_subplot(111)
ax.grid()
ax.set_xlabel("x", fontsize=15)
ax.set_ylabel("sinc(x)", fontsize=15)
ax.set_xlim(-15, 15)
ax.set_ylim(-0.4, 1.2)

# グラフをプロット
ax.plot(x, sinc(x), color="blue", label="normalized sinc function")
ax.plot(x, sinc(x, norm=False), color="red", label="unnormalized sinc function")

# 凡例を右上に表示
ax.legend(loc="upper right")

numpy.sinc()

numpy.sinc(x) は正規化された sinc関数 の値を返します。
以下のコードを実行すると、sinc関数の概形が表示されます。

# NUMPY_SINC_01

import numpy as np
import matplotlib.pyplot as plt

# 区間[-6,6]を512分割
x = np.linspace(-6, 6, 513)

# sinc(x)
sc = np.sinc(x)

# FigureとAxes
fig = plt.figure(figsize=(6.5, 5))
ax = fig.add_subplot(111)
ax.grid()
ax.set_xlabel("x", fontsize=15)
ax.set_ylabel("sinc(x)", fontsize=15)
ax.set_xlim(-6, 6)
ax.set_ylim(-0.4, 1.2)

# グラフをプロット
ax.plot(x, sc, color = "darkblue", label="sinc(x)")

# 凡例を表示
ax.legend(loc="upper left")

# 画像を保存
plt.savefig("sinc.png", bbox_inches = "tight")

sympy.functions.elementary.trigonometric.sinc()

SymPyパッケージから正規化されていない sinc関数を呼び出すことができます。

# SYMPY_SINC

# In[2]

from sympy import *

# 記号を定義
var("x")

# sinc(x)のグラフをプロット
plot(sinc(x), (x, -12, 12), ylabel = "sinc(x)")

$\mathrm{sinc}(x)$ を $9$ 次の項まで級数展開してみます。

# In[2]

# sinc関数を級数展開
sinc(x).series()

# 1 - x**2/6 + x**4/120 - x**6/5040 + x**8/362880 + O(x**10)

　sinc関数の不定積分が正弦積分 $\mathrm{Si}(x)$ になることを確認してみます。

# In[3]

# sinc関数の不定積分
integrate(sinc(x))

# Si(x)

sinc関数の積分 (正弦積分と余弦積分)

sinc関数の不定積分
 $$\begin{array}{}\text{(7)}& \mathrm{Si}(x)={\int }_0^x\frac{\sin t}{t}dt\end{array}$$
を正弦積分とよび、
 $$\begin{array}{}\text{(8)}& \mathrm{Ci}(x)=-{\int }_0^{x+\mathrm{\infty }}\frac{\cos t}{t}dt=\gamma +\log x+{\int }_0^x\frac{\cos t-1}{t}dt\end{array}$$
によって定義される関数を余弦積分といいます。正弦積分と余弦積分をまとめて三角積分とよびます。$\gamma$ はオイラーの定数で、およそ $0.57721$ の値をとります。

scipy.special.sici()

scipy.special.sici() は正弦積分と余弦積分の値を返します。

# SCIPY_SICI

import numpy as np
from scipy.special import sici
import matplotlib.pyplot as plt

# 区間[-6,6]を512分割
x = np.linspace(-10, 10, 513)

# 正弦積分
si = sici(x)[0]

# 余弦積分
ci = sici(x)[1]

# FigureとAxes
fig = plt.figure(figsize=(6.5, 5))
ax = fig.add_subplot(111)
ax.grid()
ax.set_xlabel("x", fontsize=15)
ax.set_ylabel("si(x), ci(x)", fontsize=15)
ax.set_xlim(-6, 6)
ax.set_ylim(-2, 2)

# グラフをプロット
ax.plot(x, si, color = "blue", label="si(x)")
ax.plot(x, ci, color = "red", label="ci(x)")

# 凡例を表示
ax.legend(loc="upper left")

# 画像を保存
plt.savefig("sici.png", bbox_inches = "tight")

【参考文献】https://ja.wikipedia.org/wiki/Sinc関数