interpolación trigonométrica

Si el número de puntos N es impar, digamos N=2K+1 , al aplicar la fórmula de Lagrange para la interpolación polinómica a la formulación polinómica en el plano complejo se obtiene que la solución se puede escribir en la forma

dónde

${\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}}}}.}$

El factor ${\displaystyle e^{-iKx+iKx_{k}}}$en esta fórmula se compensa el hecho de que la formulación del plano complejo también contiene potencias negativas de${\displaystyle e^{ix}}$y por lo tanto no es una expresión polinómica en${\displaystyle e^{ix}}$La corrección de esta expresión se puede verificar fácilmente observando que${\displaystyle t_{k}(x_{k})=1}$y eso${\displaystyle t_{k}(x)}$es una combinación lineal de las potencias derechas de${\displaystyle e^{ix}}$. Al utilizar la identidad

el coeficiente${\displaystyle t_{k}(x)}$puede escribirse en la forma

Si el número de puntos N es par, digamos N=2K , al aplicar la fórmula de Lagrange para la interpolación polinómica a la formulación polinómica en el plano complejo se obtiene que la solución se puede escribir en la forma

dónde

Aquí, las constantes${\displaystyle \alpha _{k}}$se puede elegir libremente. Esto se debe a que la función de interpolación ( 1 ) contiene un número impar de constantes desconocidas. Una elección común es exigir que la frecuencia más alta tenga la forma de una constante por${\displaystyle \cos(Kx)}$, es decir el${\displaystyle \sin(Kx)}$El término desaparece, pero en general se puede elegir la fase de la frecuencia más alta.${\displaystyle \varphi _{K}}$Para obtener una expresión para${\displaystyle \alpha _{k}}$, obtenemos usando ( 2 ) que ( 3 ) se puede escribir en la forma

${\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})}}.}$

Esto produce

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

y

${\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})}}.}$

Tenga en cuenta que debe tenerse cuidado para evitar infinitos causados ​​por ceros en los denominadores.

Simplificar aún más utilizando ( 4 ) sería un enfoque obvio, pero obviamente es complejo. Un enfoque mucho más sencillo es considerar el núcleo de Dirichlet.

${\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}},}$

dónde${\displaystyle N>0}$es extraño. Se puede ver fácilmente que${\displaystyle D(x,N)}$es una combinación lineal de las potencias derechas de${\displaystyle e^{ix}}$y satisface

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

Dado que estas dos propiedades definen de forma única los coeficientes${\displaystyle t_{k}(x)}$En ( 5 ), se deduce que

${\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}}}$

Aquí, la función sinc evita cualquier singularidad y se define por

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

Para${\displaystyle N}$Incluso, definimos el núcleo de Dirichlet como

${\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}}.}$

Una vez más, se puede ver fácilmente que${\displaystyle D(x,N)}$es una combinación lineal de las potencias derechas de${\displaystyle e^{ix}}$, no contiene el término${\displaystyle \sin {\tfrac {1}{2}}Nx}$y satisface

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

Utilizando estas propiedades, se deduce que los coeficientes${\displaystyle t_{k}(x)}$en ( 6 ) están dados por

${\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}}}$

Tenga en cuenta que${\displaystyle t_{k}(x)}$no contiene el${\displaystyle \sin {\tfrac {1}{2}}Nx}$también. Finalmente, tenga en cuenta que la función${\displaystyle \sin {\tfrac {1}{2}}Nx}$desaparece en todos los puntos${\displaystyle x_{m}}$Por lo tanto, siempre se pueden sumar múltiplos de este término, pero comúnmente se omite.

Una implementación en MATLAB de lo anterior se puede encontrar aquí y se presenta de la siguiente manera:

función P = triginterp ( xi,x,y ) % TRIGINTERP Interpolación trigonométrica. % Entrada: % xi puntos de evaluación para el interpolante (vector) % x nodos de interpolación equidistantes (vector, longitud N) % y valores de interpolación (vector, longitud N) % Salida: % P valores del interpolante trigonométrico (vector) N = longitud ( x ); % Ajustar el espaciado de la variable independiente dada. h = 2 / N ; escala = ( x ( 2 ) - x ( 1 )) / h ; x = x / escala ; xi = xi / escala ; % Evaluar el interpolante. P = zeros ( tamaño ( xi )); para k = 1 : N P = P + y ( k ) * trigcardinal ( xi - x ( k ), N ); finfunción tau = trigcardinal ( x,N ) ws = advertencia ( 'off' , 'MATLAB:divideByZero' ); % La forma es diferente para N par e impar. if rem ( N , 2 ) == 1 % impar tau = sin ( N * pi * x / 2 ) ./ ( N * sin ( pi * x / 2 )); else % par tau = sin ( N * pi * x / 2 ) ./ ( N * tan ( pi * x / 2 )); end advertencia ( ws ) tau ( x == 0 ) = 1 ; % valor fijo en x=0

El caso especial en el que los puntos x _n están igualmente espaciados es especialmente importante. En este caso, tenemos

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

La transformación que asigna los puntos de datos y _n a los coeficientes a _k , b _k se obtiene a partir de la transformada discreta de Fourier (DFT) de orden  N.

${\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}\,}$

(Debido a la forma en que se formuló el problema anteriormente, nos hemos limitado a un número impar de puntos. Esto no es estrictamente necesario; para un número par de puntos, se incluye otro término coseno correspondiente a la frecuencia de Nyquist ).

El caso de la interpolación solo de coseno para puntos igualmente espaciados, que corresponde a una interpolación trigonométrica cuando los puntos tienen simetría par , fue tratado por Alexis Clairaut en 1754. En este caso, la solución es equivalente a una transformada discreta de coseno . La expansión solo de seno para puntos igualmente espaciados, que corresponde a simetría impar, fue resuelta por Joseph Louis Lagrange en 1762, para la cual la solución es una transformada discreta de seno . El polinomio de interpolación completo de coseno y seno, que da origen a la DFT, fue resuelto por Carl Friedrich Gauss en un trabajo inédito alrededor de 1805, momento en el que también derivó un algoritmo de transformada rápida de Fourier para evaluarlo rápidamente. Clairaut, Lagrange y Gauss se ocuparon del estudio del problema de inferir la órbita de planetas , asteroides , etc., a partir de un conjunto finito de puntos de observación; dado que las órbitas son periódicas, una interpolación trigonométrica fue una elección natural. Véase también Heideman et al. (1984).