Interpolação trigonométrica: Séries de Fourier, DFT e FFT¶

• Até agora usamos polinômios para aproximar e interpolar
• Outra possibilidade: senos e cossenos
• Extremamente versátil e poderoso
• Usado em eqs. diferenciais, análise matemática, funções aleatórias, etc, etc, etc ...

Polinômio trigonométrico¶

\begin{align} S_N(x) = \frac{a_0}{2} +& a_1 \cos\left(\frac{2\pi x}{L_0}\right) + a_2\cos\left(\frac{4\pi x}{L_0}\right) + \cdots + a_N\cos\left(\frac{2\pi n x}{L_0}\right) + \\ &b_1 \sin\left(\frac{2\pi x}{L_0}\right) + b_2\sin\left(\frac{4\pi x}{L_0}\right) + \cdots + b_N\sin\left(\frac{2\pi N x}{L_0}\right) \end{align}

Nesta equação, $T_0$ é o período da função.

O que temos é uma base $\phi_n(t)$ com $2N+1$ termos: \begin{align} \phi_0(x) &= \frac{1}{2}\\ \phi_n(x) &= \cos\left(\frac{2\pi n x}{L_0}\right) \qquad 0\le n \le N \\ \phi_{n+N}(x) &= \sin\left(\frac{2\pi n x}{L_0}\right) \qquad 1\le n \le N \end{align}

Como calcular os coeficientes $a_i$ e $b_i$???¶

Ortogonalidade de senos e cossenos

\begin{align} \int_{-L_0/2}^{L_0/2} \cos\left(\frac{2\pi n x}{L_0}\right) \cdot \cos\left(\frac{2\pi m x}{L_0}\right) \: dx &= \left\{\begin{matrix}0 \quad n\neq m\\ \frac{L_0}{2} \quad n=m\end{matrix}\right.\\ \int_{-L_0/2}^{L_0/2} \sin\left(\frac{2\pi n x}{L_0}\right) \cdot \sin\left(\frac{2\pi m x}{L_0}\right) \: dx &= \left\{\begin{matrix}0 \quad n\neq m\\ \frac{L_0}{2} \quad n=m\end{matrix}\right.\\ \int_{-L_0/2}^{L_0/2} \cos\left(\frac{2\pi n x}{L_0}\right) \cdot \sin\left(\frac{2\pi m x}{T_0}\right) \: dx &= 0\\ \end{align}

Como calcular os coeficientes?¶

\begin{align} a_n &= \frac{2}{L_0}\int_{-L_0/2}^{L_0/2} y(x) \cdot \cos\left(\frac{2\pi n x}{L_0}\right)\: dx \\ b_m &= \frac{2}{L_0}\int_{-L_0/2}^{L_0} y(x) \cdot \sin\left(\frac{2\pi m x}{L_0}\right)\: dx \\ \end{align}

Exemplo: $y(x) = |x|, -1 \le x \le 1$ ($L_0 = 2$)¶

\begin{align} a_0 &= 2\int_0^1 x\:dx \qquad a_0 = 1\\ a_n &= 2\int_0^1 \cos\left(\frac{2\pi n x}{2}\right)\:dx \qquad a_n = \frac{2}{\pi^2 n^2}\cdot\left[(-1)^n - 1\right] \\ b_n &= 0 \end{align}

Usando o SymPy¶

Vamos generalizar?¶

Imagine que queremos aproximar $u(x)$ usando uma base ortogonal $\phi_n(x)$ com $1 \le n \le N$ no domínio $a \le x \le b$: $$ u(x) \approx u^\delta(x) = \sum_{n=1}^N \hat{u}_n \phi_n(x) $$ O erro da aproximação vale: $$ \varepsilon(x) = u^\delta(x) - (x) \quad\longrightarrow\quad \varepsilon^2 = \int_a^b w(x)\left[ u^\delta(x) - u(x) \right]^2\:dx $$

Nesta equação, $w(x)$ é uma função de peso que agora vamos considerar constante e igual a 1.

Mínimos quadrados.¶

$$ \varepsilon^2 = f(\hat{u}_1, \ldots, \hat{u}_N) $$ Para minimizar o erro, fazemos como fizemos na aula anterior: $$ \frac{\partial \varepsilon^2}{\partial \hat{u}_k} = 0 \quad 1 \le k \le N $$

Com isso chegamos ao seguinte sistema de equações lineares: $$ \sum_{j=1}^N \hat{u}_j \int_a^b w(x)\phi_n(x)\phi_k(x)\:dx = \int_a^b w(x)\hat{u}_k u(x)\:dx \qquad 1\le k \le N $$

Base ortogonal:¶

$$ \int_a^b w(x) \phi_j(x) \phi_k(x) \:dx = \left\{\begin{matrix} 0 & \quad & se\: j\ne k \\ \gamma(k) & \quad & se\: j=k\end{matrix}\right. $$ Quando $\gamma(k) = 1$ para $1\le k \le N$, a base é ortonormal.

Para a função de peso $w(x)=1$ e a base formada com senos e cossenos, recuperamos as equações para $a_n$ e $b_n$ acima.

Sempre que possível, use bases ortogonais

Quando não for possível, ainda ajuda usar bases parcialmente ortogonais.

Mínimos quadrados discreto¶

Neste caso, com uma base $\phi_n(x)$ com $1\le n \le N$ e com a função $u(x)$ conhecida nos pontos $x_k$, $1 \le k \le Q$, podemos usar o método dos mínimos quadrados: $$ \varepsilon^2 = \sum_{k=1}^Q w_k\left[u^\delta(x_k) - u(x_k)\right]^2 = \sum_{k=1}^Q w_k\left[\left(\sum_{n=1}^N \hat{u}_n \phi_n(x)\right) - u_k\right]^2 $$

Chegamos ao sistema apresentado na aula anterior:

$$ \left( \begin{matrix} \sum_{k=1}^Q w_k\phi_1(x_k)\cdot\phi_1(x_k) & \cdots & \sum_{k=1}^Q w_k\phi_1(x_k)\cdot\phi_N(x_k) \\ \sum_{k=1}^Q w_k\phi_2(x_k)\cdot\phi_1(x_k) & \cdots & \sum_{k=1}^Q w_k\phi_2(x_k)\cdot\phi_N(x_k) \\ \vdots & \ddots & \vdots \\ \sum_{k=1}^Q w_k\phi_N(x_k)\cdot\phi_1(x_k) & \cdots & \sum_{k=1}^Q w_k\phi_N(x_k)\cdot\phi_N(x_k) \\ \end{matrix}\right) \cdot \left(\begin{matrix} \hat{u}_1 \\ \hat{u}_2 \\ \vdots \\ \hat{u}_N \end{matrix}\right) = \left(\begin{matrix} \sum_{k=1}^Q w_k u_k \phi_1(x_k) \\ \sum_{k=1}^Q w_k u_k \phi_2(x_k) \\ \vdots \\ \sum_{k=1}^Q w_k u_k \phi_N(x_k)\end{matrix}\right) $$

Aqui introduzimos a possibilidade de haver um peso $w_k$ para cada ponto.

Voltando para os polinômios trigonométricos...¶

Ortogonalidade discreta¶

Já vimos que a base dos polinômios trigonométricos é ortogonal, de modo que para $k\ne n$, $$ \int_{-L_0/2}^{L_0/2} \phi_n(x) \cdot \phi_k(x) \: dx = 0 $$ Mas podemos aproximar esta integral, chamando $\Delta x = L_0 / Q$ $$ \int_{-L_0/2}^{L_0/2} \phi_n(x) \cdot \phi_k(x) \: dx \approx \Delta x\cdot \sum_{i=0}^Q \phi_n(x_i)\cdot \phi_k(x_i) \approx 0 \qquad x_i = -\frac{L_0}{2} + i\cdot\Delta x\quad i = 0, \ldots, Q-1 $$

Ortogonalidade discreta¶

Admitindo que temos $Q$ pontos $x_j = j\cdot \Delta x$ para $0\le j\le Q$ e $\Delta x = L_0/Q$, demonstra-se que

$$ \sum_{j=0}^{Q-1} \cos\left(\frac{2\pi n x_j}{L_0}\right) = 0 \qquad \sum_{j=0}^{Q-1} \sin\left(\frac{2\pi n x_j}{L_0}\right) = 0 $$ quando $n$ é multiplo de $Q$. Por outro lado

$$ \sum_{j=0}^{Q-1} \left[\cos\left(\frac{2\pi n x_j}{L_0}\right)\right]^2 = \frac{Q}{2} \qquad \sum_{j=0}^{Q-1} \left[\sin\left(\frac{2\pi n x_j}{L_0}\right)\right]^2 = \frac{Q}{2} $$ quando $n$ não é múltiplo de $Q/2$

Ortogonalidade discreta¶

Lembrando que \begin{align} 2\cos\theta\cos\phi &= \cos(\theta-\phi) + \cos(\theta+\phi)\\ 2\cos\theta\sin\phi &= \sin(\theta+\phi) - \sin(\theta-\phi)\\ \end{align}

Chega-se à relação de ortogonalidade desejada: $$ \sum_{i=0}^Q \phi_n(x_i)\cdot \phi_k(x_i)=0 \qquad k\ne n,k+n \text{ e } k-n \text{ não são múltiplos de}\: Q/2 $$

Transformada de Fourier real discreta¶

Com as relações de ortogonalidade, chega-se à seguinte relação para os coeficientes $a_n$ e $b_n$:

\begin{align} a_n &= \frac{2}{Q}\cdot\sum_{j=0}^{Q-1} y_j\cos\left(\frac{2\pi n x_j}{L_0}\right) \qquad n=0, 1, \ldots, N\\ b_n &= \frac{2}{Q}\cdot\sum_{j=0}^{Q-1} y_j\sin\left(\frac{2\pi n x_j}{L_0}\right) \qquad n=1, \ldots, N \end{align}

Outra abordagem para a transformada de Fourier real discreta¶

\begin{align} a_n &= \frac{2}{L_0}\int_{-L_0/2}^{L_0/2} y(x) \cdot \cos\left(\frac{2\pi n x}{L_0}\right)\: dx \\ b_m &= \frac{2}{L_0}\int_{-L_0/2}^{L_0} y(x) \cdot \sin\left(\frac{2\pi m x}{L_0}\right)\: dx \\ \end{align}

Aproximando as integrais: \begin{align} a_n &= \frac{2}{L_0}\int_{-L_0/2}^{L_0/2} y(x) \cdot \cos\left(\frac{2\pi n x}{L_0}\right)\: dx \approx \frac{2}{Q\cdot \Delta x} \sum_{j=0}^{Q-1} y(x_j) \cos\left(\frac{2\pi n x_j}{L_0}\right) \cdot \Delta x =\\ & = \frac{2}{Q} \sum_{j=0}^{Q-1} y_j \cos\left(\frac{2\pi n x_j}{L_0}\right) \end{align}

Onde $\Delta x = L_0 / Q$

E assim obtemos o mesmo resultado que no slide anterior!

Exemplo: $y(x) = |x|, -1 \le x \le 1$ ($L_0 = 2$)¶

Problemas digitais:¶

O que nós temos: $$ y_j, \Delta x, \quad\longrightarrow\quad x_{j+1}-x_j = \Delta x =\frac{L_0}{Q} \qquad 0\le j < Q $$ não temos os valores de $x_n$ exatos. Então podemos assumir $x_j = j \Delta x$: \begin{align} a_n &= \frac{2}{Q} \sum_{j=0}^{Q-1} y_j \cos\left(\frac{2\pi n x_j}{L_0}\right) = \frac{2}{Q} \sum_{j=0}^{Q-1} y_j \cos\left(\frac{2\pi n j \Delta x}{Q\cdot\Delta x}\right) = \frac{2}{Q} \sum_{j=0}^{Q-1} y_j \cos\left(\frac{2\pi n j}{Q}\right) \\ b_n &= \frac{2}{Q} \sum_{j=0}^{Q-1} y_j \sin\left(\frac{2\pi n x_j}{L_0}\right) = \frac{2}{Q} \sum_{j=0}^{Q-1} y_j \sin\left(\frac{2\pi n j}{Q}\right) \end{align}

Forma complexa da série de Fourier¶

$$ c_n = \left\{ \begin{matrix} \frac{a_0}{2} & n = 0 \\ \frac{1}{2}\left(a_n - i b_n\right) & n > 0\\ \frac{1}{2}\left(a_n + i b_n\right) & n < 0 \end{matrix}\right. $$ $$ f(x) = \sum_{n=-N}^N c_n \exp\left(\frac{2\pi i n x}{L_0} \right) \qquad c_n = \frac{1}{L_0}\int_0^{L_0} f(x) \exp\left(-\frac{2\pi i n x}{L_0} \right)\:dx $$ Se $f(x)$ for uma função real, para $n$ positivo: $$ \begin{matrix} c_n + c_{-n} = a_n\\ c_n - c_{-n} = -i b_n \end{matrix} \longrightarrow c_n = c^*_{-n} $$

DFT - Transformada de Fourier discreta¶

$$ c_n = \frac{1}{L_0}\int_0^{L_0} f(x) \exp\left(-\frac{2\pi i\cdot n x}{L_0} \right)\:dx \approx \frac{1}{\Delta x \cdot Q}\sum_{j=0}^{Q-1} y_j \exp\left(-\frac{2\pi i \cdot n j\Delta x}{Q\cdot \Delta x} \right)\cdot\Delta x $$

Portanto, $$ c_n = \frac{1}{Q}\sum_{j=0}^{Q-1} y_j \exp\left(-\frac{2\pi i\cdot n j}{Q} \right) $$

DFT - Discrete Fourier Transform¶

$$ c_n = X_n = \frac{1}{N}\sum_{k=0}^{N-1} x_k \exp\left(-\frac{2\pi n k}{N} \right) $$

Chamando $$ w_k = w^k = \exp\left(-\frac{2\pi k}{N}\right) $$ Repare que isto pode ser escrito como uma multiplicação de matrizes:

$$ \left(\begin{matrix} X_0 \\ X_1 \\ \vdots \\ X_{N-1}\end{matrix} \right) = \left(\begin{matrix} 1 & 1 & \cdots & 1\\ 1 & w & \cdots & w^{N-1} \\ \vdots & \vdots & \ddots &\vdots \\ 1 & w^N & \cdots & w^{(N-1)^2}\\ \end{matrix}\right) \cdot \left(\begin{matrix} x_0 \\ x_1 \\ \vdots \\ x_{N-1}\end{matrix} \right) $$

FFT - Transformada Rápida de Fourier¶

Um algorítmo rápido para $$ \sum_{j=0}^{Q-1} y_j \exp\left(-\frac{2\pi i\cdot n j}{Q} \right) $$

Se $0\le n\le Q$, o custo computacional é $\mathcal{O}(Q^2)$.

FFT: custo cai para $\mathcal{O}(Q\cdot\log Q)$

Originalmente, só funcionava com $Q$ igual a potências de 2.

FFTW - Fastest Fourier Transform of the West¶

Exemplos¶

A DFT é periódica (aliasing)!¶

Sejam os pontos $y_j$, $0 \le j < N$:

$$ Y_k = \frac{1}{N}\sum_{j=0}^{N-1} y_j\exp\left(-\frac{2\pi i j k}{N}\right) $$ $$ Y_{N+k} = \frac{1}{N}\sum_{j=0}^{N-1} y_j\exp\left[-\frac{2\pi i j (k+N)}{N}\right] = \frac{1}{N}\sum_{j=0}^{N-1} y_je^{-2\pi i j}\exp\left(-\frac{2\pi i j k}{N}\right) $$ $$ Y_{N+k} = \frac{1}{N}\sum_{j=0}^{N-1} y_j\exp\left(-\frac{2\pi i j k}{N}\right) = Y_k $$

Exemplo de Aliasing¶

O que acontece com frequências negativas (quando $y_j$ é real)?¶

E se o $y_j$ forem reais? $$ Y_{-k} = \frac{1}{N}\sum_{j=0}^{N-1} y_j\exp\left(\frac{2\pi i j k}{N}\right) = Y^*_k $$

Mas lembre-se que $Y$ também é periódico:

$$ Y_{-k} = Y_{N-k} = Y^*_k $$

Qual o número de frequencias relevantes?¶

• Q números reais
• fft gera Q números complexos - 2 $\times$ vezes mais informação
• A frequência 0 sempre vai corresponder a um número real (a média)!

Q é um número ímpar¶

• $a_0$
• $(Q-1)/2\:\longrightarrow\:a_n$
• $(Q-1)/2\:\longrightarrow\:b_n$

N é um número par¶

• $a_0$
• $Q/2\:\longrightarrow\:a_n$
• $Q/2\:\longrightarrow\:b_n$ mas $b_{1+Q/2} = 0$

$$ N_\text{relevante} = \frac{Q}{2} + 1 $$

rfft: especializado para números reais¶

Interpolação trigonométrica¶

Exemplo¶

E usando fft (não rfft)?¶

Transformação inversa¶

Resample¶

Usando myresample¶

Outros formatos¶

O mesmo para ifft, rfft, etc!

Lembra que a DFT era uma multiplicação de matrizes?¶

Transformada inversa¶

Outros tipos de transformadas¶

• Cosseno: dct, idct
• Mais genérica: r2r (ver <www.fftw.org>)
• Tudo o que foi feito até aqui, vale para problemas multidimensionais!

Conclusões¶

Algumas questões

• Será que esta série converge para a função original?
• Em que condições esta série converge?
• Quão rápido se dá esta convergência?
• O que acontece quando $N\longrightarrow\infty$?

Respostas genéricas

• A série converge para condições bem amplas, inclusive se houver discontinuidades na função.
• Esta convergência se dá para $N\longrightarrow\infty$. Para valor de $N$ finito, as coisas são mais complicadas (como veremos adiante).
• Para funções periódicas e suaves, não existe aproximação melhor!!!

Problemas¶

Problema 1¶

Seja a seguinte função definida no domínio $-1 < x \le 1$:

$$ y(x) = \left\{\begin{matrix} -1 & x < 0 \\ 1 & x > 0\\ 0 & x=0\end{matrix}\right. $$

1. Ache os coeficientes $a_n$ e $b_n$ da série de Fourier
2. Plote a série de Fourier para diferentes números de termos da série de Fourier. Percebe algo estranho?
3. Qual a relação desta função com a função $y = |x|$ que estudamos anteriormente?
4. Escolha diferentes números de pontos e interpole esta função nestes pontos

Problema 2¶

Verifique a ortogonalidade discreta da base trigonométrica, ou seja,

$$ \sum_{i=0}^Q \phi_n(x_i)\cdot \phi_k(x_i)=0 \qquad k\ne n, k+n \text{ e } k-n \text{ não são múltiplos de}\: Q/2 $$

para

\begin{align} \phi_0(x) &= \frac{1}{2}\\ \phi_n(x) &= \cos\left(\frac{2\pi n x}{L_0}\right) \qquad 0\le n \le N \\ \phi_{n+N}(x) &= \sin\left(\frac{2\pi n x}{L_0}\right) \qquad 1\le n \le N \end{align}

Problema 3¶

Nos exemplos iniciais usando a DFT, ao se escolher os pontos, foi utilizado o seguinte código:

L₀ = 1.0; Q = 16; freq=1
x = range(0.0, L₀, length=Q+1)[1:Q]

O que acontece se utilzarmos a expressão mais simples a seguir?

x = range(0.0, L₀, length=Q)

O que está acontecendo? Você consegue explicar?

Problema 4¶

Implemente a DFT usando a definição. Verifique como aumenta o custo computacional quando o número de pontos aumenta.

Problema 5¶

(Desafio) procure na net ou em algum livro como implementar a FFT. Tente!

Problema 6¶

A transformada de Fourier está diretamente relacionada com a série de Fourier e é dada por:

$$ X(\omega) = \frac{1}{2\pi} \int_{-\infty}^\infty x(t)e^{-i\omega t}\: dt $$

1. Calcula a transformada de Fourier para a função $$ x(t) = \left\{\begin{matrix} 1 & -1 < x < 1 \\ 0 & x \ge 1\end{matrix}\right. $$
2. Discretize a função e use a FFT para calcular a transformada.
3. Tente verificar o que acontece se você adiciona zeros à esquerda e à direita da função discretizada.

Problema 7¶

Faça alguns exemplos do pacote ApproxFun (https://github.com/JuliaApproximation/ApproxFun.jl)