三角関数補間

点の数Nが奇数、例えばN=2K+1の場合、複素平面上の多項式表現にラグランジュの多項式補間公式を適用すると、解は次の形式で表すことができます。

どこ

${\displaystyle t_{k}(x)=e^{-iKx+iKx_{k}}\prod _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K}{\frac {e^{ix}-e^{ix_{m}}}{e^{ix_{k}}-e^{ix_{m}}}}.}$

その要因 ${\displaystyle e^{-iKx+iKx_{k}}}$この式では、複素平面の定式化には負のべき乗も含まれているという事実を補正します。${\displaystyle e^{ix}}$したがって、多項式表現ではない${\displaystyle e^{ix}}$この表現の正しさは、以下の点を観察することで容易に検証できる。$t_k(x_k)=1$そしてそれは$t_k(x)$は、右のべき乗の線形結合です。${\displaystyle e^{ix}}$. IDを使用すると

係数$t_k(x)$次のように書くことができます

点の数Nが偶数、例えばN=2Kの場合、複素平面上の多項式表現にラグランジュの多項式補間公式を適用すると、解は次の形式で表すことができます。

どこ

ここで定数は$\displaystyle \alpha _{k}}$は自由に選択できます。これは、補間関数（ 1 ）に奇数個の未知定数が含まれているためです。一般的な選択は、最高周波数が定数倍の形であることを要求することです。${\displaystyle \cos(Kx)}$つまり、${\displaystyle \sin(Kx)}$項は消滅するが、一般に最高周波数の位相は次のように選択できる。$\displaystyle \varphi _{K}}$式を取得するには$\displaystyle \alpha _{k}}$（ ２）を用いると、 （ ３）は次の形式で書けることがわかる。

${\displaystyle t_{k}(x)={\frac {\cos {\tfrac {1}{2}}{\Biggl (}2Kx-\alpha _{k}+\displaystyle \sum \limits _{m=0,\,m\neq k}^{2K-1}x_{m}{\Biggr )}+\sum \limits _{m=-(K-1)}^{K-1}c_{k}e^{imx}}{2^{N}\sin {\tfrac {1}{2}}(x_{k}-\alpha _{k})\displaystyle \prod \limits _{m=0,\,m\neq k}^{2K-1}\sin {\tfrac {1}{2}}(x_{k}-x_{m})}}.}$

これにより、

${\displaystyle \alpha _{k}=\sum _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K-1}x_{m}-2\varphi _{K}}$

そして

${\displaystyle t_{k}(x)={\frac {\sin {\tfrac {1}{2}}(x-\alpha _{k})}{\sin {\tfrac {1}{2}}(x_{k}-\alpha _{k})}}\prod _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K-1}{\frac {\sin {\tfrac {1}{2}}(x-x_{m})}{\sin {\tfrac {1}{2}}(x_{k}-x_{m})}}.}$

分母にゼロが含まれることで無限大が生じるのを避けるため、注意が必要であることに注意してください。

（4）を用いることでさらに簡略化することは明白なアプローチだが、明らかに複雑である。より単純なアプローチは、ディリクレ核を考えることである。

${\displaystyle D(x,N)={\frac {1}{N}}+{\frac {2}{N}}\sum _{k=1}^{(N-1)/2}\cos(kx)={\frac {\sin {\tfrac {1}{2}}Nx}{N\sin {\tfrac {1}{2}}x}},}$

どこ${\displaystyle N>0}$奇妙だ。${\displaystyle D(x,N)}$は、右のべき乗の線形結合です。${\displaystyle e^{ix}}$そして満たす

${\displaystyle D(x_{m},N)={\begin{cases}0{\text{ for }}m\neq 0\\1{\text{ for }}m=0\end{cases}}.}$

これら2つの特性は係数を一意に定義するため${\displaystyle t_{k}(x)}$（ ５ ）において、以下が導かれる。

${\displaystyle {\begin{aligned}t_{k}(x)&=D(x-x_{k},N)={\begin{cases}{\dfrac {\sin {\tfrac {1}{2}}N(x-x_{k})}{N\sin {\tfrac {1}{2}}(x-x_{k})}}{\text{ for }}x\neq x_{k}\\[10mu]\lim \limits _{x\to 0}{\dfrac {\sin {\tfrac {1}{2}}Nx}{N\sin {\tfrac {1}{2}}x}}=1{\text{ for }}x=x_{k}\end{cases}}\\&={\frac {\mathrm {sinc} \,{\tfrac {1}{2}}N(x-x_{k})}{\mathrm {sinc} \,{\tfrac {1}{2}}(x-x_{k})}}.\end{aligned}}}$

ここで、sinc関数は特異点を防ぎ、次のように定義されます。

${\displaystyle \mathrm {sinc} \,x={\frac {\sin x}{x}}.}$

のために${\displaystyle N}$さらに、ディリクレ核を次のように定義します。

${\displaystyle D(x,N)={\frac {1}{N}}+{\frac {1}{N}}\cos {\tfrac {1}{2}}Nx+{\frac {2}{N}}\sum _{k=1}^{(N-1)/2}\cos(kx)={\frac {\sin {\tfrac {1}{2}}Nx}{N\tan {\tfrac {1}{2}}x}}.}$

また、容易にわかるように、${\displaystyle D(x,N)}$は、右のべき乗の線形結合です。${\displaystyle e^{ix}}$、という用語を含みません${\displaystyle \sin {\tfrac {1}{2}}Nx}$そして満たす

${\displaystyle D(x_{m},N)={\begin{cases}0{\text{ for }}m\neq 0\\1{\text{ for }}m=0\end{cases}}.}$

これらの性質を用いると、係数は次のようになる。${\displaystyle t_{k}(x)}$（６）は次のように与えられる。

${\displaystyle {\begin{aligned}t_{k}(x)&=D(x-x_{k},N)={\begin{cases}{\dfrac {\sin {\tfrac {1}{2}}N(x-x_{k})}{N\tan {\tfrac {1}{2}}(x-x_{k})}}{\text{ for }}x\neq x_{k}\\[10mu]\lim \limits _{x\to 0}{\dfrac {\sin {\tfrac {1}{2}}Nx}{N\tan {\tfrac {1}{2}}x}}=1{\text{ for }}x=x_{k}.\end{cases}}\\&={\frac {\mathrm {sinc} \,{\tfrac {1}{2}}N(x-x_{k})}{\mathrm {sinc} \,{\tfrac {1}{2}}(x-x_{k})}}\cos {\tfrac {1}{2}}(x-x_{k})\end{aligned}}}$

ご了承ください${\displaystyle t_{k}(x)}$含まれていません${\displaystyle \sin {\tfrac {1}{2}}Nx}$同様に。最後に、関数に注意してください。${\displaystyle \sin {\tfrac {1}{2}}Nx}$すべての点で消滅する${\displaystyle x_{m}}$したがって、この項の倍数は常に加算できますが、通常は省略されます。

上記のMATLAB実装はここにあり、以下のとおりです。

function P = triginterp ( xi,x,y ) % TRIGINTERP 三角関数による補間。% 入力: % xi 補間関数の評価点 (ベクトル) % x 等間隔の補間ノード (ベクトル、長さ N) % y 補間値 (ベクトル、長さ N) % 出力: % P 三角関数による補間値 (ベクトル) N = length ( x ); % 指定された独立変数の間隔を調整します。h = 2 / N ; scale = ( x ( 2 ) - x ( 1 )) / h ; x = x / scale ; xi = xi / scale ; % 補間関数を評価します。P = zeros ( size ( xi )); for k = 1 : N P = P + y ( k ) * trigcardinal ( xi - x ( k ), N ); endfunction tau = trigcardinal ( x,N ) ws = warning ( 'off' , 'MATLAB:divideByZero' ); % 偶数と奇数のNで形式が異なります。if rem ( N , 2 ) == 1 % 奇数tau = sin ( N * pi * x / 2 ) ./ ( N * sin ( pi * x / 2 )); else % 偶数tau = sin ( N * pi * x / 2 ) ./ ( N * tan ( pi * x / 2 )); end warning ( ws ) tau ( x == 0 ) = 1 ; % x=0 の値を固定

点x _nが等間隔に配置されている特殊なケースは特に重要です。この場合、次のようになります。

${\displaystyle x_{n}=2\pi {\frac {n}{N}},\qquad 0\leq n<N.}$

データ点y _nを係数a _k、b _kにマッピングする変換は、次数Nの離散フーリエ変換(DFT)から得られます。 

${\displaystyle Y_{k}=\sum _{n=0}^{N-1}y_{n}\ e^{-i2\pi nk/N}\,}$

${\displaystyle y_{n}=p(x_{n})={\frac {1}{N}}\sum _{k=0}^{N-1}Y_{k}\ e^{i2\pi nk/N}\,}$

（上記の問題の定式化の都合上、ここでは奇数個の点に限定して議論を進めています。これは必ずしも必要ではありません。偶数個の点の場合は、ナイキスト周波数に対応する別のコサイン項を含める必要があります。）

等間隔の点に対するコサインのみの補間（点が偶対称である場合は三角関数補間に相当）は、1754 年にAlexis Clairautによって扱われました。この場合、解は離散コサイン変換に相当します。奇対称に対応する等間隔の点に対するサインのみの展開は、 1762 年にJoseph Louis Lagrangeによって解決され、その解は離散サイン変換です。DFT を生み出す完全なコサインとサインの補間多項式は、1805 年頃の未発表の研究でCarl Friedrich Gaussによって解決され、その時点で彼はそれを迅速に評価するための高速フーリエ変換アルゴリズムも導出しました。Clairaut、Lagrange、Gauss は皆、有限個の観測点から惑星、小惑星などの軌道を推測する問題の研究に関心を持っていました。軌道は周期的であるため、三角関数による補間が自然な選択であった。Heideman et al. (1984)も参照のこと。