In [1]:
from sympy import *
import numpy as np

import matplotlib.pyplot as plt

def L2_inner_product(f, g, tau):
    # you can define t here to not rely on a global definition
    t = symbols('t', real=True)
    return (integrate(f*conjugate(g), (t, 0, tau))/tau).simplify()

def L2_norm(f, tau):
    return sqrt(L2_inner_product(f, f, tau))

def L2_coefficient(f, g, tau):
    return (L2_inner_product(f, g, tau)/L2_norm(g, tau)**2).simplify()

Function to Approximate

In [2]:
t = symbols('t', real=True)
f1 = Piecewise(
    (1, t < pi),   # t is zero before time pi
    (-1, True)    # otherwise, t is 1
)
plot(f1, (t, -2*pi, 4*pi))
Out[2]:
<sympy.plotting.plot.Plot at 0x7f81f47770a0>
In [3]:
f2 = Piecewise(
    (0, t < pi/2),   # t is zero before time pi
    (1, t < pi),      # otherwise, t is 1
    (0, True)      # otherwise, t is 1
)
plot(f2, (t, -pi, 4*pi))
Out[3]:
<sympy.plotting.plot.Plot at 0x7f817171f490>
In [4]:
f = f1 # choose function to evaluate

Fourier Series using $\{e^{-ikt}\}$, $k=[-\infty, \dots, \infty]$ as Basis

In [5]:
k = symbols('k', real=True, integer=True)
psi_k = exp(-I*k*t)
tau = 2*pi

ak = L2_coefficient(f, psi_k, tau)
In [6]:
# this is just a fancy way to make sympy print the left hand side of the equation "a_k = "
# otherwise if will just print the right hand side

Eq(symbols('\psi_k'), psi_k)
Out[6]:
$\displaystyle \psi_{k} = e^{- i k t}$
In [7]:
Eq(symbols('a_k'), ak)
Out[7]:
$\displaystyle a_{k} = \begin{cases} \frac{i \left(1 - \left(-1\right)^{k}\right)}{\pi k} & \text{for}\: k \neq 0 \\0 & \text{otherwise} \end{cases}$

$approx = a_0*\psi_0 + a_1*\psi_1 + a_{-1}*\psi_{-1} + \ldots$

In [8]:
n = 5  # onto basis functions psi_{-10} ... psi_{10}
start = -2*pi
stop = 4*pi

approx = 0
for k_i in range(-n, n+1):
    approx += (ak*psi_k).subs(k, k_i)

Eq(symbols('approx'), approx)
Out[8]:
$\displaystyle approx = - \frac{2 i e^{5 i t}}{5 \pi} - \frac{2 i e^{3 i t}}{3 \pi} - \frac{2 i e^{i t}}{\pi} + \frac{2 i e^{- i t}}{\pi} + \frac{2 i e^{- 3 i t}}{3 \pi} + \frac{2 i e^{- 5 i t}}{5 \pi}$
In [9]:
fig = plot(approx, f, (t, start, stop), legend=True, show=False)
fig[0].line_color='r'
fig[0].label = 'p_' + str(n)
fig[1].label = 'f'
fig.show()
In [10]:
def do_plot():
    # how to convert symbolic functions into numerical functions
    f_func = lambdify(t, f)
    approx_func = lambdify(t, approx)
    
    # generate data
    t_vals = np.linspace(-0.5*np.pi, 2.5*np.pi, 1000)
    f_vals = f_func(t_vals)
    
    approx_vals =  np.real(approx_func(t_vals))
    plt.plot(t_vals, f_vals)
    plt.plot(t_vals, approx_vals)
    
    
    #plt.plot(0*tau, 0.5, 'ro')
    #plt.plot(tau, 0.5, 'ro')

do_plot()
In [11]:
approx
Out[11]:
$\displaystyle - \frac{2 i e^{5 i t}}{5 \pi} - \frac{2 i e^{3 i t}}{3 \pi} - \frac{2 i e^{i t}}{\pi} + \frac{2 i e^{- i t}}{\pi} + \frac{2 i e^{- 3 i t}}{3 \pi} + \frac{2 i e^{- 5 i t}}{5 \pi}$
In [12]:
approx.expand(complex=True)
Out[12]:
$\displaystyle \frac{4 \sin{\left(t \right)}}{\pi} + \frac{4 \sin{\left(3 t \right)}}{3 \pi} + \frac{4 \sin{\left(5 t \right)}}{5 \pi}$

Power Spectrum

In [13]:
# there isn't a great way to do scatter plots in sympy, so we will plot our 
# power spectrum using the matplotlib library, imported at the top of this
# notebook
k_values = range(-n, n+1)
power = abs(ak)**2

# this is building a list of power values evaluating the function
# using list comprehension, a short way to write for loops
# https://www.programiz.com/python-programming/list-comprehension
power_values = [ power.subs(k, ki) for ki in k_values]
plt.bar(k_values, power_values)
plt.xlabel('k')
plt.ylabel('power: $|a_k|^2$')
plt.title('power spectrum')
plt.grid()
In [14]:
plt.bar(k_values, power_values)
Out[14]:
<BarContainer object of 11 artists>

Using $\{1, \cos(2\pi k t/\tau), \sin(2\pi k t/\tau\}$ as Basis

In [15]:
psi_1k = cos(k*t*2*pi/tau)
Eq(symbols('\psi_1k'), psi_1k)
Out[15]:
$\displaystyle \psi_{1k} = \cos{\left(k t \right)}$
In [16]:
alpha_k = L2_coefficient(f, psi_1k, tau)
Eq(symbols('alpha_k'), alpha_k)
Out[16]:
$\displaystyle \alpha_{k} = 0$
In [17]:
psi_2k = sin(k*t*2*pi/tau)
Eq(symbols('\psi_2k'), psi_2k)
Out[17]:
$\displaystyle \psi_{2k} = \sin{\left(k t \right)}$
In [18]:
beta_k = L2_coefficient(f, psi_2k, tau)
Eq(symbols('beta_k'), beta_k)
Out[18]:
$\displaystyle \beta_{k} = \begin{cases} \frac{2 \left(1 - \left(-1\right)^{k}\right)}{\pi k} & \text{for}\: k \neq 0 \\\text{NaN} & \text{otherwise} \end{cases}$
In [19]:
psi_1 = 1
a_0 = L2_coefficient(f, psi_1, tau)
Eq(symbols('a_0'), a_0)
Out[19]:
$\displaystyle a_{0} = 0$
In [20]:
approx = a_0
for k_i in range(1, n+1):
    approx += (alpha_k*psi_1k + beta_k*psi_2k).subs(k, k_i)
    print(approx)
4*sin(t)/pi
4*sin(t)/pi
4*sin(t)/pi + 4*sin(3*t)/(3*pi)
4*sin(t)/pi + 4*sin(3*t)/(3*pi)
4*sin(t)/pi + 4*sin(3*t)/(3*pi) + 4*sin(5*t)/(5*pi)
In [21]:
fig = plot(approx, f, (t, start, stop), legend=True, show=False)
fig[0].line_color='r'
fig[0].label = 'p_' + str(n)
fig[1].label = 'f'
fig.show()
In [22]:
approx
Out[22]:
$\displaystyle \frac{4 \sin{\left(t \right)}}{\pi} + \frac{4 \sin{\left(3 t \right)}}{3 \pi} + \frac{4 \sin{\left(5 t \right)}}{5 \pi}$