import math,cmath
import numpy as np
from matplotlib.pyplot import *

def cnCreneau(n):
    if n==0:
        return 0
    else:
        return -1j/(n*math.pi)*(1-(-1)**n)
            

def plotSpectre(cn,nmax):
    indices = np.arange(start=1,stop=nmax,step=2)
    n = indices.size
    spectre = np.zeros(n)
    for k in range(n):
        spectre[k] = 20*math.log10(abs(cn(indices[k]))) 
    stem(indices,spectre)
    xlabel('n')
    ylabel('AdB')
    grid()
figure(figsize=(8,4))
plotSpectre(cnCreneau,20)
            

def serie(cn,nmax,t):
    s = 0.0
    for k in range(1,nmax):
        z = cn(k)*cmath.exp(1j*2*math.pi*k*t)
        s += 2*z.real
    return s
            

def plotSomme(cn,nmax,step):
    t = np.arange(start=0.0,stop=2.0,step=step)
    n = t.size
    u = np.zeros(n)
    for k in range(n):
        u[k] = serie(cn,nmax,t[k])
    plot(t,u)
    xlabel('t')
    ylabel('u')
    grid()
figure(figsize=(8,4))
plotSomme(cnCreneau,20,0.001)
            

figure(figsize=(8,4))
plotSomme(cnCreneau,100,0.001)
            

figure(figsize=(8,4))
plotSomme(cnCreneau,1000,0.0001)
            

import numpy.fft
def calculerSignal(cn,nmax):
    indices = np.arange(start=0,stop=nmax,step=1)
    n = indices.size
    spectre = np.zeros(n,dtype=numpy.complex)
    for k in range(n):
        spectre[k] = cn(indices[k])
    u = numpy.fft.ifft(spectre)*nmax
    return u
u = calculerSignal(cnCreneau,1000)
figure(figsize=(8,4))
plot(u)
grid()
            

def filtrage(H,cn_in):
    def cn_out(n):
        return H(1.0*n)*cn_in(n)
    return cn_out
                

def H(f):
    x = f/10
    s = 1j*x
    return 1.0/(1+s)
                

cnSortie = filtrage(H,cnCreneau)
figure(figsize=(8,4))
plotSomme(cnSortie,100,0.001)
                

def H(f):
    x = f/0.1
    s = 1j*x
    return 1.0/(1+s)
cnSortie = filtrage(H,cnCreneau)
figure(figsize=(8,4))
plotSomme(cnSortie,100,0.001)
                

def butterworth(f,fc,p):
    s = 1j*f/fc
    h = 1.0
    n = 2*p
    for k in range(p):
        c = 2*math.sin(math.pi/n*(k+0.5))
        h = h*(s**2+c*s+1)
    return 1.0/h
                

def H(f):
    return butterworth(f,1,1)
def reponseFrequentielle(H,start,stop):
    freq = np.logspace(start=start,stop=stop,num=1000)
    n = freq.size
    gdb = np.zeros(n)
    phi = np.zeros(n)
    for k in range(n):
        h = H(freq[k])
        gdb[k] = 20*math.log10(abs(h))
        phi[k] = cmath.phase(h)
    return [np.log10(freq),gdb,phi]
[logf,gdb,phi] = reponseFrequentielle(H,-2,2)
figure(figsize=(8,4))
plot(logf,gdb)
xlabel('log(f)')
ylabel('GdB')
grid()
                

figure(figsize=(8,4))
plot(logf,phi)
xlabel('log(f)')
ylabel('phi')
grid()
                

cn_out = filtrage(H,cnCreneau)
figure(figsize=(8,4))
plotSpectre(cn_out,20)
                

figure(figsize=(8,4))
plotSomme(cn_out,100,0.001)
grid()
                

def H(f):
    return butterworth(f,1,2)
cn_out = filtrage(H,cnCreneau)
figure(figsize=(8,4))
plotSpectre(cn_out,20)
                

figure(figsize=(8,4))
plotSomme(cn_out,100,0.001)
grid()