何かあれば GitHub のリポジトリに issue を作るか ryukau@gmail.com までお気軽にどうぞ。
Update: 2024-05-06
3 次補間の特性を比較します。ついでに線形補間も比較します。
ここでは補間したい位置の前後の数サンプルだけがわかっていれば使える方法のみを扱います。線形補間や Catmull-Rom 補間はこの種類の補間です。また等間隔でサンプリングされていることにします。等間隔でサンプリングされているときは補間の計算を簡略化できます。
補間したい曲線に関わるすべてのサンプルが必要になる方法としてはベジエ補間や
natural cubic
補間があげられます。ただし音の補間については、すべてのサンプルが事前にわかっているのであればスプライン補間よりも
scipy.signal.resample
のような離散フーリエ変換によるリサンプリングのほうが適しています。画像では離散フーリエ変換によるリサンプリングをつかうとリンギングが問題になるので、
PCHIP などの単調な補間のほうがいいかもしれません。
scipy.interpolate
などを参考にして以下の補間を実装しました。
以下のリンク先に Python3 の実装を掲載しています。
以下のリンク先に C++ の実装を掲載しています。
t は範囲が [0.0, 1.0] の補間係数です。
C++20 では std::lerp
が使えるようになります。
def linear(y0, y1, t):
return y0 + t * (y1 - y0)サンプルが等間隔で並んでいるときは Catmull-Rom 補間が使えます。以下の実装はウェーブテーブルの帯域制限と位相方向の補間で紹介したものです。
def cubic(y0, y1, y2, y3, t):
t2 = t * t
c0 = y1 - y2
c1 = (y2 - y0) * 0.5
c2 = c0 + c1
c3 = c0 + c2 + (y3 - y1) * 0.5
return c3 * t * t2 - (c2 + c3) * t2 + c1 * t + y1ディレイの実装で紹介したラグランジュ補間の実装をさらに整理しました。
def lagrange3(y0, y1, y2, y3, t):
u = 1 + t
d0 = y0 - y1
d1 = d0 - (y1 - y2)
d2 = d1 - ((y1 - y2) - (y2 - y3))
return y0 - u * (d0 + (1 - u) / 2 * (d1 + (2 - u) / 3 * d2))N 次のラグランジュ補間は以下の v
の式を展開すれば実装が求められます。 d の項の係数はパスカルの三角形と一致しそうです。
d(0, i) = y[i] - y[i + 1]
d(n, i) = d(n - 1, i) - d(n - 1, i + 1)
u = (N - 1) / 2. N is odd.
s(n) = d(n, n) + (n + 1 - u) / (n + 2) * s(n + 1)
s(N - 1) = d(N - 1, N - 1)
v = y[0] - u * s(0)Catmull-Rom との違いはオーバーシュートしないことです。 LFO の補間にはいいかもしれません。
def pchip(y0, y1, y2, y3, t):
m0 = y1 - y0
m1 = y2 - y1
m2 = y3 - y2
dk0 = 0.0 if m0 * m1 <= 0 else 2 * (m0 * m1) / (m0 + m1)
dk1 = 0.0 if m1 * m2 <= 0 else 2 * (m1 * m2) / (m1 + m2)
t2 = t * t
c0 = y1 - y2
c1 = dk0
c2 = c0 + c1
c3 = c0 + c2 + dk1
return c3 * t * t2 - (c2 + c3) * t2 + c1 * t + y1Akima 補間の特長は、 1 階微分しても滑らかでありながら、補間したい位置の前後の数サンプルだけがわかっていれば計算できることです。
def akima(y0, y1, y2, y3, y4, y5, t):
m0 = y1 - y0
m1 = y2 - y1
m2 = y3 - y2
m3 = y4 - y3
m4 = y5 - y4
w2 = abs(m1 - m0)
w3 = abs(m2 - m1)
w4 = abs(m3 - m2)
w5 = abs(m4 - m3)
b2 = m1 if w2 + w4 == 0 else (w4 * m1 + w2 * m2) / (w2 + w4)
b3 = m2 if w3 + w5 == 0 else (w5 * m2 + w3 * m3) / (w3 + w5)
c2 = 3.0 * m2 - 2.0 * b2 - b3
d2 = b2 + b3 - 2.0 * m2
return ((d2 * t + c2) * t + b2) * t + y2Uniform B-spline はコントロールポイントを通りません。音にかけるとローパスフィルタのようになります。
def uniformBSpline(y0, y1, y2, y3, t):
t2 = t * t
t3 = t2 * t
b0 = 1 / 6 - t / 2 + t2 / 2 - t3 / 6
b1 = 2 / 3 - t2 + t3 / 2
b2 = 1 / 6 + t / 2 + t2 / 2 - t3 / 2
b3 = t3 / 6
return b0 * y0 + b1 * y1 + b2 * y2 + b3 * y3Centripetal Catmull-Rom は Catmull-Rom のサンプリング間隔を
alpha で変えられるようにした補間方法です。 2
次元以上の曲線の描画でコントロールポイントの間隔が近いときに
alpha
を調整してループやふくらみなどの問題を回避することができます。
以下の実装は Wikipedia の記事に掲載されているコードをもとにして、 0 除算を避けるように改変しています。
def centripetalCatmullRom(y0, y1, y2, y3, t, alpha=0.5):
# 0 =< alpha =< 1.
if t <= 0:
return y1
if t >= 1:
return y2
t0 = 0
t1 = t0 + np.power(abs(y1 - y0), alpha)
if t0 == t1:
t1 += np.finfo(np.float64).eps
t2 = t1 + np.power(abs(y2 - y1), alpha)
if t1 == t2:
t2 += t1 * np.finfo(np.float64).eps
t3 = t2 + np.power(abs(y3 - y2), alpha)
if t2 == t3:
t3 += t2 * np.finfo(np.float64).eps
# ここに到達したときに t0 < t1 < t2 < t3 でなければ 0 除算が起きる。
t = t1 + t * (t2 - t1)
A1 = (t1 - t) / (t1 - t0) * y0 + (t - t0) / (t1 - t0) * y1
A2 = (t2 - t) / (t2 - t1) * y1 + (t - t1) / (t2 - t1) * y2
A3 = (t3 - t) / (t3 - t2) * y2 + (t - t2) / (t3 - t2) * y3
B1 = (t2 - t) / (t2 - t0) * A1 + (t - t0) / (t2 - t0) * A2
B2 = (t3 - t) / (t3 - t1) * A2 + (t - t1) / (t3 - t1) * A3
return (t2 - t) / (t2 - t1) * B1 + (t - t1) / (t2 - t1) * B2振幅、位相、群遅延特性をプロットして比較します。以下は周波数特性を計算するコード例です。実際にプロットに使ったコードは実装の節の Python3 の実装を参照してください。
import numpy as np
import scipy.signal as signal
def lagrange3(y0, y1, y2, y3, t):
u = 1 + t
d0 = y0 - y1
d1 = d0 - (y1 - y2)
d2 = d1 - ((y1 - y2) - (y2 - y3))
return y0 - u * (d0 + (1 - u) / 2 * (d1 + (2 - u) / 3 * d2))
impulseResponse = [
[lagrange3(*impulse, t) for impulse in np.eye(4)]
for t in np.linspace(0, 1, 9, endpoint=False)
]
sampleRate = 2
for fir in impulseResponse:
ω, response = signal.freqz(fir, fs=sampleRate) # response は複素数の周波数特性。
gain = 20 * np.log10(np.abs(response))
phase = np.unwrap(np.angle(response))
_, delay = signal.group_delay((fir, 1), fs=sampleRate)
# プロットは省略。以下では図の左上が振幅特性、左下が位相特性、右上が群遅延特性、右下がインパルス応答です。青くなるほど
t = 0 、黄色くなるほど t = 1
に近くなります。振幅特性は t = 0.5
を中心として対称性があるので重なっている曲線があります。
典型的な分数ディレイフィルタの特性がみられます。
0 < t < 1 のときの振幅特性はローパス。t = 0.5
のとき最もローパスが強くかかるが、群遅延はフラット。t = 0.5 を中心に ±0.5
の範囲で対称。
線形補間よりは低域でフラットな領域が増えています。単純に考えると FIR フィルタ係数が 2 倍長いので、振幅特性のロールオフが 2 倍ほど急峻になるはずです。
Catmull–Rom 補間とほとんど同じですが、低域でのフラットさが増して、高域での誤差が増えています。
Catmull–Rom 補間と 3 次ラグランジュ補間については周波数特性に基づいて選ぶよりも、ベンチマークをとって速いほうを使うことも考えられます。 godbolt.org でアセンブラ出力を見たところ、今回の実装では Catmull–Rom 補間のほうが 3 次ラグランジュ補間よりもコードが短くなる傾向がありました。
PCHIP
補間は入力信号によって補間の特性が変わるので、以下の周波数特性は参考程度のものです。ただし、線形補間、
Catmull–Rom 補間、 3
次ラグランジュ補間に比べると、周波数特性が悪くなるケースがあるということは言えそうです。また、群遅延特性が
t の値に応じて一様に分布しない点にも注意が必要です。
PCHIP 補間と同様に入力信号によって補間の特性が変わるので、以下の周波数特性は参考程度のものです。
ローパスが強くかかるので用途は限られそうです。
alpha = 0.5 としていますが、入力信号がインパルスのときは
0 除算を避けるコードパスを通るので、 alpha
の値を変えても同じ特性となります。音のデータは等間隔でサンプリングされていることが普通なので
centripetal Catmull-Rom
補間のようにサンプリング間隔を変えて補間の特性を変える手法は音の補間に適していないようです。エフェクタに歪みを加えるときにはいいかもしれません。
計算の速さが求められるときは線形補間で十分です。
3 次の補間を使うときは Catmull-Rom 補間か 3 次ラグランジュ補間が使えます。 Catmull–Rom 補間は 3 次ラグランジュ補間より高速ですが、低域での特性が悪くなります。
オーバーシュートが許容できないときは PCHIP 補間が使えます。
Akima 補間、 uniform B-spline 補間、 centripetal Catmull-Rom 補間は音の補間には不適です。音に癖をつけたいときはいいかもしれません。
以下のように sinc 関数から直接 FIR
フィルタ係数を計算する方法も考えられますが、 補間のたびに
sin() を FIR
のタップ数だけ繰り返し計算するのはリアルタイムの音の処理ではかなり重たいです。
def sinc4(y0, y1, y2, y3, t):
s0 = np.sinc(t + 1) / np.e
s1 = np.sinc(t)
s2 = np.sinc(t - 1)
s3 = np.sinc(t - 2) / np.e
denom = s0 + s1 + s2 + s3
return (y0 * s0 + y1 * s1 + y2 * s2 + y3 * s3) / denomnp.e
で値を割っているのは窓関数の計算です。以下のコードから出力されるインデックス
0 とインデックス 3 の値とほぼ一致します。
win = scipy.signal.get_window(("kaiser", 2 * np.pi), 6, fftbins=False)[1:-1]
print(win / win[1])特性は窓関数によって変わりますが、上のコードの sinc4 は
3 次ラグランジュ補間とほぼ一致します。
