Polytech Nice Sophia - cours Signal Son et Image pour l'Informaticien

Calculer le spectre d'un signal audio

Support écrit par Jean-Paul Stromboni, pour les élèves de troisième année SI,
nécessite un vidéo projecteur, durée 1h , octobre 2008

Voici ce que vous saurez faire après cette séance :

 Définir le spectre d'un signal audio et le moyen de l'obtenir
 Donner des propriétés de la transformation de Fourier et des
transformées qui seront utilisées dans la suite du cours
 Présenter la transformée de Fourier Discrète et l’algorithme de FFT
 Donner les paramètres de la FFT, et les résultats qu'elle donne
 Calculer et tracer le spectre d'un signal avec MATLAB

Le TD revient sur ces notions avec Matlab

Savez vous déjà répondre aux questions suivantes ?

Quelle est l’information donnée par le spectre Quelle est la résolution fréquentielle d'une
d’un signal audio ? FFT 32 points à fe= 8kHz ?

Que calcule la transformation de Fourier, et à Quel est le spectre de : r (t )  0.5 cos(2 440t )
partir de quelle information ?

Pourquoi une transformée de Fourier discrète? Préciser la fonction l’algorithme de FFT

Le signal x(t) est-il à bande limitée ? Dessiner la fenêtre temporelle

x(t )  cos (2000t )
2
d’expression :  (50t )

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4, page 1

La transformée de Fourier d'un signal continu
donne sa composition fréquentielle ou spectre
 Transformation de Fourier (ou TF) d’un signal



 2 ift
s (t ) 
TF
S( f )  s (t ) e dt

 S(f) est le spectre de s(t)
 i est tel que i 2= - 1
 On se rappellera e 2ift  cos( 2ft )  i sin( 2ft )
 Le spectre est une quantité a priori complexe,
avec un module et un argument. On distingue

 le spectre d'amplitude (qu'on appelle S( f )

spectre par abus de langage):
 le spectre de phase (peu considéré ici) : S ( f )
 Le spectre S(f) et le signal s(t) sont équivalents,
Comparer
on retrouve s(t) à partir de S(f) en utilisant la TF et TF-1
transformation inverse :

 s (t )   S ( f )e 2ift df
1
S ( f ) TF



 Définition : le signal s(t) est à bande limitée si :

f Max , tel que S ( f )  0 pour f  f Max

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 2

Quelques propriétés de la transformation de
Fourier qui pourront être réutilisées ensuite.
 La transformation de Fourier est une
transformation linéaire :
TF [as (t )  bf (t )]  aTF [ s (t )]  bTF [ f (t )]
 Pour a et b réels ou complexes

 La transformée du produit de convo-

lution est le produit des transformées
 Définition du produit de convolution
 
(h * e)(t )   h(t   )e( )d   h( )e(t   )d
 
 Transformée du produit de convolution
s (t )  (h * e)(t ) 
TF
S ( f )  H ( f )  E( f )
 Application aux enveloppes et fenêtres
s(t )  env (t )  x(t ) 
TF
S ( f )  ENV ( f ) * X ( f )

 Il existe une dualité entre transformée

de Fourier et transformée inverse :
 
x(t )   X ( f )e 2 ift
df  x( f )   X (t )e 2itf dt
 


x( f )   X (t )e  2ift dt  TF [ X (t )]


Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 3

Quelques transformées de Fourier
qui pourront être réutilisées :
1. Transformée de la distribution de Dirac
 La distribution ou impulsion de Dirac :


 (t )dt  1  (t )
 (t  0)  0 1
lim t 0  (t )  
0 t
 Transformée de l’impulsion de Dirac :
 


 (t )e  2ift dt  e 0   (t )dt  1


2. Transformée du peigne de Dirac

 Définition de la fonction peigne de Dirac

PeigneT (t )    (t  nT )

PeigneT (t )

T 0 t
T 2T 3T
 La transformée d’un peigne est un peigne
1
TF [ PeigneT (t )]  Peigne 1 ( f )
T T

1 n

T
 ( f 
n  T
)

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 4

 Quelques transformées à réutiliser:
cosinus et fenêtre rectangulaire
3. Transformée de la fonction cosinus :
e 2if t  e 2i f t
0 0

s (t )  cos(2f 0t ) 
2
 La transformée du cosinus contient 2 raies
( ( f  f 0 )   ( f  f 0 ))
S ( f )  TF [ s (t )] 
2

 Compléter :
f

4. Transformée de la fenêtre rectangulaire

 Définition : rectangle de durée T
t 1 (t / T )
s(t )   ( ) 1
T
 T T
t   , 
 2 2 T / 2 0 T /2 t
 La transformée de la fenêtre rectangulaire
est la fonction sinus cardinal ci-dessous :
sin(fT )
S ( f )  T sin c(fT )  T
fT

f
Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 5
Pour calculer le spectre d'un signal discret, on
utilise la Transformée de Fourier Discrète (TFD)
 Définition de la Transformée de Fourier Discrète
Préciser la
X ( f )  TFD[ x ( nTe )]  n x ( nTe )e
  i 2nf / f e
période de

 x(nTe)=xn, n entier relatif, est un signal discret. e 2if / f e

 X(f) est son spectre, c'est une quantité complexe,
périodique sur f, |X(f)| est doté de symétries
 TFD à court terme – fenêtre d’analyse :
 Si on réduit le calcul de la TFD aux N échantillons
x ( nTe )  x n de n  n0 à n  n0  N  1
Quelle est
 on définit une fenêtre d'analyse rectangulaire de la durée de
durée N/fe= NTe : la fenêtre ?

S ( f )  n  n
n  n0  N 1
x n e i 2nf / f e
0

 S(f) est constitué de N lobes de largeurs fe/N et 2fe/N

 Exemple (tracé par Matlab):
fe=8000, N=32, s=0.25*cos(2*pi*500*t)
5 Combien de
spectred'amplitude

lobes :
4
de largeurs
3
période de |
2
S(f)|
1 symétrie ?
la hauteur
0 max ?
0 2000 4000 6000 8000
fréquence(Hz)
Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 6
 L'algorithme de F.F.T. (pour Fast Fourier
Transform) calcule rapidement la TFD.
 Principe de la F.F.T.: calculer seulement N points
de la TFD, les X(fk), f k  kf e / N , k  0 N  1
 si N est une puissance de 2, le calcul est accéléré

fe
 d'où la résolution fréquentielle   f k 1  f k 
N

 Définition de l’algorithme de F.F.T. pour N points

de la TFD calculés et une fenêtre de N échantillons
n0  N 1 2 ink

FFT : X ( fk )  
n  n0
x(nTe )e N

2 ink
1 ( N / 2 ) 1
FFT : x(nTe )  ( )  X ( f k )e
1 N
N k  N / 2
 Application à l’exemple de la page précédente :
4
spectred'amplitude

N=
3
fe =
2

1
retrouver sur
0 l'exemple de
0 2000 4000 6000 8000 la page 6
fréquence(Hz)

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 7

Quelques Transformées de Fourier Discrètes

 la fenêtre rectangulaire est telle que 0n N

 TFD de l'impulsion hn  1, n  0, hn  0, n  0

TFD ( hn )  H ( f )  n 0 hn e 2in ( f / f e )  1
N 1

 TFD d'une constante, c'est également la TFD

de la fenêtre rectangulaire
rn  1, 0  n  N  1

TFD ( rn )  R( f )  n 0 e
N 1 2 in ( f / f )
e

1  e 2iNf / f e
R( f ) 
1  e 2if / f e
e  iNf / f e eiNf / f e  e  iNf / f e
R ( f )  if / f e
e eif / f e  e if / f e
sin(Nf / f e )  i
  i
R ( f )  e i ( N 1) f / f e car sin( ) 
e e
sin(f / f e ) 2i

sin(Nf / f e ) spectre d'amplitude

R( f )  de la fenêtre
sin(f / f e ) rectangulaire

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 8

Annexe : la TFD de la fenêtre rectangulaire
 recopier le résultat obtenu page 8 :
sin(Nf / f e )
 R ( f ) 
rn  1, 0  n  N  1  TFD

sin(f / f e )
 |R(f)| est périodique, de période fe, et s'annule en f  kf e / N

 en effet sin(  k )  sin( ) , k entier  périodicit é 

 et de plus sin(k )  0, k entier

 donc, le numérateur de |R(f)| est périodique, de période f  f e / N

et s'annule en f  kf e / N
 le dénominateur de |R(f)| est périodique, période f  f e

Nf / f e
 De plus, |R(0)|=N : lim f 0 R ( f )  lim f 0 N
f / f e

 D'où l'allure de |R(f)| sur une période entre 0 et fe

R( f ) N N

 f
fe 2 fe
-fe 0 fe
N N

 Que peut-on prévoir si N augmente ?

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 9
 Quelques transformées de Fourier Discrètes
Allure de |R(f)| sur 2 périodes, avec :
N=16, fe= 8000Hz, D=N*Te=0.002s, deltaf=500 Hz
20
spectred'amplitude

15
Que prévoir si N=8
10
prélever la FFT 16 points
5

0
-8000 -6000 -4000 -2000 0 2000 4000 6000 8000
fréquence(Hz)
 Spectre d'un signal sinusoïdal
sn  cos(2nf 0 / f e ), 0  n  N  1
e 2in ( f 0 / f e )  e 2in ( f 0 / f e )
sn  cos(2nf 0 / f e ) 
2

TFD ( sn )  S ( f )  n 0 sn e 2nf / f e
N 1

e 2in ( f  f 0 ) / f e  e 2in ( f  f 0 ) / f e
 n 0
N 1

2
1
 ( R ( f  f 0 )  R ( f  f 0 ))
2
Allure de |S(f)| sur 2 périodes isoler la période entre -
N=16, fe= 8000Hz, D=N*Te=0.002s, deltaf=500 Hz 4kHz et 4kHz
10
spectred'amplitude

5
prélever la FFT 32 points

0
-8000 -6000 -4000 -2000 0 2000 4000 6000 8000
fréquence(Hz)
Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 10
 Le problème de synchronisation de la FFT
On illustre avec le signal s suivant composé
d'une ou de deux fréquences f0 et f1 :
s=a*cos(2*pi*f0*t)+a1*cos(2*pi*f1*t)

1. si f1 n'est pas une fréquence calculée fk, risque

d'erreur sur la fréquence et sur l'amplitude
f0=500, f1=550, a=0, a1=0.25 fe
4 f1  k
spectred'amplitude

N
3

0
0 2000 4000 6000 8000
fréquence(Hz)
2. si deux composantes de fréquence sont trop proches, risque de confusion et
d'erreur

f0=500, f1=625, a=0.25, a1=0.25 fe

5 f 0  f1 
N
spectred'amplitude

0
0 2000 4000 6000 8000
fréquence(Hz)
Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 11
Influence sur la FFT de la taille N de
la fenêtre d'analyse rectangulaire
s= 0.75*cos(2*pi*440*t), D=0.04s, fe=8kHz

1
fenêtrerectangle N=32 fe=
N=
0.5
durée=
signal s

0 fe/N=
-0.5 fe/M=
-1
hauteur
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
temps(s)
spectred'amplitude/N, N=32 M=256
0.4
hauteur du
0.3
lobe
spectre/N

0.2 central :
0.1

0
0 1000 2000 3000 4000 5000 6000 7000 8000
fréquence(Hz)

fenêtrerectangle N=256
1 fe=
0.5
N=
signal s

durée=
0
fe/N=
-0.5
fe/M=
-1
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
temps(s)
spectred'amplitude/N, N=256 M=2048
0.4
hauteur du
0.3 lobe
spectre/N

0.2 central :
0.1

0
0 1000 2000 3000 4000 5000 6000 7000 8000
fréquence(Hz)

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 12

Influence de la fenêtre d'analyse utilisée :
exemple de la fenêtre de Hamming
h( n )  0.54  0.46 cos(2n /( N  1)),0  n  N  1
fenêtredeHamming N=32
1.5
Hamm
Rect
ing
N=
durée=
1
amplitude

0.5

0
0 0.5 1 1.5 2 2.5 3 3.5
temps(s) -3
x10

Comparaison des spectres d'amplitude de ces deux fenêtres :

N=32, M=256 points calculés, D=N*Te=0.004, deltaf=250
1 fe=
Rect N=
spectrenormalisé: fftshift( abs(fft)) / N

0.9 Hamming
durée=
0.8
fe/N=
0.7 maxR=
0.6 maxH=
fmax=
0.5

0.4
nblobesH
0.3 largeur
0.2
nblobesR
0.1
largeurs:
0
-4000 -3000 -2000 -1000 0 1000 2000 3000 4000
fréquence(Hz)

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 13

 Illustrer l'influence de la taille N de la
fenêtre d'analyse de Hamming sur la FFT
fenêtredeHamming N=32
1

0.5
signal s

-0.5

-1
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
temps(s)
spectred'amplitude/(N*0.53), N=32M=256
0.4
spectre/(N*0.53)

0.3

0.2

0.1

0
0 1000 2000 3000 4000 5000 6000 7000 8000
fréquence(Hz)

fenêtredeHamming N=256 détecter la

1
fenêtre de
0.5 Hamming
signal s

-0.5

-1
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
temps(s)
spectred'amplitude/(N*0.53), N=256 M=2048
0.4
spectre/(N*0.53)

0.3

0.2

0.1

0
0 1000 2000 3000 4000 5000 6000 7000 8000
fréquence(Hz)

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 14

Paramètres de la fonction Matlab fft(.)
 spectre= fft(s(n0:n0+N-1)) :
utilise N échantillons successifs n0,n0+1, …
pour calculer N points du spectre dans le
vecteur spectre :
1. N, taille de la fenêtre d'analyse :
 durée du signal analysé : D= (N-1)Te
 résolution fréquentielle :  = fe/N
2. fk=[0,fe/N,2*fe/N, … (N-1)*fe/N] contient les
fréquences des valeurs du spectre calculées
3. abs(spectre)/N donne l'amplitude exacte des
raies du spectre (sauf en cas de problème de
synchronisation, rappelez vous !!)
 spectre=fft(s(1:N),M) : calcule M points de
la TFD au lieu de N, avec M>N:
 fk= 0, fe/M, 2*fe/M, … , (M-1)*fe/M
 le spectre d'amplitude reste abs(spectre)/N
 la résolution fréquentielle reste  = fe/N
 spectre= fft(s(1:N).*window(1:N)) :
applique une fenêtre d'analyse window (comme
une enveloppe temporelle) au signal avant de
calculer la FFT, avec l'effet suivant :
 on retrouve le spectre de window, dupliqué sur
chaque composante fréquentielle du signal
analysé (cf. ci-après)
 exemple : window égale hamming(N)pour une
fenêtre de Hamming

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 15

 Calculer et tracer la F.F.T. avec Matlab
1024points dans la fenêtre
1
fe=8192;
t=[0:8191]/fe;
s=0.75*cos(2*pi*1760*t); 0.5
n0= 512; N=1024;
fenetre =zeros(size(s));

signal
fenetre(n0:n0+N-1)=1; 0
% tracer la fenêtre d'analyse
plot(t, s.*fenetre), grid
-0.5
xlabel('temps(s)')
ylabel('signal')
-1
0 0.2 0.4 0.6 0.8 1
temps(s)
Spectre d'amplitudeentre 0et fe
0.4

% spectre d'amplitude sur [0,fe[

0.3
spectred'amplitude(FFT/N)

specampl=abs(fft(s(n0:n0+N-1),N));
f=(0:N-1)*fe/N;
plot(f,specampl/N), 0.2

grid
xlabel('fréquence (Hz)') 0.1
ylabel('spectre d''amplitude')
0
0 1000 2000 3000 4000 5000 6000 7000 8000 9000
fréquence(Hz)

Spectre d'amplitudeentre -fe/2et fe/2

0.4
%% spectre dans [-fe/2,fe/2[
f=[-N/2:N/2-1]/N*fe;
0.3
plot(f,abs(fftshift(specampl))/N)
spectred'amplitude(FFT/N)

grid,
xlabel('fréquence (Hz)') 0.2
ylabel('spectre d''amplitude')
Rappeler le spectre du signal:
0.1

0
-5000 0 5000
fréquence(Hz)

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 16

On obtient le spectrogramme en glissant
la fenêtre d'analyse sur la durée du signal
fe=8000; N=4096; % script spectrogramme.m
t=[0:16000]/fe;
s=0.5*cos(2*pi*1000*t)+ 0.75*cos(4000*pi*t);
f=[-N/2:N/2-1]*fe/N; Durée de la
spec= fftshift(fft(s(1:N))) fenêtre ?
plot(f,abs(spec)), grid, figure
spectrogram(s,N,N/2,N,fe,'yaxis')
colorbar fe ?
400

350

300

250
Retrouver
200
amplitudes et
150
fréquences
100

50

0
-5000 0 5000

Noter pour le spectrogramme :

1. la fenêtre d’analyse ci-dessous est une fenêtre de Hamming
2. le code de couleurs est affiché à droite (en dB)

4000
-20
3500

-40
3000
Frequency (Hz)

-60
2500

-80
2000

1500 -100

1000 -120

500 -140

0
0.5 1 1.5
Time

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 17

Annexe : Des propriétés utiles de (t)

 rappeler la définition de la
distribution ou impulsion  
 (t )dt  1  (t )
de Dirac, son symbole et  (t  0)  0
avec ses propriétés : lim t 0  (t )  
0 t
 Rappeler la définition du 

produit de convolution (h * e)(t )   h(t   )e( )d


vue page 3 : 
  h( )e(t   )d


 Représenter ci contre
 t (t )   (t  t0 )
0 1
t

t0
 Comment expliquer cette 
propriété ? 
s(t ) (t  t0 )dt  s(t0 )

s(t ) (t  t0 )  s(t0 ) (t  t0 )
c'est le modèle mathématique de l'échantillonnage de s(t) en t=t0
 Et que vaut le produit de
convolution ci-contre ?
( t0 * e)(t )

( t0 * e)(t )    (t  t0   )e( )d  e(t  t0 )


le produit de convolution de e(t) par (t-T) décale le signal e de T

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 18
 Autre représentation du spectrogramme
 pour le signal vocal, on sait que la durée de la
fenêtre d’analyse ne doit pas dépasser 30ms (?)
 si fe=8 kHz, c’est une fenêtre de 240 échantillons.
 On calcule la FFT de la fenêtre,
 on déplace la fenêtre et on recommence
 On regroupe les résultats dans un spectrogramme,
en 3D (cf. ci-dessous) ou en 2D (cf. Goldwave)
 Quelle est ici la résolution fréquentielle ?
 Comment obtenir une fenêtre de 20ms, sachant que
fe=22050Hz ? Donner la résolution fréquentielle.
 Voici le spectrogramme de piano_c3.wav tracé par
WaveLab : retrouver les informations de fréquence
fondamentale, durée du signal, enveloppe …

Calculer le spectre d'un signal audio, cours S.S.I.I. n°4 , page 19