Interpolação Trigonométrica.

Consideramos o problema de interpolação no intervalo [0,2p[ com funções trigonométricas (senos e co-senos, ou mais simplesmente através da exponencial complexa). Dado um conjunto de n nós de interpolação t0, ... , tn-1 em [0,2p[ pretende-se interpolar os pontos (t0, f0), ... , (tn-1, fn-1).
Através da mudança de variável x = eit, este problema está relacionado com a interpolação polinomial

 q(t) = a0+ a1x + .... + an-1xn-1 = a0+ a1eit + .... + an-1ei(n-1)t

e atribuindo xk = exp(itk), obtemos um sistema através da matriz de Vandermonde:

é
ê
ê
ê
ê
ë
1  x0  ...   x0n-1
1  x1  ...   x1n-1
... ...
1  xn-1  ... xn-1n-1
ù
ú
ú
ú
ú
û
 =  é
ê
ê
ê
ê
ë
1  exp(it0)  ...   exp(i(n-1)t0)
1  exp(it1)  ...   exp(i(n-1)t1)
... ...
1  exp(itn-1)  ... exp(i(n-1)tn-1)
ù
ú
ú
ú
ú
û

Interessa-nos especialmente o caso de nós igualmente espaçados tk=2pk/n, com k=0, ..., n-1.
Nessa situação obtemos como elementos da matriz xkm = exp(imtk) = exp(2pikm/n) para k,m=0, ..., n-1. Em particular, a matriz de Vandermonde é simétrica, com a primeira linha e coluna igual a 1 (caso k=0 ou m=0). Obtemos assim um sistema
 

é
ê
ê
ê
ê
ë
1  1  ...   1
exp(2pi/n)  ...   exp(2pi(n-1)/n)
... ...
exp(2pi(n-1)/n)  ...  exp(2pi(n-1)2/n)
ù
ú
ú
ú
ú
û
é
ê
ê
ê
ê
ë
a0
a1
...
an-1
ù
ú
ú
ú
ú
û
 =  é
ê
ê
ê
ê
ë
f0
f1
...
fn-1
ù
ú
ú
ú
ú
û.

que abrevidamente escrevemos na forma Wa = f.



Proposição: Tem-se a igualdade W*W=nI ,
(em que W* é a matriz adjunta - conjugada e transposta de W).
demonstração:
Um elemento (k,m) da matriz W*W é dado por
 
[W*W]km =  n-1
S
r =0
 e -2p ikr/n  e2p irm/n  n-1
S
r =0
e2p ir(m-k)/n
No caso em que m=k (elementos da diagonal) obtemos o somatório de 1's que dá n.
No caso em que temos m e k diferentes obtemos zero, porque designando Z=e2p i (m-k)/n, Z é diferente de 1, e podemos obter
 
[W*W]km =  n-1
S
r =0
 e2p ir(m-k)/n  n-1
S
r =0
Zr  Zn-1
  Z -1 
 = 0
pois Zn = exp(2pi(m-k)n)=1.



Corolário: A matriz de Vandermonde W é invertível com W-1=(1/n)W*. Portanto, a solução do sistema é dada por a = (1/n)W*f .

Definição: Dado um vector v = (v0, ..., vn-1) designamos por Fv =W*v,
a Transformada de Fourier Discreta de v, ou seja, é o vector cujas componentes são dadas por

[Fv]m  n-1
S
r = 0 
 ve -2pi rm/n,     (m = 0, ..., n-1) 

Observação: A Inversa da Transformação de Fourier Discreta pode ser obtida de forma semelhante. Basta reparar que WFv = WW*v = nv.
Portanto considerando a inversa é dada por F-1u = (1/n) Wu,  o que dá uma expressão semelhante à anterior, excepto que o sinal na exponencial é positivo e há uma divisão por n. Alguns autores podem considerar a raiz de 1/n em ambas as expressões ou colocar a divisão por n na definição de Fv.

A noção de Transformação de Fourier Discreta (DFT) tem várias aplicações, e o caso da interpolação é apenas um aspecto. 
O cálculo da DFT envolve O(n2) operações, mas através de uma sucessiva decomposição é possível definir um algoritmo que reduz esse número de operações para O(n log(n)). Trata-se de um algoritmo de Cooley e Tukey (de 1965, apesar de ter sido usado 160 anos antes, num trabalho de Gauss sobre trajectórias de asteróides). Esse algoritmo, ou variantes, são denominados Transformação de Fourier Rápida (FFT).


Valores Reais - Interpolação com senos e co-senos:
Iremos simplificar os resultados anteriores, para o caso em que fk são reais e q(t) também.
No caso em que os valores de vr são reais verifica-se que [Fv]n-m são os conjugados de [Fv]m porque

[Fv]n-m  n-1 
S
r = 0 
 vr exp(-2pi r (n-m)/n ) =   n-1
S
r = 0
 vr exp(+2pi r m/n ) = [Fv]m*
já que exp(-2pi r) =1. O valor [Fv]0  é obviamente real (reduz-se à soma dos vr ).
De forma análoga, reparamos que para pontos igualmente espaçados t0 , t1 , ..., tn-1 da forma tk=2pk/n, obtemos

exp(i (n-m) tk) = exp(i(n-m) 2pk/n) = exp(-im 2pk/n) = exp(-i mtk).

Número ímpar de pontos:
No caso em que temos n=2p+1 pontos, o valor nos nós tk , é dado por

fk = q(tk) = a0 + a1exp(itk) + .... + apexp(iptk)+ ap+1exp(i(p+1)tk) + .... + a2pexp(2iptk)

Em termos do vector solução, vimos que a = Fv  é da forma a = ( a0 , a1, ..., ap , ap*, ..., a1* ).
Também vimos que exp(i(2p+1-m)tk) = exp(-imtk), logo

exp(i(p+1)tk) = exp(-iptk), ..., exp(2iptk) = exp(-itk) .

Portanto,
fk = q(tk) = a0 + a1exp(itk) + .... + apexp(iptk)+ ( ap*)exp(-iptk) + .... + (a1*)exp(-itk)
= a0 + [a1exp(itk) + (a1exp(itk))*] + .... + [apexp(iptk) + (apexp(iptk))*]
e como z + z* = 2 Re(z), fica
= a0 + 2 Re[ a1exp(itk) +.... +apexp(iptk) ].

A função Q(t) = a0+ 2 Re[ a1eit + .... + apeipt] é real e, como vimos, é interpoladora nos 2p+1 pontos, trata-se do "polinómio" trigonométrico de grau p que se escreve através de senos e co-senos da seguinte forma:

 Q(t) = a0+ 2 [ Re(a1)cos(t) + .... + Re(ap)cos(pt) ] - 2 [ Im(a1)sin(t)+ .... + Im(ap)sin(pt) ]

porque Re[ak exp(ikt) ] = Re(ak)cos(kt) - Im(ak)sin(kt) .

Assim, no caso real, é possível obter uma função real interpoladora em 2p+1 pontos usando a parte real obtida pela interpolação com p+1 funções base complexas 1, eit, ..., eipt, ou de forma equivalente,  usando 2p+1 funções base reais, 1, cos(t), ..., cos(pt), ..., sin(t), ..., sin(pt).
Em termos de cálculo, basta avaliar a primeiras p+1 componentes dada pela Transformação de Fourier Discreta  (a0 , a1, ..., ap).