Décomposition Empirique des Modes et Transformation Hilbert-Huang pour l'Analyse de Signaux sous MATLAB

  1. Fondements Théoriques et Processus

Étapes Clés

  1. Décomposition EMD : séparer le signal en composantes intrinsèques (IMF).
  2. Transformée de Hilbert : obtenir le signal analytique pour chaque IMF.
  3. Calcul de la fréquence instantanée : dériver la phase pour l'obtenir.
  4. Construction du spectre temps-fréquence 3D : cartographie fréquence-temps-amplitude.

Définitions Mathématiques

  • Conditions IMF : écart ≤ 1 entre les extrema et les passages par zéro, moyenne locale nulle.
  • Transformée de Hilbert : définie par l'intégrale singulière.
  • Fréquence instantanée : \( \omega(t) = \frac{d}{dt} \theta(t) \), où \( \theta(t) = \angle H[x(t)] \) est la phase analytique.
  1. Implémentation sous MATLAB

Préparation de l'Environnement


% Vérification de la version MATLAB (nécessite R2018a ou ultérieur)
versionActuelle = version('-release');
assert(startsWith(versionActuelle, 'R2018') || startsWith(versionActuelle, 'R2019'), ...
    'Version MATLAB R2018a ou supérieure requise.');

% Ajout du chemin pour la boîte à outils HHT
addpath(genpath('HHT_Toolbox')); % Télécharger depuis MathWorks File Exchange

Code Complet Modifié


%% Préparation des Données
fe = 1000; % Fréquence d'échantillonnage (Hz)
temps = 0:1/fe:1-1/fe; % Vecteur temps (s)
freqA = 50; freqB = 120; % Fréquences du signal
amplitudeA = 0.5; amplitudeB = 0.2;
bruit = 0.15 * randn(size(temps));
signal = amplitudeA * sin(2*pi*freqA*temps) + amplitudeB * sin(2*pi*freqB*temps) + bruit;

%% Décomposition EMD
imfComposantes = emd(signal, 'MaxNumIMF', 6, 'Display', 0); % Limite à 6 IMF

%% Transformée de Hilbert
signalAnalytique = hilbert(signal); % Appliqué au signal original pour illustration

%% Calcul de la Phase et de l'Amplitude Instantanées
phaseInstantanee = angle(signalAnalytique);
phaseDepliee = unwrap(phaseInstantanee); % Dépliage de phase
amplitudeInstantanee = abs(signalAnalytique);

%% Estimation de la Fréquence Instantanée
pasTemps = diff(temps);
variationPhase = diff(phaseDepliee);
freqInstantanee = [phaseDepliee(1); variationPhase ./ pasTemps]; % En rad/s

%% Construction du Spectre 3D
plageFrequence = linspace(0, fe/2, 500); % Plage de fréquences pour le spectre
resolutionTemporelle = 0.01;
tempsEchantillonnes = 0:resolutionTemporelle:temps(end);
matriceSpectre = zeros(length(tempsEchantillonnes), length(plageFrequence));

for idx = 1:size(imfComposantes, 2)
    [tempsInterp, freqInterp] = tfridge(amplitudeInstantanee, freqInstantanee, fe);
    matriceSpectre(:, idx) = interp1(tempsInterp, freqInterp, tempsEchantillonnes, 'linear', 0);
end

%% Visualisation
figure;

subplot(2,2,1);
plot(temps, signal, 'Color', [0.2 0.2 0.2], 'LineWidth', 1.5);
hold on;
plot(temps, imfComposantes(:,1), 'r-', 'LineWidth', 0.8);
plot(temps, imfComposantes(:,2), 'b--', 'LineWidth', 0.8);
legend('Signal Original', 'IMF1', 'IMF2');
title('Signal et Premières Composantes IMF');
xlabel('Temps (s)'); ylabel('Amplitude');

subplot(2,2,2);
imagesc(temps, 1:size(imfComposantes,2), freqInstantanee');
set(gca, 'YDir', 'normal');
title('Distribution des Fréquences Instantanées');
xlabel('Temps (s)'); ylabel('Indice IMF');
colorbar('Label', 'Fréquence (rad/s)');

subplot(2,2,3);
surf(tempsEchantillonnes, 1:size(imfComposantes,2), matriceSpectre', 'EdgeColor', 'none');
title('Spectre Temps-Fréquence 3D');
xlabel('Temps (s)'); ylabel('Indice IMF'); zlabel('Amplitude');
view(45, 30);

subplot(2,2,4);
hs = hht(imfComposantes, fe);
imagesc(hs.time, hs.f, hs.power');
set(gca, 'YDir', 'normal');
title('Spectre de Hilbert');
xlabel('Temps (s)'); ylabel('Fréquence (Hz)');
colorbar('Label', 'Densité Spectrale');

3. Optimisation des Paramètres

Paramètre Valeur Recommandée Impact sur l'Analyse
MaxNumIMF 6-10 Nombre de composantes extraites ; affecte la résolution fréquentielle.
SiftRelativeTol 0.15-0.25 Tolérance relative dans le processus de tamisage ; influence la précision.
MaxSift 80-120 Nombre maximal d'itérations pour la convergence ; contrôle le temps de calcul.

4. Astuces d'Application Technique

Traitement des Effets de Bord


function signalPadded = extensionBord(signal, longueurPad)
    % Méthode de prolongation par réflexion miroir
    debutPad = flipud(signal(1:longueurPad));
    finPad = flipud(signal(end-longueurPad+1:end));
    signalPadded = [debutPad; signal; finPad];
end

Stratégies de Réduction du Bruit


% Débruitage par seuillage d'ondelettes (utilisation de wavelet toolbox)
[c, l] = wavedec(signal, 5, 'sym8'); % Ondelettes Symlets au lieu de Daubechies
seuil = wthrmngr('dw1ddenoLVL', 'penalhi', c);
signalLisse = waverec(wthresh(c, seuil, 's'), l, 'sym8');

5. Interprétation des Résultats

  • IMF1 : composante de haute fréquence (bande 0-40 Hz), souvent associée au bruit.
  • IMF2 : oscillation principale à 50 Hz correspondant au signal utile.
  • Résidu : tendance basse fréquence (environ 0.05 Hz).

Le spectre 3D montre une concentration d'énergie autour de 50 Hz et 120 Hz, confirmant la décomposition correcte.

6. Remarques Pratiques

  1. Entrelacement modal : survient si le rapport de fréquences > 0.5 ; privilégier l'EEMD pour une meilleure séparation.
  2. Effets de bord : utiliser l'extension miroir ou appliquer une fenêtre de Hamming sur les extrémités.
  3. Dérive fréquentielle : pour les signaux longs, segmenter en blocs de 1024 échantillons avec chevauchement.

Étiquettes: MATLAB EMD HHT IMF Analyse fréquentielle

Publié le 22 juillet à 23h57