Continuous Time Fourier Series
==============================

The CT Fourier Series considers periodic signals. A signal $x(t)$ is
periodic with period $T$ in case:

.. math::
   \forall t:\quad x(t+T) = x(t)

Note that in case $x$ is periodic with period $T$ it is also periodic
with period $2T$ (or $nT$). This leads us to the definition of the
**fundamental period** $T_0$ being the smallest value such that
$x(t+T_0)=x(t)$. The **fundamental frequency** then is $f_0=1/T_0$.


The Complex Exponential Functions
---------------------------------

We have already introduced the complex exponential functions:

.. math::
   x(t) = e^{j \omega t}

This funtion is periodic with period

.. math::
   T = \frac{2\pi}{\omega}

If we fix $\omega$, say at value $\omega_0$ then the fundamental
period is $T_0=2\pi/\omega_0$. The functions $x(t)=\exp(j k\omega_0
t)$ for integer $k>0$ then are all $T_0$ periodic, in fact the
fundamental period is:

.. math::
   T_k = \frac{2\pi}{k\omega_0} = \frac{T_0}{k}

Consider only the imaginary part of the complex exponential:

.. math::
   \sin(k \omega_0 t) = \sin\left(\frac{2 k \pi}{T_0} t\right)


.. ipython:: python
      
   from matplotlib.pylab import *
   T0 = 10
   t = linspace(0,T0,1000)
   for k in range(1,5):
       subplot(4,1,k)
       plot(t, sin(2*k*pi*t / T0))
   @savefig sines.png
   show()
   

The Fourier Series
------------------

Fourier (a French mathematician) showed that nearly every periodic
function (with period $T_0$) can be written as a linear combination
(superposition) of the complex exponential functions $\exp(j k
\omega_0 t)$ for $k=0,1,2,\cdots$ where $\omega_0 = 2\pi / T_0$:

.. math::
   x(t) = \sum_{k=-\infty}^{\infty} a_k e^{ j k \omega_0 t}

Note that the coefficients determining the weight of every frequency
$k\omega_0$ in general are complex valued.


..  As an example consider the
    periodic function $x(t)$ with period $T_0$:

    .. math::
       x(t) = \begin{cases}
       1 &: 0\leq t <T_0/2\\
       0 &: T_o/2 \leq t < T_0
       \end{cases}

    Note that we only defined the periodic function on one of its
    periods. The function has the entire real axis as its domain.

    This function can be written as the sum of complex exponentials with:

    .. math::
       a_0 &= 1/2\\
       a_k &= \frac{sin(\frac{\pi}{2}k)}{\pi k}\quad \text{for } k\geq 1


    .. ipython:: python
       :suppress:

       from matplotlib.pylab import *
       T0 = 1

       def f(t):
	   return ((t-floor(t))<=0.5).astype('float')


       t = linspace(-0.2,1.7,1000)
       x = f(t)

       def phi(t, k):
	   if k==0:
	       return 1/2*ones_like(t)
	   return sin(pi/2*k)/(pi*k) * exp(1j*k*2*pi/T0*t)

       K = 9
       xa = zeros_like(t)

       for k in range(1,K,2):
	   xak = phi(t,k)
	   xa = xa + xak
	   subplot(K,2,2*k+1)
	   plot(t,x,lw=1)
	   plot(t,xak,lw=2)
	   subplot(K,2,2*k+2)
	   plot(t,x,lw=1)
	   plot(t,xa,lw=2)

       axis([-0.2, 1.7, -0.1, 1.1])
       @savefig fourierblock.png
       show()


The **Fourier coefficients** are given by:

.. math::
   a_k = \frac{1}{T_0} \int_{\langle T_0 \rangle} x(t) e^{-j k \omega_0 t} dt


Summarizing:

.. list-table::
  :widths: 10 10
  :header-rows: 1
  :class: tablefullwidth

  * - Synthesis

    - Analysis

  * - 
      .. math::
	 x(t) = \sum_{k=-\infty}^{\infty} a_k e^{ j k \omega_0 t}

    -
      .. math::
	 a_k = \frac{1}{T_0} \int_{\langle T_0 \rangle} x(t) e^{-j k \omega_0 t} dt


	     


Properties I
------------

In signal processing we are most often working with real valued
signals. For such a **real valued signal** we have:

.. math::
   a_k = a^\star_{-k}

This is not too difficult to prove. First the synthesis equation can
be rewritten such that we 'unroll the summation loop':

.. math::
   x(t) = a_0 + \sum_{k=1}^{\infty}
   \left( a_{-k} e^{ - j k \omega_0 t} + a_k e^{ j k \omega_0 t} \right)

Using Eulers formula and writing $a_k=A_k+jB_k$ where $A_k$ and $B_k$
are real numbers, we get:

.. math::
   x(t) = a_0 + \sum_{k=1}^{\infty}
   \left( (A_{-k}+jB_{-k}) (\cos(k\omega_0 t) - j\sin(k\omega_0 t)) +
          (A_k+jB_k) (\cos(k\omega_0 t) + j\sin(k\omega_0 t)) \right)


Collecting real and imaginary parts leads to:

.. math::
   x(t) = a_0 + \sum_{k=1}^{\infty} \left(
   (A_{-k} + A_k)\cos(k\omega_0 t) + (B_{-k}-B_k)\sin(k\omega_0 t) \right) + \\
   j \sum_{k=1}^{\infty} \left(
     (B_{-k}+B_k)\cos(k\omega_0 t) + (-A_{-k} + A_k)\sin(k\omega_0 t) \right)


The imaginary part of the above equation can only be zero (what we
need for a real signal) in case:

.. math::
   A_{-k} = A_k\\
   B_{-k} = - B_k

Both of these make $a_{-k} = a^\star_k$. 

So we can write:

.. math::
   x(t) = a_0 + 2 \sum_{k=1}^{\infty} \left(
   A_k \cos(k\omega_0 t) -  B_k \sin(k\omega_0 t) \right) 

Remember from high school trigonometry that a summation of a sine and
cosine of the same frequency results in a shifted sine:

.. math::
   x(t) = a_0 + 2 \sum_{k=1}^{\infty} 
   \sqrt{A_k^2+B_k^2} \sin(k\omega_0 t + \phi)

with the phase (shift) $\phi$ given by:

.. math::
   \phi = \arctan\left(\frac{-B_k}{A_k}\right)
   
This last synthesis equation shows that a periodic real valued
function can be written as the sum of sine functions of frequency
$k\omega_0$ for $k>0$.
   
For an **even signal**, i.e. $x(-t)=x(t)$ we have:

.. math::
   a_k \text{is real}

For an **odd signal**, i.e. $x(-t)=-x(t)$ we have:

.. math::
   a_k \text{is imaginary}



Fourier Series Examples
-----------------------

Showing examples requires the calculation of integrals in order to get
at the values of $a_k$ (the analysis equation). Most computer
scientists are not that fond of integral calculations so we try 
to do it here using Sympy, a symbolic math package.

..  ipython:: python

    import matplotlib.pylab as plt
    import sympy as sp

    t, T0 = sp.symbols("t, T0")
    x = sp.Piecewise(
           (1, t % T0<=T0/2),
	   (0, True))
    sp.plot(x.subs(T0,1), (t,-.2, 1.7), xlabel="t", ylabel="x(t)", line_width=3);
    @savefig sp_block.png
    plt.show()

Next we define a function that calculates the $a_k$:

.. ipython:: python

   import numpy as np
	     
   def ak(k):
       return 1/1 * sp.integrate( x.subs(T0,1) * sp.exp(-1j * k * 2 * sp.pi / 1 * t), (t, 0, 1))

   #%timeit print(sp.N(ak(0))) # swwitched off because it takes around 30 s on my machine...
   
But this functions is very very slow... (probably due to the fact that
the piecewise function makes integration hard when Sympy doesn't
recognize the fact that it just has to integrate the exponential
function over half the period. Using this knowledge expicitly we get:

.. ipython:: python

   import numpy as np
	     
   def ak(k):
       return 1/1 * sp.integrate( sp.exp(-1j * k * 2 * sp.pi / 1 * t), (t, 0, 1/2))


   for k in range(-5,6):
	     print(k, sp.N(ak(k))) # almost instantanuous


.. admonition::

   Plot the frequency components and their sum for several values of K (k=0,...,K)



Properties II
-------------


.. admonition:: RVDB
		
   Differentiation


Convergence of the Fourier Series
---------------------------------

.. admonition:: RVDB

   - difference between point wise and integral (??) convergence

   - Gibbs phenomenon