# coding: utf-8

# # Fourier Interpolation

# In[2]:

import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt


# Let's fix the number of points and the number of functions as `n`.
# 
# Make sure `n` is odd.

# In[49]:

n = 5

assert n % 2 == 1


# In[50]:

x = np.linspace(0, 2*np.pi, n, endpoint=False)


# Next, fix the values of $k$ in $\cos(kx)$ as `cos_k` and in $sin(kx)$ as `sin_k` so that there are exactly $n$ altogether:

# In[51]:

cos_k = np.arange(0, n//2 + 1, dtype=np.float64)
sin_k = np.arange(1, n//2 + 1, dtype=np.float64)


# In[52]:

print(cos_k)
print(sin_k)


# Next, build the generalized Vandermonde matrix.
# 
# Make sure to order the matrix by increasing $k$:

# In[53]:

V = np.zeros((n,n))
V[:, ::2] = np.cos(cos_k*x[:, np.newaxis])
V[:, 1::2] = np.sin(sin_k*x[:, np.newaxis])
V


# Now let's try and do interpolation with this. Here are a few functions:

# In[54]:

if 1:
    def f(x):
        return x
elif 0:
    def f(x):
        return np.abs(x-np.pi)
elif 1:
    def f(x):
        return (x<=np.pi).astype(np.int32).astype(np.float64)


# Find the coefficients as `coeffs`:

# In[55]:

coeffs = la.solve(V, f(x))


# In[56]:

plot_x = np.linspace(0, 2*np.pi, 1000)


# In[57]:

interpolant = 0 * plot_x
for i, k in enumerate(cos_k):
    interpolant += coeffs[2*i] * np.cos(k * plot_x)
for i, n in enumerate(sin_k):
    interpolant += coeffs[2*i+1] * np.sin(n * plot_x)


# In[58]:

pt.plot(plot_x, interpolant)
pt.plot(plot_x, f(plot_x), "--", color="gray")
pt.plot(x, f(x), "or")


# * For $f(x)=x$, why does it do what it does at the interval end?
# * What happens if we increase the number of points?
# * What if the function is actually periodic (e.g. the function with `abs` above?)?
# * What if the function has jumps?

# In[59]:

# Answers
#
# * Because we're interpolating with periodic functions--so the interpolant is forced to be periodic.
# * We observe a distinct "overshoot". This overshoot is called "Gibbs phenomenon".
# * Periodic: no Gibbs at interval end.
# * Gibbs can also happen in the middle of the interval.


# ## Computing the coefficients

# In[60]:

B = V.T.dot(V)


# In[61]:

B[np.abs(B)<1e-12] = 0


# In[62]:

B


# * What do you observe?
# * How could this be useful for computing the inverse?
# * What is the normalization?
# * This is a pretty special matrix. What is the cost of computing $Ax$ with it?
# * The so-called [Fast Fourier transform](https://en.wikipedia.org/wiki/Fast_Fourier_transform) (FFT) computes a version of $Ax$ (with complex exponentials) in $O(n\log n)$!

# In[19]:

# Answers:
# 
# * V's columns are orthogonal. (though not normalized)
# * The transpose of $V$ (with appropriate normalization) its inverse. This makes Fourier coefficients cheap to compute.
# * The normalization is n for the first entry, and n/2 for all the ones after that.
# * Computing Ax costs n**2 operations.