Dans ce chapitre, vous allez apprendre à tracer le spectre de données expérimentales.
Les données expérimentales correspondent ici au parcours dans un plan d'une petite bille chargée dans un champ magnétique.
La position (x et y) de cette bille a été relevée toutes les millisecondes pendant 4 secondes.
Vous devriez maintenant être habitué à la procédure : importez les données expérimentales Trajectoire.dat.
%Importation des données, il y a trois colonnes
data = load('Trajectoire.dat')
t = data(:,1);% la première colonne est une colonne qui correspond aux instants d'échantillonnage
x = data(:,2); %La deuxième est une position horizontale
y
= data(:,3); %la troisième est une position verticale
Afin de se faire une petite idée de ce que vous avez comme donnée, visualisons tout cela :
%représentation des données
figure(1)
plot(x(t<1/10),y(t<1/10)) %La condition t<10 limite la quantité de points tracés
%partie cosmétique
title('Trajectoire en 1/10e de s')
xlabel('x')
ylabel('y')
% On utilise une figure
différente pour tracer la trajectoire car les axes sont différents
figure(2)
plot(t,x,'k')
hold on
plot(t,y,'r')
hold off
legend({'x(t)' 'y(t)'})
Vous devriez obtenir :
Vous obtenez quelque chose de semblable à une particule qui suit une trajectoire aléatoire. Mais si vous regardez de plus près, juste ce qui se passe sur x et sur y :
Vous pouvez donc en conclure que la particule oscille à la fois horizontalement (
est périodique) et verticalement ( aussi). Il semble par ailleurs qu'il y ait des oscillations rapides et des oscillations lentes, ce que vous pouvez voir en zoomant un peu :
Détail des variations de et en fonction du temps. Pour , une variation de période 5 ms (20 Hz) se superpose à une autre 1ms (100 Hz). Pour , c'est similaire, sauf que les fréquences sont différentes. Pouvez vous les deviner?
Donc, vous pouvez vous attendre à ce que la représentation fréquentielle traduise cela. C'est-à-dire que le signal ait des composantes fortes sur les vecteurs de Fourier dont la fréquence est la même que celle que nous devinons sur le signal.
Il vous faut donc construire les vecteurs de Fourier, afin de pouvoir projeter les données dessus. Ce qui nous permettra ainsi d'obtenir le spectre.
%nous construisons un premier vecteur de Fourier.
N = length(x);% N est la longueur du vecteur de données expérimentales
k = [0:N-1]';% k est un indice variant de 0 à N-1, ici nous construisons un vecteur colonne
F0 = exp(1i*2*pi*0/N*k);
Une fois le vecteur construit, vous pouvez vérifier qu'il s'agit bien d'une suite de 1, donc un vecteur constant (car l'exponentielle complexe de fréquence 0 vaut toujours 1).
Puis, vous pouvez projeter les données expérimentales sur ces vecteurs.
%Premier terme de la Transformée de Fourier/Projection de x sur le premier vecteur de Fourier
tfx(1) = sum(x .* conj(F0))% On utilise conj pour conjugé.
%La même chose pour y:
tfy(1) = sum(y .* conj(F0))% On utilise conj pour conjugé.
La console répond :
tfx = -23.168
tfy = 8006.1
Notez que si vous divisez ces valeurs par N, vous obtiendrez la valeur moyenne (vous pouvez vérifier que y est bien de valeur moyenne - ou de composante constante - 2).
>> tfy/N
ans = 2.0010
Bon, eh bien il ne vous reste plus qu'à calculer toutes les composantes fréquentielles possibles.
Pour cela, il faut projeter les données sur tous les vecteurs de Fourier :
clear tfx tfy freq% permet le nettoyage des tableaux tfx, tfy et freq
for n = 0:N-1
Fn = exp(1i*2*pi*n/N*k); %Création du n-ieme vecteur de Fourier qu'on appelle Fn
p = n+1; %position
dans le tableau, les indices matlab commencent à 1 et pas à 0
freq(p) = n/N;%fréquence de Fourier associée
tfx(p) = dot(x,Fn); % La même chose que sum(x .* conj(Fn)). matlab a une fonction
toute prète qui fait ça. ("dot product" en anlgais signigie produit scalaire)
tfy(p) = dot(y,Fn);
end
tfx(1:3)'%affiche les trois premiers termes le ' est pour la lisibilité
La console devrait vous donner les trois premiers termes :
ans =
-23.1676 - 0.0000i
3.6870 + 17.7100i
6.0581 - 8.2493i
Nous retrouvons bien le -23 trouvé pour le premier terme de la TF de x, mais qu'en est-il des suivants ?
Ce sont des nombres complexes !! Un peu une surprise, mais pas tant que cela, puisque nous utilisons des vecteurs complexes, les composantes peuvent l'être également.
Pour l'instant, nous nous préoccupons juste du module de ces composantes, module qui s'obtient avec abs :
abs(tfx(1:3)')
ans =
23.168
18.090
10.235
Pour finir, vous pouvez visualiser ces composantes. Comme il s'agit de composantes complexes, nous allons simplement en visualiser le module (le spectre).
%Tracé du spectre
figure(3)
Te = 1e-3;
stem(freq/Te,abs(tfx)/N);%stem est un type de plot souvent utilisé pour représenter les spectres
xlabel('Fréquence en Hz')
ylabel('Composantes de Fourier (normalisées)')
Vous devriez avoir :
Il se passe clairement des choses bizarres sur les très hautes fréquences (cela est dû au fait que le signal soit échantillonné et réel). Ne vous en préoccupez pas pour le moment et zoomez sur la zone 0-120 Hz :
xlim([0 120]) %limite la zone de traçage du graph
Vous remarquez donc que le signal contient une forte composante à la fréquence 18 Hz, mais également une composante vers 83 Hz. Vous pouvez vérifier sur le signal lui-même qu'il s'agit bien des fréquences que nous avions observées sur le signal lui-même.
Si vous faites de même avec le spectre de y, vous obtiendrez :
Que pouvons-nous dire ?
Autant les composantes à 20 et 115 Hz étaient à peu près visibles dans le signal de départ (et encore, ce n'était pas évident), autant celle à 60 était bien invisible.
Voilà pourquoi la représentation fréquentielle est parfois déterminante, elle permet de voir des choses pas forcément visibles dans la représentation temporelle.
Vous pouvez vérifier que vous obtenez exactement les mêmes résultats en utilisant la fonction toute faite de Matlab qui calcule les transformées de Fourier : fft.
tfx_matlab = fft(x);% calcul de la TF via la fonction built-in de matlab
plot(freq/Te,abs(tfx))
hold on
plot(freq/Te,abs(tfx_matlab) + 100,'r') % le +100 est là pour décaler les courbes qui sont sinon parfaitement superposées
hold
off
fft de Matlab. fft est alors : freq = [0:N-1]/(N*Te).abs).