sin、cos を再帰的に求める / 加法定理による方法 & IIR フィルターによる方法 / Python / それに意味のないとき #算数

参考記事: 三上直樹,「知っ得! 軽量sin/cos計算アルゴリズム第4回その2: 高速に正弦波が生成できる再帰的計算法」,『インターフェース』, CQ出版, 2017年10月号, pp.158-159

sin(0 Δθ), sin(1 Δθ), sin(2 Δθ), sin(3 Δθ), ...... のように sin を順次求めることを考える。同時に cos も求める。

これは、リソースの乏しいマイコンで二相発振器などを構成するときに用いる方法であって、PC で計算する意味はない。手順の確認だけである。

加法定理による方法

計算手順は下のとおりである。

まず sin(0)、cos(0)、sin(Δθ)、cos(Δθ) の 4 つだけを組み込み函数か何かでちゃんと計算しておく。
最初に求めるべき sin(0 Δθ) は計算済みである。cos も同様である。
次に求めるべき sin(1 Δθ) も計算済みである。cos も同様である。
次に求めるべき sin(2 Δθ) は sin(Δθ + Δθ) であるので、加法定理により sin(Δθ) cos(Δθ) + cos (Δθ) sin(Δθ) と表現できる。各要素はすべて計算済みであるので掛け算と足し算だけで求まる。cos も同様である。
次に求めるべき sin(3 Δθ) は sin(2 Δθ + Δθ) であるので、加法定理により sin(2 Δθ) cos(Δθ) + cos (2 Δθ) sin(Δθ) と表現できる。各要素はすべて計算済みであるので掛け算と足し算だけで求まる。cos も同様である。以下同様に掛け算、足し算だけで sin、cos が順次求まる。

この方法は、誤差を含んだ sin、cos の値を基にして次の sin、cos の値を決定するので振幅がどんどんずれていってしまう。それを防ぐ手段として sin、cos の値を振幅 sqrt(sin^2 + cos^2) で割って強制的に振幅を 1 にすることにする。すなわち求まった sin、cos に 1/(sqrt(sin^2 + cos^2)) を乗じるということである。

ただし 1/sqrt(x) の計算自体に時間がかかる。しかし振幅がずれるといっても 1 から大きくずれることはない。そこで 1/sqrt(x) を、x=1 を中心にテイラー展開して、その 2 項目まで (1.5 - 0.5 * x) を代わりに使って振幅を補正することにする。

(見る/隠す) ソースコード (.py) と実行結果

import math
import matplotlib.pyplot as plt

def SinCos(Δθ):
    Δθ    = Δθ
    sin0  = math.sin(0)
    cos0  = math.cos(0)
    sinΔθ = math.sin(Δθ)
    cosΔθ = math.cos(Δθ)
    sin   = [None, sin0]
    cos   = [None, cos0]

    def sincos():
        sin[0] = sin[1]
        cos[0] = cos[1]

        sin[1] = sin[0] * cosΔθ + cos[0] * sinΔθ
        cos[1] = cos[0] * cosΔθ - sin[0] * sinΔθ

        adj = 1.5 - 0.5 * (sin[0]**2 + cos[0]**2)
        return sin[0] * adj, cos[0] * adj

    return sincos


########
# test #
########
Δθ      = math.pi/8
sin_cos = SinCos(Δθ) # 初期値をセットした状態でクロージャを開いて、
θList   = []
sinList = []
cosList = []
for i in range(50):
    sin, cos = sin_cos() # 外に出した函数を実行し続ける。実行するたびに sin、cos の値が更新される。

    θList.append(i * Δθ)
    sinList.append(sin)
    cosList.append(cos)
    #print("θ: %.2f, sin: %.6f, cos: %.6f, amp: %.16f" % (i*Δθ, sin, cos, math.sqrt(sin**2 + cos**2)))

plt.plot(θList,sinList,
         θList,cosList)
plt.xlabel("radian")
plt.ylabel("amplitude")
plt.grid(color="lightblue")
plt.show()

実行結果:

IIR フィルターによる方法

記事に紹介されている計算式を Python で再現しただけである。理窟は全然わからないが、ものすごく単純である。加法定理による方法と違って振幅は 1 からずれないとの由。

計算手順は下のとおりである。

a = 2 cos(Δθ)、c = -cos(Δθ)、s = sin(Δθ) の 3 つだけを組み込み函数か何かでちゃんと計算しておく。
u[n] = a * u[n-1] - u[n-2] という漸化式を考える。n は 0, 1, 2, ......。初期値は u[0] = 0、u[1] = 1 とする。

sin、cos は、この漸化式 u[n] からそれぞれ sin(n Δθ) = s * u[n]、cos(n Δθ) = u[n+1] + c * u[n] で求まる。n は 0, 1, 2, ......。

(見る/隠す) ソースコード (.py) と実行結果

import math
import matplotlib.pyplot as plt

def SinCos(Δθ):
    A = 2 * math.cos(Δθ)
    C = -math.cos(Δθ)
    S = math.sin(Δθ)
    u = [0, 1]
    def sincos():
        sin = S * u[0]
        cos = u[1] + C * u[0]
        u[0], u[1] = u[1], A * u[1] - u[0]
        return sin, cos

    return sincos

########
# test #
########
Δθ      = math.pi/8
sin_cos = SinCos(Δθ) # 初期値をセットした状態でクロージャを開いて、
θList   = []
sinList = []
cosList = []

for i in range(50):
    sin, cos = sin_cos() # 外に出した函数を実行し続ける。実行するたびに sin、cos の値が更新される。

    θList.append(i * Δθ)
    sinList.append(sin)
    cosList.append(cos)
    #print("θ: %.2f, sin: %.6f, cos: %.6f, amp: %.16f" % (i*Δθ, sin, cos, math.sqrt(sin**2 + cos**2)))


plt.plot(θList,sinList,
         θList,cosList)
plt.xlabel("radian")
plt.ylabel("amplitude")
plt.grid(color="lightblue")
plt.show()

実行結果: