Autocorrélation d'un signal

2.a. Définition

Soit x(t) un signal. La fonction d'autocorrélation temporelle ([1]) est définie par :

$$C\left(\tau \right)=\underset{T\to \infty }{lim}\frac{1}{T}{\int }_{-T/2}^{T/2}x\left(t\right)x(t-\tau )dt\text{(1)}$$

Il s'agit donc de la moyenne temporelle du produit du signal par lui-même décalé d'un temps τ. La fonction d'autocorrélation est paire; on peut donc l'étudier pour τ>0.

Les signaux réels sont limités dans le temps. L'intégrale définissant la fonction d'autocorrélation est alors calculée pour une durée T finie assez grande.

Soit x_k un signal numérique obtenu par échantillonnage avec une période T_e. Soit M le nombre de points utilisés pour calculer la moyenne (T=MT_e). La fonction d'autocorrélation discrète est définie par :

$$C_n=\frac{1}{M}\sum _{k=i}^{i+M-1}{x}_k{x}_{k-n}\text{(2)}$$

La moyenne est ainsi calculée pour les points d'indices i à i+M-1. Soit N le nombre de points de la fonction d'autocorrélation. L'indice n varie de 0 à N-1.

2.b. Fonction de calcul

import numpy
                    

La fonction suivante effectue le calcul de l'autocorrélation pour un signal donné par un tableau numpy x. L'indice de départ i doit être supérieur à N.

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        
                    

2.c. Exemple

On considère comme exemple un signal aléatoire obtenu en ajoutant un bruit gaussien à un signal périodique :

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()
    
                    

figA.pdf

Voyons son autocorrélation calculée sur une durée 1000T_e :

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")
                    

figB.pdf

Comme on le voit sur cet exemple, la fonction d'autocorrélation permet d'obtenir la période d'un signal, même lorsqu'il est fortement perturbé par le bruit.

Par ailleurs, la valeur de l'autocorrélation pour τ=0 est :

print(C[0])
--> 0.6816976624914507

Pour un signal de valeur moyenne nulle, il s'agit de la fluctuation quadratique moyenne :

$$C\left(0\right)=\underset{T\to \infty }{lim}\frac{1}{T}{\int }_{-T/2}^{T/2}x{\left(t\right)}^{2}dt\text{(3)}$$

La fonction d'autocorrélation est très utile pour analyser un signal sonore :

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)")
                    

figC.pdf

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

figD.pdf

L'instant du premier maximum donne la période du signal.

3.a. Énoncé

La transformée de Fourier du signal x(t) est :

$$X\left(f\right)={\int }_{-\infty }^{\infty }x\left(t\right)exp(-i2\pi ft)dt\text{(4)}$$

Le signal s'exprime en fonction de sa transformée de Fourier par la transformée de Fourier inverse :

$$x\left(t\right)={\int }_{-\infty }^{\infty }X\left(f\right)exp\left(i2\pi ft\right)df\text{(5)}$$

La densité spectrale de puissance du signal est définie par :

$$S\left(f\right)=\underset{T\to \infty }{lim}\frac{1}{T}|X\left(f\right){|}^2\text{(6)}$$

Le théorème de Wiener-Khinchine ([2]) énonce que la densité spectrale de puissance est la transformée de Fourier de la fonction d'autocorrélation :

$$S\left(f\right)={\int }_{-\infty }^{\infty }C\left(\tau \right)exp(-i2\pi f\tau )d\tau \text{(7)}$$

Inversement, la fonction d'autocorrélation est la transformée de Fourier inverse de la densité spectrale :

$$C\left(\tau \right)={\int }_{-\infty }^{\infty }S\left(f\right)exp\left(i2\pi f\tau \right)df\text{(8)}$$

3.b. Exemple

On calcule la transformée de Fourier discrète du signal échantillonné précédent, puis son spectre de puissance.

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()])
                    

figE.pdf

La transformée de Fourier inverse permet d'obtenir la fonction d'autocorrélation :

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

figF.pdf