import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML 
HTML('../style/course.css') #apply general CSS

"""
Explanation: Outline
Glossary
2. Mathematical Groundwork
Previous: 2.4 The Fourier Transform
Next: 2.6 Cross-correlation and auto-correlation




Import standard modules:
End of explanation
"""


import math

from IPython.display import HTML
HTML('../style/code_toggle.html')

"""
Explanation: Import section specific modules:
End of explanation
"""


import math

from matplotlib import rcParams
rcParams['text.usetex'] = True

#def trianglewave(x, T):
#    """
#    This is a sawtooth, though
#    """
#    return np.mod(x/T,1.)*np.logical_and(x>=0,x<=T)

def trianglewave(x, T):
    """
    T is the period.
    """
    return np.abs(2.*(np.mod(x/T,1.)-0.5))-0.5

def boxcar(x,a,b,amp):
    return amp*np.logical_and(x>=a,x<=b)
    
def plottriboxconv(a, b, period):

    # limits of boxcar Play arround with this
#    a = -0.1
#    b = 0.1
    
    # Plotting range
    xrange = [-2., 2.]

    # Create functions
    xpoints = 1000
   
    # Resolution element
    dx = (xrange[1]-xrange[0])/float(xpoints)

    x = np.linspace(xrange[0], xrange[1], xpoints)
    y = boxcar(x, a, b, 1.)

    # boxcar will be normalised to 1. amp = 1./(b-a) works in the limit of many points, but here we do
    # numberofpixelsinbox*dx*amplitude = y.sum *dx*amplitude = 1
    # to take into account numerical effects
    amp = float(xpoints)/((xrange[1]-xrange[0])* y.sum())
    y = boxcar(x, a, b, 1./(b-a))
    ycorr = boxcar(x, a, b, amp)
    z = trianglewave(x, period)

    result = np.convolve(ycorr,z,'same')
    result = dx*result
        
    # Start the plot, create a figure instance and a subplot
    fig = plt.figure()
    ax1  = fig.add_subplot(311)
    fig.tight_layout()
    plt.subplots_adjust(hspace = 0.6)
    
    # Axis ranges
    ax1.axis([xrange[0]+(b-a), xrange[1]-(b-a), z.min()-0.1*(z.max()-z.min()), z.max()+0.1*(z.max()-z.min())])

    # Plot a grid
    ax1.grid(True)

    # Insert lines at x=0 and y=0
    ax1.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax1.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    ax1.plot(x,z,'b-')

    plt.title("Triangle wave", fontsize=14,color='black')
    
    ax2  = fig.add_subplot(312, sharex=ax1) 

    # Axis ranges
    ax2.axis([xrange[0]+(b-a), xrange[1]-(b-a), ycorr.min()-0.1*(ycorr.max()-ycorr.min()), \
              ycorr.max()+0.1*(ycorr.max()-ycorr.min())])

    # Plot a grid
    ax2.grid(True)

    # Insert lines at x=0 and y=0
    ax2.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax2.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    e1 = int(math.ceil(xpoints*(a-xrange[0])/(xrange[1]-xrange[0])))
    ax2.plot(x[:e1],y[:e1],'b-')
    ax2.plot([a, a],[0., amp],'b--')
    e2 = int(math.floor(xpoints*(b-xrange[0])/(xrange[1]-xrange[0])))
    ax2.plot(x[e1:e2],y[e1:e2],'b-')
    e3 = xpoints
    ax2.plot(x[e2:],y[e2:],'b-')
    ax2.plot([b, b],[0., amp],'b--')

    plt.title("Rectangle function", fontsize=14,color='black')
    
    ax3  = fig.add_subplot(313, sharex=ax2) 

    # Axis ranges: mask out border effects
    rmin = result.min()
    rmax = result.max()
    
    # Just to make the result a bit more beautiful if the function is very flat
    if (rmax - rmin) < 0.1:
        rmin=rmin-0.1
        rmax=rmax+0.1

    ax3.axis([xrange[0]+(b-a), xrange[1]-(b-a), rmin-0.1*(rmax-rmin), rmax+0.1*(rmax-rmin)])

    # Plot a grid
    ax3.grid(True)

    # Insert lines at x=0 and y=0
    ax3.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax3.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    plr1 =  int(xpoints*(b-a)/(xrange[1]-xrange[0]))
    plr2 =  int(xpoints*(1-(b-a)/(xrange[1]-xrange[0])))
    
    ax3.plot(x[plr1:plr2],result[plr1:plr2],'b-')

    plt.title("Triangle wave filtered with rectangle function", fontsize=14,color='black')
        
# first two arguments give the position of the rectangle, third the period of the Triangle
plottriboxconv(-0.1, 0.1, 1.0)
# <a id='math:fig:trifilt'></a><!--\label{math:fig:trifilt}-->

"""
Explanation: 2.5 Convolution<a id='math:sec:convolution'></a>
The convolution is an operation connecting two functions, with the result of a mutual broadening. In signal processing, the convolution is often used to represent instrumental broadening of a signal. For any observation, the signal received is "filtered" by an instrumental function. The signal is smeared out. The mathematical description for this effect is the convolution of the function representing the original signal with the instrumental function. In this section, we give a detailed description. Note, the Convolution Theorem, an essential theorem in Fourier analysis and aperture synthesis is discussed in Section 2.7.4 &#10142;.
2.5.1 Definition of Convolution<a id='math:sec:definition_of_convolution'></a>
The convolution $\circ$ is an operation acting on two complex-valued functions.
<a id='math:eq:5_001'></a><!--\label{math:eq:5_001}-->$$
\circ: \left{f\,|\, f:\mathbb{R}\rightarrow \mathbb{C}\right}\,\times\, \left{f\,|\, f:\mathbb{R}\rightarrow \mathbb{C}\right} \rightarrow \left{f\,|\, f:\mathbb{R}\rightarrow \mathbb{C}\right}\
(f\circ g)(x) \,=\, \int_{-\infty}^{+\infty} f(x-t)\,g(t)\,dt
$$
or, in more than one dimension
<a id='math:eq:5_002'></a><!--\label{math:eq:5_002}-->$$
\circ: \left{f\,|\, f:\mathbb{R}^n\rightarrow \mathbb{C}\right}\,\times\, \left{f\,|\, f:\mathbb{R}^n\rightarrow \mathbb{C}\right} \rightarrow \left{f\,|\, f:\mathbb{R}^n\rightarrow \mathbb{C}\right} \, \quad n \in \mathbb{N}\
\begin{align}
(f\circ g)(x_1,\ldots,x_n ) \,&=\, (f\circ g)({\bf x})\
\,&=\, \int_{-\infty}^{+\infty} \ldots \int_{-\infty}^{+\infty} f(x_1-t_1, \ldots , x_n-t_n)\,g(t_1, \ldots, t_n) \,d^nt\
\,&=\, \int_{-\infty}^{+\infty} f({\bf x}-{\bf t})\,g({\bf t}) \,d^nt\end{align}
$$
2.5.2 Properties of Convolution<a id='math:sec:properties_of_convolution'></a>
The following rules apply:
<a id='math:eq:5_003'></a><!--\label{math:eq:5_003}-->$$
\forall\,f,g\in \left{h\,|\, h:\mathbb{R}\rightarrow \mathbb{C}\right}, \quad a \in \mathbb{C}\
\begin{align}
f\circ g \,&=\, g \circ f&\qquad (\text{commutativity})\
(f\circ g)\circ h \,&=\, f \circ (g\circ h)&\qquad (\text{assiociativity})\
f \circ (g + h) \,&=\, (f \circ g) + (f\circ h) &\qquad (\text{distributivity})\
(a\, g)\circ h \,&=\, a \, (g\circ h)&\qquad (\text{assiociativity with scalar multiplication})\
\end{align}
$$
Because (in one dimenstion, to keep it short)
<a id='math:eq:5_002'></a><!--\label{math:eq:5_002}-->$$
\begin{split}
(f\circ g)(x) \,&=\, \int_{-\infty}^{+\infty} f(x-t)\,g(t)\,dt\
&\underset{t^{\prime} = x - t}{=}\, \int_{\infty}^{-\infty} f(t^{\prime})\,g(x-t^{\prime})\,\frac{dt}{dt^{\prime}}dt^{\prime}\
&=\, - \int_{-\infty}^{+\infty} f(t^{\prime})\,g(x-t^{\prime})(-1)\,dt'\
&= (g\circ f)(x)\
((f\circ g)\circ h)(x) \,&=\, \int_{-\infty}^{+\infty} \left[\int_{-\infty}^{+\infty} f((x-t^{\prime})-t)\,g(t)\,dt\right]\,h(t^\prime)\,dt^{\prime}\
&=\, \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(x-t -t^{\prime})\,g(t)\,h(t^\prime)\,dt\,dt^{\prime}\
&=\, \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f((x-t) -t^{\prime})\,h(t^\prime)\,g(t)\,dt^{\prime}\,dt\
&=\, \int_{-\infty}^{+\infty} \left[\int_{-\infty}^{+\infty} f((x-t) -t^{\prime})\,h(t^\prime)\,dt^{\prime}\right]\,g(t)\,dt\
&=\, (f\circ (g\circ h))(x)
\end{split}\qquad \rm 
$$
The last two rules can be easily verified.
2.5.3 Convolution Theorem<a id='math:sec:convolution_theorem_link'></a>
See Section 2.7.4 &#10142;
2.5.4 Convolution Examples<a id='math:sec:convolution_examples'></a>
As said, the convolution is often used to represent an instrumental function. We want to demonstrate this. Let us assume a simple function, the triangle wave and a rectangle function (scaled to an area of 1). If we convolve them we get this:
End of explanation
"""


# first two arguments give the position of the rectangle, third the period of the Triangle
plottriboxconv(-0.5, 0.5, 1.0)
# <a id='math:fig:trifilt'></a><!--\label{math:fig:trifilt}-->

"""
Explanation: Figure 2.5.1: Rectangle-filtered triangle wave
One might assume that one is observing a (co-)sine function. But it can get worse:
End of explanation
"""


from matplotlib import rcParams
rcParams['text.usetex'] = True

def noisycosinewave(x, amplitude, T, sigma):
    """
    T is the period, sigma is the dispersion, amplitude the amplitude
    """
    return amplitude*np.cos(2.*math.pi*x/T)+np.random.normal(scale=sigma, size=x.size)

def boxcar(x,a,b,amp):
    return amp*np.logical_and(x>=a,x<=b)
    
def plotcosboxconv(a, b, period, sigma):

    # limits of boxcar Play arround with this
#    a = -0.1
#    b = 0.1
    
    # Plotting range
    xrange = [-2., 2.]

    # Create functions
    xpoints = 1000
   
    # Resolution element
    dx = (xrange[1]-xrange[0])/float(xpoints)

    x = np.linspace(xrange[0], xrange[1], xpoints)
    y = boxcar(x, a, b, 1.)

    # boxcar will be normalised to 1. amp = 1./(b-a) works in the limit of many points, but here we do
    # numberofpixelsinbox*dx*amplitude = y.sum *dx*amplitude = 1
    # to take into account numerical effects
    amp = float(xpoints)/((xrange[1]-xrange[0])* y.sum())
    y = boxcar(x, a, b, 1./(b-a))
    ycorr = boxcar(x, a, b, amp)
    z = noisycosinewave(x, 1., period, sigma)
    c = np.cos(2.*math.pi*x/period)
    
    result = np.convolve(ycorr,z,'same')
    result = dx*result
        
    # Start the plot, create a figure instance and a subplot
    fig = plt.figure()
    
    ax1  = fig.add_subplot(411)
    fig.tight_layout()
    plt.subplots_adjust(hspace = 0.8)
    
    # Axis ranges
    ax1.axis([xrange[0]+(b-a), xrange[1]-(b-a), c.min()-0.1*(c.max()-c.min()), c.max()+0.1*(c.max()-c.min())])

    # Plot a grid
    ax1.grid(True)

    # Insert lines at x=0 and y=0
    ax1.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax1.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    ax1.plot(x,c,'b-')

    plt.title("Original function (cos)", fontsize=14,color='black')

    ax1  = fig.add_subplot(412)
 
    # Axis ranges
    ax1.axis([xrange[0]+(b-a), xrange[1]-(b-a), z.min()-0.1*(z.max()-z.min()), z.max()+0.1*(z.max()-z.min())])

    # Plot a grid
    ax1.grid(True)

    # Insert lines at x=0 and y=0
    ax1.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax1.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    ax1.plot(x,z,'b-')

    plt.title("Noise added", fontsize=14,color='black')
    
    ax2  = fig.add_subplot(413, sharex=ax1) 

    # Axis ranges
    ax2.axis([xrange[0]+(b-a), xrange[1]-(b-a), ycorr.min()-0.1*(ycorr.max()-ycorr.min()), \
              ycorr.max()+0.1*(ycorr.max()-ycorr.min())])

    # Plot a grid
    ax2.grid(True)

    # Insert lines at x=0 and y=0
    ax2.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax2.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    e1 = int(math.ceil(xpoints*(a-xrange[0])/(xrange[1]-xrange[0])))
    ax2.plot(x[:e1],y[:e1],'b-')
    ax2.plot([a, a],[0., amp],'b--')
    e2 = int(math.floor(xpoints*(b-xrange[0])/(xrange[1]-xrange[0])))
    ax2.plot(x[e1:e2],y[e1:e2],'b-')
    e3 = xpoints
    ax2.plot(x[e2:],y[e2:],'b-')
    ax2.plot([b, b],[0., amp],'b--')

    plt.title("Rectangle function", fontsize=14,color='black')
    
    ax3  = fig.add_subplot(414, sharex=ax2) 

    # Axis ranges: mask out border effects
    rmin = result.min()
    rmax = result.max()
    
    # Just to make the result a bit more beautiful if the function is very flat
    if (rmax - rmin) < 0.1:
        rmin=rmin-0.1
        rmax=rmax+0.1

    ax3.axis([xrange[0]+(b-a), xrange[1]-(b-a), rmin-0.1*(rmax-rmin), rmax+0.1*(rmax-rmin)])

    # Plot a grid
    ax3.grid(True)

    # Insert lines at x=0 and y=0
    ax3.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax3.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    plr1 =  int(xpoints*(b-a)/(xrange[1]-xrange[0]))
    plr2 =  int(xpoints*(1-(b-a)/(xrange[1]-xrange[0])))
    
    ax3.plot(x[plr1:plr2],result[plr1:plr2],'b-')

    plt.title("Noisy function filtered with rectangle function", fontsize=14,color='black')
        
# first two arguments give the position of the rectangle, third the period of the Triangle
plotcosboxconv(-0.1, 0.1, 1.0, 2.5)
# <a id='math:fig:filtnoise'></a><!--\label{math:fig:filtnoise}-->

"""
Explanation: Figure 2.5.2: Triangle wave filtered with rectangle function, another example.
This example illustrates that the process of filtering can destroy information about our signal. However, filtering can also be useful. Given noisy observations of a function, a rectangle function can be used to filter out the noise. This is illustrated in the subsequent example.
End of explanation
"""


from matplotlib import rcParams
rcParams['text.usetex'] = True

def gausshermetian(x, amp, mu, sigma, h3, h4):
    """
    T is the period, sigma is the dispersion, amplitude the amplitude
    """
    y = (x-mu)/sigma
    return amp*np.exp(-0.5*y**2)*(1+h3*(2*np.sqrt(2.)*y**3-3*np.sqrt(2.)*y)/np.sqrt(6.)+h4*(4*y**4-12*y**2+3)/np.sqrt(24))

#amplitude*np.cos(2.*math.pi*x/T)+np.random.normal(scale=sigma, size=x.size)

def boxcar(x,a,b,amp):
    return amp*np.logical_and(x>=a,x<=b)
    
def plotskewedgaussobs(pos1, pos2, boxwidth, sigma, h3, h4):

    # limits of boxcar Play arround with this
#    a = -0.1
#    b = 0.1
    
    # Plotting range
    xrange = [-2., 2.]

    # Create functions
    xpoints = 1000
   
    # Resolution element
    dx = (xrange[1]-xrange[0])/float(xpoints)

    x = np.linspace(xrange[0], xrange[1], xpoints)
    y = boxcar(x, pos1-boxwidth/2., pos1+boxwidth/2, \
               1./boxwidth)+0.5*boxcar(x, pos2-boxwidth/2., pos2+boxwidth/2, 1./boxwidth)

    # boxcar will be normalised to 1. amp = 1./(b-a) works in the limit of many points, but here we do
    # numberofpixelsinbox*dx*amplitude = y.sum *dx*amplitude = 1
    # to take into account numerical effects
    z = gausshermetian(x, 1., 0., sigma, h3, h4)
    
    result = np.convolve(y,z,'same')
    result = dx*result
        
    # Start the plot, create a figure instance and a subplot
    fig = plt.figure()
      
    ax1  = fig.add_subplot(311) 
    fig.tight_layout()
    plt.subplots_adjust(hspace = 0.7)
 
    # Axis ranges
    ax1.axis([xrange[0]+boxwidth, xrange[1]-boxwidth, y.min()-0.1*(y.max()-y.min()), y.max()+0.1*(y.max()-y.min())])

    # Plot a grid
    ax1.grid(True)

    # Insert lines at x=0 and y=0
    ax1.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax1.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    ax1.plot(x,y,'b-')

    plt.title("Original function, impulse", fontsize=14,color='black')
    
    ax2  = fig.add_subplot(312, sharex=ax1)
 
    # Axis ranges
    ax2.axis([xrange[0]+boxwidth, xrange[1]-boxwidth, z.min()-0.1*(z.max()-z.min()), z.max()+0.1*(z.max()-z.min())])

    # Plot a grid
    ax2.grid(True)

    # Insert lines at x=0 and y=0
    ax2.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax2.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    ax2.plot(x,z,'b-')

    plt.title("Instrumental function", fontsize=14,color='black')

    ax3  = fig.add_subplot(313, sharex=ax2) 

    # Axis ranges: mask out border effects
    rmin = result.min()
    rmax = result.max()
    
    ax3.axis([xrange[0]+boxwidth, xrange[1]-boxwidth, rmin-0.1*(rmax-rmin), rmax+0.1*(rmax-rmin)])

    # Plot a grid
    ax3.grid(True)

    # Insert lines at x=0 and y=0
    ax3.axhline(0.,linewidth=1, color = 'k', linestyle='dashed')
    ax3.axvline(0.,linewidth=1, color = 'k', linestyle='dashed')
    
    # Plot function
    plr1 =  int(xpoints*boxwidth/(xrange[1]-xrange[0]))
    plr2 =  int(xpoints*(1-boxwidth/(xrange[1]-xrange[0])))
    
    ax3.plot(x[plr1:plr2],result[plr1:plr2],'b-')

    plt.title("Image: original function filtered with instrumental function", fontsize=14,color='black')
    
# first two arguments give the position of the rectangle, third the period of the Triangle
plotskewedgaussobs(0.0, 1.0, 0.01, 0.1, 0.2, 0.1)
# <a id='math:fig:stars'></a><!--\label{math:fig:stars}-->

"""
Explanation: Figure 2.5.3: Filtering a sine wave completely hidden in the noise with a rectangle function.
Note, that while the signal is not visible in the noisy data, it is partially recovered in the output of the filter.
Representing instrumental functions, it is important to differentiate between the response function in a certain direction and the image of an impulse, which is the reverse of the response function. The function used to represent a measurement via convolution is the image of an impulse function at the origin. This becomes evident, when we convolve two asymmetric functions.
End of explanation
"""


