Synthèse spectrale à phases aléatoires

Problématique

En essais vibratoires ou en simulation de charges aléatoires (transport, sismique, turbulence…), on dispose d'un profil de Densité Spectrale de Puissance (DSP) cible — issu d'une mesure sur site ou d'une norme (MIL-STD-810, IEC 60068…). La question est :

Comment générer un signal temporel \(x(t)\) qui soit statistiquement conforme à ce profil DSP, tout en restant aléatoire ?

La réponse est la synthèse spectrale à phases aléatoires, dont la formulation rigoureuse est due à Shinozuka & Deodatis (1991).

Fondement théorique — Méthode de représentation spectrale

Représentation d'un processus stationnaire

Un processus stochastique stationnaire \(x(t)\) de DSP bilatérale \(S_{xx}(f)\) peut s'écrire comme une somme de sinusoïdes à amplitudes déterministes et phases aléatoires uniformes (Rice, 1954 ; Shinozuka, 1972) :

\[x(t) = \sum_{k=1}^{N/2-1} A_k \cos\left(2\pi f_k t + \varphi_k\right)\]

avec :

\(f_k = k \cdot \Delta f\) — fréquences discrètes, \(\Delta f = 1/T\)
\(\varphi_k \sim \mathcal{U}(0, 2\pi)\) — phases indépendantes, uniformément distribuées
\(A_k\) — amplitudes déduites du profil DSP cible

Lien amplitude–DSP

La puissance portée par la raie \(f_k\) dans la bande \([f_k, f_k + \Delta f]\) vaut \(S(f_k) \cdot \Delta f\).

Cette puissance est répartie symétriquement entre \(+f_k\) et \(-f_k\) (spectre bilatéral), d'où :

\[\boxed{A_k = \sqrt{\frac{S(f_k) \cdot \Delta f}{2}}}\]

DSP unilatérale vs bilatérale

En pratique, les profils DSP sont donnés en DSP unilatérale \(G(f)\) définie pour \(f > 0\) : \(\(G(f) = 2 \cdot S(f) \quad (f > 0)\)\) La formule d'amplitude devient alors \(A_k = \sqrt{G(f_k) \cdot \Delta f / 2} = \sqrt{S(f_k) \cdot \Delta f}\). Le code ci-dessous utilise la convention \(S(f)\) bilatérale → facteur \(1/2\) conservé.

Propriété fondamentale

Grâce à l'équidistribution des phases, le signal synthétisé est :

Stationnaire au sens large
Gaussien (par le théorème central limite, somme de nombreux cosinus)
De DSP égale au profil cible en espérance

Implémentation via la TFD inverse (IFFT)

La somme de sinusoïdes s'exprime naturellement comme une transformée de Fourier inverse. En construisant le spectre complexe :

\[\hat{X}_k = A_k \cdot e^{j\varphi_k}\]

et en imposant la symétrie hermitienne \(\hat{X}_{N-k} = \hat{X}_k^*\) (nécessaire pour que \(x(t)\) soit réel), l'IFFT donne directement le signal temporel en \(\mathcal{O}(N \log N)\).

Organigramme

flowchart TD
    A["<b>Profil DSP cible</b>\ndsp_in(:,1) → fréquences\ndsp_in(:,2) → S(f) en g²/Hz"]
    B["<b>Interpolation log-log</b>\nlog S = interp1(log f, log S, log f_pos)\n→ S(fk) sur grille uniforme"]
    C["<b>Amplitudes</b>\nAk = √(S(fk)·Δf / 2)"]
    D["<b>Phases aléatoires</b>\nφk ~ 𝒰(0, 2π)\nrng(graine) pour reproductibilité"]
    E["<b>Spectre complexe unilatéral</b>\nÂk = Ak · exp(j·φk)"]
    F["<b>Spectre bilatéral hermitien</b>\n[0, Â₁…Â_{N/2-1}, 0, Â*_{N/2-1}…Â*₁]"]
    G["<b>IFFT × N</b>\nx(t) = real(ifft(ACC_full)) × N"]
    H["<b>Signal temporel x(t)</b>\nstationnaire, gaussien\nDSP ≈ profil cible"]

    A --> B --> C --> D --> E --> F --> G --> H

    style A fill:#e8f4f8,stroke:#2196F3
    style D fill:#fff3e0,stroke:#FF9800
    style H fill:#e8f5e9,stroke:#4CAF50

Explication ligne par ligne du code MATLAB

fs = 5000;   % fréquence d'échantillonnage (Hz)
T  = 600;    % durée du signal (s)
N  = T*fs;   % nombre d'échantillons = 3 000 000
df = 1/T;    % résolution fréquentielle = 1/600 ≈ 0,0017 Hz
graine = 42;

df = 1/T est la conséquence directe du principe d'incertitude temps-fréquence : plus le signal est long, plus les raies sont fines. Avec \(T = 600\) s et \(F_e = 5000\) Hz, le signal couvre \([0{,}0017 ; 2500]\) Hz avec une résolution de 1,7 mHz.

k     = (1 : N/2-1);
f_pos = k * df;

On construit les fréquences positives uniquement (hors DC \(k=0\) et hors Nyquist \(k=N/2\)), soit les raies \(f_1 = \Delta f, \; f_2 = 2\Delta f, \ldots, f_{N/2-1}\).

log_S = interp1(log10(dsp_in(:,1)), log10(dsp_in(:,2)), log10(f_pos), 'linear', 'extrap');
S = 10.^log_S;

L'interpolation se fait en espace log-log car les profils DSP normatifs (MIL-STD, IEC…) sont des droites en échelle log-log (pentes en dB/octave). Une interpolation linéaire en espace direct fausserait la forme du spectre. 'extrap' prolonge le profil aux extrémités par extrapolation linéaire log-log.

A = sqrt(S * df / 2);

Application directe de la relation amplitude–DSP : \(A_k = \sqrt{S(f_k) \cdot \Delta f / 2}\).

Le facteur \(1/2\) vient de la répartition de la puissance entre fréquences positives et négatives dans le spectre bilatéral.

rng(graine);
phi = 2 * pi * rand(1, N/2-1);

Génération de \(N/2-1\) phases indépendantes \(\varphi_k \sim \mathcal{U}(0, 2\pi)\). La graine fixe (rng(42)) rend la simulation reproductible : deux exécutions du code donnent exactement le même signal — indispensable pour la validation et le débogage.

ACC_pos = A .* exp(1j * phi);

Construction du spectre complexe unilatéral : \(\hat{X}_k = A_k \cdot e^{j\varphi_k}\).

ACC_full = [0, ACC_pos, 0, conj(fliplr(ACC_pos))];

Construction du spectre bilatéral hermitien de longueur \(N\) :

Indice MATLAB	Contenu	Signification
`1`	`0`	DC (\(f=0\)) — nul pour un signal sans composante continue
`2 … N/2`	`ACC_pos`	Raies positives \(f_1 \ldots f_{N/2-1}\)
`N/2+1`	`0`	Nyquist (\(f = F_e/2\)) — nul pour éviter une asymétrie
`N/2+2 … N`	`conj(fliplr(ACC_pos))`	Raies négatives = miroir conjugué

La symétrie hermitienne \(\hat{X}_{N-k} = \hat{X}_k^*\) garantit que real(ifft(...)) est réel sans erreur numérique.

acc = real(ifft(ACC_full)) * N;

MATLAB normalise ifft par \(1/N\), donc on multiplie par \(N\) pour retrouver les amplitudes correctes. Le real(...) élimine les résidus imaginaires d'arrondi numérique (\(\sim 10^{-14}\)).

t = (0 : N-1)' / fs;
acc_reelle = timeseries(acc', t);

Mise en forme pour le bloc From Workspace de Simulink, qui attend un objet timeseries pour injecter le signal dans un modèle de simulation.

Vérification de la conformité spectrale

Après génération, il est indispensable de vérifier que la DSP du signal synthétisé correspond au profil cible. On utilise la méthode de Welch :

=== "MATLAB"

```matlab
[S_check, f_check] = pwelch(acc, hann(8192), 4096, 8192, fs);

figure;
loglog(dsp_in(:,1), dsp_in(:,2), 'r--', 'LineWidth', 2); hold on;
loglog(f_check, S_check, 'b', 'LineWidth', 1);
xlabel('Fréquence (Hz)'); ylabel('DSP (g²/Hz)');
legend('Profil cible', 'Signal synthétisé');
grid on;
```

=== "Python"

```python
from scipy.signal import welch
import numpy as np

f_check, S_check = welch(acc, fs=fs, window='hann', nperseg=8192, noverlap=4096)

import matplotlib.pyplot as plt
plt.loglog(dsp_in[:,0], dsp_in[:,1], 'r--', lw=2, label='Profil cible')
plt.loglog(f_check, S_check, 'b', lw=1, label='Signal synthétisé')
plt.xlabel('Fréquence (Hz)'); plt.ylabel('DSP (g²/Hz)')
plt.legend(); plt.grid(True, which='both')
plt.show()
```

Convergence de l'estimateur de Welch

La DSP estimée convergera d'autant mieux vers le profil cible que le signal est long. Avec \(T = 600\) s et des segments de 8192 points à 5000 Hz, on obtient \(K \approx 730\) segments → variance de l'estimateur très faible.

Origines théoriques

Ce code s'inscrit directement dans la méthode de représentation spectrale (SRM) formalisée par Shinozuka. La filiation est la suivante :

Auteur(s)	Année	Contribution
Rice	1954	Représentation d'un bruit comme somme de sinusoïdes à phases aléatoires
Shinozuka	1972	Application à la simulation de processus stochastiques pour le génie civil
Shinozuka & Jan	1972	Extension aux champs aléatoires multi-variables
Shinozuka & Deodatis	1991	Formalisation complète de la SRM, preuve de convergence vers la gaussianité
Deodatis	1996	Extension aux processus non-stationnaires (SRM évolutive)

Le code présenté est une implémentation directe de Shinozuka & Deodatis (1991), avec deux particularités propres à l'usage industriel :

L'interpolation log-log du profil DSP — adaptée aux spectres normatifs exprimés en dB/octave
L'injection dans Simulink via timeseries — pour simulation dans la boucle avec un modèle de structure

Références

Shinozuka, M. & Deodatis, G. (1991). Simulation of Stochastic Processes by Spectral Representation. Applied Mechanics Reviews, 44(4), 191–204.
Shinozuka, M. (1972). Monte Carlo solution of structural dynamics. Computers & Structures, 2(5–6), 855–874.
Rice, S.O. (1954). Mathematical analysis of random noise. Bell System Technical Journal, 23–24.
Deodatis, G. (1996). Simulation of Ergodic Multivariate Stochastic Processes. Journal of Engineering Mechanics, 122(8), 778–787.
MATLAB Documentation — pwelch, ifft, timeseries
MIL-STD-810H — Environmental Engineering Considerations and Laboratory Tests — Méthode 514.8 (vibrations)