import numpy
                    

def autocorrel(x,N,i,M):
    C = numpy.zeros(N)
    for k in range(i,i+M):
        for n in range(N):
            C[n] += x[k]*x[k-n]
    return C/M        
                    

import math
import random
from matplotlib.pyplot import *

def signal(t):
    return math.sin(2*math.pi*t)+0.5*math.sin(4*math.pi*t)+0.3*random.gauss(0,0.8)
Ns = 4000
T = 8.0
dt = T/Ns
t = numpy.zeros(Ns)
x = numpy.zeros(Ns)
for k in range(Ns):
    t[k] = k*dt
    x[k] = signal(t[k])
figure(figsize=(12,4))
plot(t,x)
xlabel('t')
ylabel('x')
grid()
    
                    

N = 1000
temps = numpy.arange(N)*dt
C = autocorrel(x,N,N,Ns-N)
figure(figsize=(12,4))
plot(temps,C)
xlabel("t")
ylabel("C")
                    

import scipy.io.wavfile as wave
rate,data = wave.read('../../tfd/spectreson3/piano_la3.wav')
n = data.size
duree = 1.0*n/rate
dt = 1.0/rate
Ns = 4000
x1 = data[0:Ns]
t1 = numpy.arange(Ns)*dt
figure()
plot(t1,x1)
xlabel("t (s)")
                    

N=1000
C = autocorrel(x1,N,N,Ns-N)
figure(figsize=(12,4))
plot(numpy.arange(N)*dt,C)
xlabel("t")
ylabel("C")
                     

import numpy.fft
tfd = numpy.fft.fft(x)
freq = numpy.arange(Ns)*1.0/T
S = numpy.square(numpy.absolute(tfd))/Ns
figure(figsize=(8,5))
plot(freq,S)
xlabel("f")
ylabel("S")
axis([0,20,0,S.max()])
                    

c = numpy.fft.ifft(S)
figure(figsize=(12,4))
plot(t,c)
xlabel("t")
ylabel("C")