Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Calcul scientifique Python Discussion :

Intégration de spectre et écart-type


Sujet :

Calcul scientifique Python

  1. #1
    Futur Membre du Club
    Intégration de spectre et écart-type
    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

  2. #2
    Membre chevronné
    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.

  3. #3
    Futur Membre du Club
    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

  4. #4
    Membre chevronné
    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).

  5. #5
    Futur Membre du Club
    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

  6. #6
    Membre éprouvé
    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

  7. #7
    Futur Membre du Club
    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ù.

  8. #8
    Membre éprouvé
    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

  9. #9
    Futur Membre du Club
    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

  10. #10
    Membre éprouvé
    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.

  11. #11
    Futur Membre du Club
    Ah oui ça me rappelle de vieux souvenirs. Merci beaucoup pour tes différentes réponses

  12. #12
    Rédacteur



    int(N/2) donne la division entière et peut être (j'ai testé) remplacée par N//2 en python 3.