[Transformée de Fourier][Numpy][Scipy]

Bonjour
Je souhaite réaliser une transformé de Fourier.
Mon jeu de donnée est le suivant:
Pièce jointe 276080

En colonne 0, le temps (ne commence pas à 0)
En colonne 2, l'élévation de surface

J'ai ouvert le fichier avec pandas, j'ai une dataFrame.
Je souhaite obtenir la transformé de Fourier afin de déterminer la période des vagues suite à l'élévation de surface.
Sauf que en passant par l'outil Scipy et Numpy, j'obtiens toujours une fréquence de 0 pour le pic de la transformée de fourier, ce qui est faux (on devrait être entre 5 et 8 secondes, les temps caractéristiques d'une houle de vent).

Quelqu'un saurait comment à partir du jeu de donnée, obtenir un graphique de la transformée de Fourier ?

PS : Voici mon code :

Code:

---

1 2 3 4 5 6 7 8 9

df=pd.read_csv(file,sep=';',header=None) df.columns=['Time','Nothing','HS','Other1','Other2'] N=len(df) T=df.loc[0,'Time']-df.loc[len(df)-1,'Time'] #Intervalle de temps freq=np.linspace(0.0, 1, N) #Vecteur fréquence   HSf=scipy.fft(df['HS']) #Transformée de Fourier plt.plot(freq,abs(HSf))

---

J'avais installé les modules avant.
Je comprend que l'erreur vient du fait le premier élément de HSf (le pic), est associé au premier élément de mon vecteur fréquence (0 du coup), mais je ne sais pas comment changer ça.

Merci de votre aide.

Salut,

Je te propose ceci:

Code:

---

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 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58

import numpy as np import pandas as pd import matplotlib.pyplot as plt   def fourier(t, signal): """Fast Fourier Transform of time series.       Parameters     ----------     t : array-like         Array of time values is seconds     signal : array-like         Array of the signal amplitudes. Must be of the same length as the         t serie.       Returns     -------     freq : numpy array         Array of frequencies in Hz.     ampl : numpy array         Array of amplitudes in [unit]^2.s.     """ N, dt = len(t), t[1]-t[0] signal = signal - np.mean(signal) dft = np.fft.fft(signal) freq = np.fft.fftfreq(t.shape[-1])[0:int(N/2)] 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   # Import data df = pd.read_csv('test.txt', sep=';', header=None, usecols=(0, 2)) df.columns = ['Time', 'Elevation']   # Perform spectral analysis freq, ampl = fourier(df.Time.values, df.Elevation.values)   # Quick checks to see if the spectral analysis seems to produce expected # results (assumption: wave elevation following a Rayleigh distribution) m0 = np.trapz(ampl * freq**0, freq) m2 = np.trapz(ampl * freq**2, freq) wave_hs = 4 * np.sqrt(m0) wave_tz = np.sqrt(m0 / m2)   # Some plots fig, axes = plt.subplots(2)   axes[0].plot(df.Time, df.Elevation) axes[0].grid(True) axes[0].set_xlabel('Time (s)') axes[0].set_ylabel('Elevation (m)')   axes[1].plot(freq, ampl) axes[1].grid(True) axes[1].set_xlabel('Frequencies (Hz)') axes[1].set_ylabel('Amplitudes ()')   plt.show()

---

Attention, je ne garanti pas l'exactitude des résultats. D'une fois sur l'autre je n'arrive pas à me souvenir s'il faut normaliser le spectre par le nombre d'échantillon, et s'il faut refaire la symétrie avant de faire les intégrations qui te seront nécessaire pour en déduire quelque chose.

A noter, le signal temporelle de l'élévation de la surface libre ne correspond pas au Hs. Le Hs est une donnée statistique à calculer. Ici je trouve 0.4m, ce qui me fait penser que j'aurais peut-être du refaire la symétrie du spectre (car Hs == Hauteur, donc double amplitude et vue le tracé du signal on devrait être plus proche de 0.8m...).

EDIT: Correction de la fonction fourier pour sortir non pas le spectre des amplitude mais le spectre de l'énergie. Et puis le calcul de la période était faux (pas de 2*pi)

Julien

Ok merci !
J'ai trouvé une autre solution assez proche de la votre pour déterminer la période caractéristique.
J'ai effectivement un petit souci avec le Hs, je vais chercher voir où peut être l'erreur (on m'a donné en résultat Hs, et ça ne colle pas).

A bientôt !

Je viens d'éditer mon post précédent. Il s'avère que je faisais une erreur dans le calcul du spectre. Il manquait un carré dans le calcul de l'amplitude afin de sortir une amplitude homogène à une énergie. Une technique toute simple pour vérifier son spectre est de calculer le moment d'ordre 0:

Code:

---

m0 = np.trapz(ampl * freq**0, freq)

---

Et de comparer sa racine avec l'écart type du signal temporel. Le théorème de Parceval dit qu'on devrait trouver la même chose (l'écart type en spectral doit être identique à l'écart type en temporel). Ce qui est désormais le cas. Je trouve donc un Hs de1.96m (donc double amplitude) et un Tz de 12.74s.

Bonne journée,

Julien

Une petite question :
comment peux-tu trouver des amplitudes de 10 ou 12 avec une élévation de 1.5 ?
Je précise que je n'ai pas regardé en détail le programme mais cela m'interpelle.

Citation:

---

Envoyé par marco056

Une petite question :
comment peux-tu trouver des amplitudes de 10 ou 12 avec une élévation de 1.5 ?
Je précise que je n'ai pas regardé en détail le programme mais cela m'interpelle.

---

Salut Marco.

Je trouve une amplitude de 1.96m (hauteur significative, correspondant grossièrement à la moyenne du tier des vagues les plus hautes, aussi appelé H1/3). 12s c'est la période entre deux vagues, ou plutôt période entre deux croisements successifs à amplitude=0m vers le haut (ou vers le bas, pareil). Si c'est bien là ton interrogation.

J

Salut Julien,
Merci pour l'édit !
Effectivement en utilisant ta première fonction, je n'avais pas trop les mêmes résultats avec ceux attendus et j'avais trouvé une manière de faire pour trouver Tp :

Code:

---

1 2 3 4 5

dft = np.fft.fft(signal) freq = np.fft.fftfreq(N,dt) #dt est le pas de temps du signal idx=np.argmax(np.abs(dft)) #Donne l'index auquel on a le max de la transformé de fourier freq[idx] #Donne la fréquence associé au max de fourier

---

Avec cette méthode j'obtenais un Tp qui était déjà plus en accord. Mais j'avais la symétrie du coup je pouvais avoir des Tp négatifs (bon suffit juste de prendre la valeur absolue.
J'étais en train de chercher à trouver Hs mais de manière plus laborieuse (passer par l'histogramme des hauteurs de vagues et faire une courbe de densité par dessus), je vais regarder la nouvelle façon.

Merci !

En fait Julien, je reviens sur ta période :
La période que tu obtiens est la période entre deux 0.
Or la période caractéristique de la vague correspond à la fréquence pour laquelle on a max(S(f)).
Avec S le spectre en m²s
Du coup je pense je vais garder ma méthode pour Tp (j'espère qu'elle est bonne).

Par contre pour Hs ça m'a l'air bon ça correspond plutôt aux données résultats que j'ai.

Quelle valeur de Tp tu penses être juste?

J'ai calculé Tz mais pour Tp il faudrait faire np.sqrt(m0 / m4). J'avoue que je trouve une valeur un peu trop élevée au vue du signal, mais je ne vois pas de loup. Comme la calcul-tu? En prenant juste la période du spectre correspondant au pic le plus élevé?

EDIT: Bon j'avais pas les yeux en face des trous. Mon spectre sortais en m2 au lieu de m2.s... Bref, maintenant j'ai un Tz de 4.98s, ce qui cest bien mieux.

Oui je pense utiliser dft (cf deux-trois messages plus haut).
Avec ce que j'ai écrit j'obtiens une période de 9.1 s

Je trouve aussi une valeur aberrante avec la formule sqrt(m0/m4) : 78 s

Citation:

---

Envoyé par nat8546

Je trouve aussi une valeur aberrante avec la formule sqrt(m0/m4) : 78 s

---

Oui, je ne retrouve plus la référence à cette formule, mais il me semblait qu'elle était valide pour une distribution de Rayleigh. Mais ça n'a pas l'air d'être juste. Normalement le Tz est juste un peu plus faible que le Tp (de 2-3s). Ta méthode est bonne pour choper le Tp. Tu peux avec le Tz et le Tp essayer d'évaluer le gamma. Ou alors comme c'est du wind sea, utiliser une approxiamtion en fonction du Hs et tu Tp.

J

En regardant de plus près avec les résultats que je devrais obtenir, j'ai quelques écarts, parfois grossier.
La méthode est assez bonne mais on m'a dit que je devrais utiliser une fenêtre de Hann sur mon échantillon, permettant d'éviter la discontinuité aux extrémités.

Je vais m'essayer à cette méthode ajd.

Merci pour l'aide !