Intégration de spectre et écart-type

**FelixFelix** · 23/03/2020, 15h29

Bonjour,

Je suis doctorant en aérodynamique et je cherche à exploiter des données obtenus en soufflerie. J'essaie d'étudier le taux de turbulence que l'on trouve étonnant. Je possède un fichier avec la tension d'une sonde (250 000 échantillons à 25 kHz, avec un filtre passe-bas pour éviter l'aliasing à 11.3 kHz), cette tension peut être convertie en vitesse avec une simple loi, j'ai donc de manière équivalente 250 000 échantillons de vitesse.

Je cherche à étudier l'influence des très basses fréquences sur la turbulence que l'on définit ainsi: Tu = Ecart-Type(U) / U_moyen, U étant la vitesse. Pour faire ça j'essaie de faire la chose suivante:
-Effectuer la PSD de la partie fluctuante de la vitesse (pour les 250 000 échantillons, on appelle cette composante u[i] = U[i] - U_moyen)
-Intégrer cette PSD sur toute la gamme de fréquence sauf la plus basse (environ 6Hz seront coupés), pour retomber sur l'écart-type(U) sans les basses fréquences (sauf erreur, intégrer la PSD de u[i] = U[i] - U_moyen doit me faire retomber sur la valeur de l'écart-type)

Dans une premier temps, j'essaie d'intégrer l'écart-type en prenant toute la gamme de fréquence pour vérifier qu'on obtient bien la même valeur que par la commande np.std(U), mais ceci ne marche pas. Sauriez-vous d'où vient l'erreur? Je n'arrive pas à savoir si elle vient d'un souci de raisonnement de ma part ou simplement d'implémentation.

Je mets ci-dessous ma boucle principale simplifiée (j'ai 36 points avec 250 000 échantillons de vitesse)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
 
for i in range (nbdepoints):    
 
    #Initialisation des variables temporaires
    tension_FC=[]
    vitesse_FC=[]
 
    #Récupération de la position et des valeurs de tension
    position=np.genfromtxt(filchaud,dtype=str,usecols=1, unpack=True, skip_header=3, max_rows=1)
    tension_FC =np.genfromtxt(filchaud,dtype=str,usecols=(1), unpack=True, skip_header=2,max_rows=nb_echantillons)
    position=str(position)
 
 
    #Conversion virgules en point et tension en vitesses
    position=convertisseur1(position)
    tension_FC=convertisseur(tension_FC)
    vitesse_FC=convertisseurvitesse(tension_FC)
    for j in range(nb_echantillons):
        U_inst[i][j]=vitesse_FC[j]
 
    #alimentation des données à garder
    U[i]=np.mean(vitesse_FC)
    u_rms[i]=np.std(vitesse_FC)
    position_longi[i]=position+position_depart
    Tu[i]=u_rms[i]/U[i]
 
 
    if i==36:                                                                                                    ######## c'est une version simplifiée où je n'étude que le dernier point pour aller plus vite
        f1,PSD1=signal.welch(vitesse_FC,frequence_acqui, nperseg=4096)        
 
        #On calcule u' adimensionné par U                                                           ######## vitesse fluctuante appelée u ci-dessus
        u_prime_ad=[]
        for j in range(len(vitesse_FC)):
            u_prime_ad.append((vitesse_FC[j] - U[i]))
 
        f2,PSD2=signal.welch(u_prime_ad,frequence_acqui, nperseg=4096) 
 
 
#Intégraton par la méthode des trapèzes
a=np.sqrt(np.trapz(PSD2,f2))  #Sur toutes les fréquences                                    ##### Logiquement je pensais obtenir ici la même valeur que par np.stp(Vitesse)
a1=np.sqrt(np.trapz(PSD2[1:],f2[1:])) #On enlève la première bande de fréquences
a2=np.sqrt(np.trapz(PSD2[2:],f2[2:])) #On enlève les deux premières

Merci d'avance

**lg_53** · 23/03/2020, 16h08

Bonjour Felix

Tu es sur un forum Python. Alors la PSD d'une fréquence ... Ca doit être super intéressant, mais je n'ai juste absolument aucune idée de ce que ça signifie.

Par contre, côté programmation, ce que je comprend, c'est que tu essaies de calculer une standard déviation (un peu modifiée), et donc a juste titre tu recodes cette algo, et tu vérifies que si tu ne fais pas ta petite modif alors tu retombes bien sur la std classique, très bien.

Mais dans ce cas, on n'a pas votre fichier d'input. Donc prenez un petit vecteur, que vous vous donnez en entrée et montrer nous comment vous y appliquer votre fonction. Car là entre le fait qu'on ait pas votre fichier et que votre algo de calcul est mélangé à la lecture du fichier, pas facile de s'y retrouver et de savoir quel doit être le résultat qu'on est censé obtenir.

**FelixFelix** · 23/03/2020, 16h28

Bonjour,

Déjà merci pour la réponse. Je me doute bien que tout le monde ne connaît pas les PSD (c'est une courbe donnant le niveau d'énergie pour chaque fréquence d'un signal en gros), même si j'ai pris soin d'aller la partie calcul scientifique! L colonne de droite ci-dessous représente la tension
3,553608E+0 2,347472E+0
3,553284E+0 2,347629E+0
3,553608E+0 2,347413E+0
3,553284E+0 2,347713E+0
3,553284E+0 2,347351E+0
3,553608E+0 2,347543E+0
3,553931E+0 2,347591E+0
3,553608E+0 2,347390E+0
3,553931E+0 2,347705E+0
3,553284E+0 2,347383E+0
3,553608E+0 2,347727E+0
3,553608E+0 2,347356E+0
3,553608E+0 2,347872E+0
3,553931E+0 2,347253E+0
3,553284E+0 2,347756E+0
3,553931E+0 2,347361E+0
3,553608E+0 2,347770E+0
3,553608E+0 2,347462E+0
3,553931E+0 2,347587E+0
3,553284E+0 2,347525E+0
3,553608E+0 2,347482E+0
3,553608E+0 2,347695E+0
3,553284E+0 2,347233E+0
3,554254E+0 2,347759E+0
3,553284E+0 2,347286E+0
3,553608E+0 2,347729E+0
3,552638E+0 2,347129E+0

Pour chaque point, je convertis cette tension en vitesse, puis j'en déduis les valeurs moyennes et écart-types (que j'appelle rms même si ce n'est pas hyper pertinent comme nom)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
vitesse_FC=convertisseurvitesse(tension_FC)
U[i]=np.mean(vitesse_FC)
u_rms[i]=np.std(vitesse_FC)

Ensuite, je calcule ma partie fluctuante (l'équivalent d'un signal AC en électricité et je fais sa PSD)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
#On calcule u' adimensionné par U
 u_prime_ad=[]
 for j in range(len(vitesse_FC)):
       u_pime_ad.append((vitesse_FC[j] - U[i]))
 
f2,PSD2=signal.welch(u_prime_ad,frequence_acqui, nperseg=4096)

Puis j'intègre cette PSD par la méthode des trapèzes

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

a=np.sqrt(np.trapz(PSD2,f2))

Le calcul direct de l'écart-type me donne environ 0.26 m/s, le calcul avec la méthode PSD donne 0.06 m/s

**lg_53** · 24/03/2020, 22h18

Bonsoir

Nous n'avons pas la fonction convertisseur_vitesse.
Donc soit vous nous donnez cette fonction, soit encore mieux, vous nous donné directement le vecteur vitesse_FC. Car après tout votre problème démarre là (indépendamment de la manière dont vous construisez vitesse_FC).

Le problème est pourquoi

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

ma_fonction_std(vitesse_FC)

est différent de

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

np.std(vitesse_FC)

n'est ce pas ?

Où la fonction ma_fonction_std est définie comme suit :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
def ma_fonction_std(vitesse_FC): 
    Ui=np.mean(vitesse_FC)
 
    u_prime_ad=[]
    for j in range(len(vitesse_FC)):
         u_pime_ad.append((vitesse_FC[j] - Ui))
 
    f2,PSD2=signal.welch(u_prime_ad,frequence_acqui, nperseg=4096)
 
    return np.sqrt(np.trapz(PSD2,f2))

on est bien d'accord ?
avec un import signal au début du code je présume ?

Il faut nous donner ce qu'on appelle dans notre jargon un MWE (Minimum Working Exemple).

**FelixFelix** · 25/03/2020, 10h25

Bonjour,

En fait comme le code est un peu long je pensais plutôt vous faciliter la tache en allégeant au maximum la partie que je copiais sur le site, au temps pour moi!
Hier j'ai résolu, ou plutôt contourné le problème. Comme je voulais tester l'influence des basses fréquences, j'ai trouvé un morceau de code permettant d'appliquer un filtre passe-haut sur mon signal et ça a bien fonctionné. Ca n'est pas exactement ce que je voulais faire au début mais ça me permet de répondre à ma question physique, donc au bout du compte j'ai ma réponse, même si ça me frustre un peu de ne pas savoir d'où vient le problème, je ne vais pas vous prendre plus de temps. Merci pour vos réponses en tout cas.

Félix

**Julien N** · 25/03/2020, 17h13

Salut !

N'abandonnons pas. Le principe est juste, l'écart type d'une série est égal à l'écart type du spectre de cette série (Parceval?). Et l'écart type d'un spectre c'est bien la racine de l'integral du spectre.

Si le résultat est différent cela peut être dû à une valeur moyenne qui traîne (mais j'ai cru voir que vous l'aviez enlevée), ou à un problème d'unité. Êtes vous certain de l'unité du spectre ? Je le suis fait avoir quelques fois...

J

**FelixFelix** · 25/03/2020, 17h38

Hello!

A priori oui, si je prends ma vitesse fluctuante, je suis en m/s, mon écart-type est donc en m/s aussi. Pour mon spectre, si je ne me gourre pas, là c'est des m²/s²/Hz. En intégrant, ça donne des m²/s², donc m/s à la racine. C'est sûr que le problème vient des très basses fréquences mais je ne sais pas d'où.

**Julien N** · 26/03/2020, 08h44

Salut,

Quand j'ai besoin de faire ce genre d'analyses, je me repose sur cette fonction que je conserve précieusement :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
import numpy as np
 
def fourier(t, signal):
    """Fast Fourier Transform of a time serie"""
    N, dt = len(t), t[1]-t[0]
    signal = signal - np.mean(signal)
    dft = np.fft.fft(signal)
    freq = np.fft.fftfreq(N, d=dt)[0:int(N/2)]
    ampl = 2 * np.abs(dft[0:int(N/2)]**2) / (N / dt)
    return freq, ampl

Prenons un exemple simple :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Create dummy signal of period 5s and amplitude 1 whatever
t = np.linspace(0, 100, 1000)
s = np.sin((2 * np.pi / 5) * t)
 
# Compute standard deviation of the signal s
s_std = np.std(s)
 
# Determine the PSD of the signal
freq, ampl = fourier(t, s)
 
# Compute spretrum's zero-order momentum and deduce the standard deviation
m0 = np.trapz(ampl, freq)
psd_std = np.sqrt(m0)
 
print(s_std, psd_std)

Ce qui chez moi donne la même chose, à savoir 0.7067.

J

**FelixFelix** · 26/03/2020, 11h27

Intéressante ton idée de faire la FFT. Je l'ai faite en utilisant ton code, et effectivement je retrouve la valeur RMS en intégrant. J'ai remarqué que la FFT donnée par ton code me donnait une résolution fréquentielle beaucoup plus grande (bandes de fréquences de 0.1 Hz). Du coup j'ai modifié un peu les paramètres de ma PSD pour m'approcher de ça, et effectivement ça modifie pas mal de le résultat (Pour le jargon, U est le signal total, u_prime le signal à valeur moyenne nulle):

u_rms = 0.26344079913034624
u_rms_uprime = 0.2634407991303463

***PSD nperseg = 4096 (résolution fréquentielle : 6 Hz)
Integration uprime PSD Totale = 0.06072585635955464
Integration U PSD Totale = 0.060725856359554624

***PSD nperseg = 194304 (résolution fréquentielle: environ 0.1 Hz)
Integration uprime PSD Totale = 0.2134084120231004
Integration U PSD Totale = 0.213408412023101

***FFT
Integration uprime FFT Totale = 0.2634407981687146
Integration U FFT Totale = 0.2634407981687146

La résolution spectrale semble avoir une grosse influence. Ça n'explique pas tout mais en tout cas ça m'apporte pas mal d'informations et ça me donne des pistes à creuser. Merci beaucoup pour cette idée!

PS: Est-ce que tu voudrais bien m'expliquer ton petit algorithme stp? En particulier les deux lignes indiquées ci-dessous, je suis un peu perdu notamment à cause des [0:int(N/2)].

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
def fourier(t, signal):
    """Fast Fourier Transform of a time serie"""
    N, dt = len(t), t[1]-t[0]
    signal = signal - np.mean(signal)
    dft = np.fft.fft(signal)
    freq = np.fft.fftfreq(N, d=dt)[0:int(N/2)]      => Ca
    ampl = 2 * np.abs(dft[0:int(N/2)]**2) / (N / dt) => Et ça
    return freq, ampl

**Julien N** · 26/03/2020, 17h35

Pas de soucis. Quand tu fais une analyse de Fourier, ton spectre est symétrique. Tu as deux fois l'information. En faisant [0:int(N/2)] je ne récupère qu'une des deux parties.