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

**nat8546** · 16/05/2017, 15h51

Bonjour
Je souhaite réaliser une transformé de Fourier.
Mon jeu de donnée est le suivant:
test.txt

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 :

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

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.

**Julien N** · 17/05/2017, 10h50

Salut,

Je te propose ceci:

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
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

**nat8546** · 17/05/2017, 11h24

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 !

**Julien N** · 17/05/2017, 14h18

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 :

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

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

**marco056** · 17/05/2017, 14h26

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.

**Julien N** · 17/05/2017, 14h47

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

**nat8546** · 17/05/2017, 14h56

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 :

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

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 !

**nat8546** · 17/05/2017, 15h21

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.

**Julien N** · 17/05/2017, 15h29

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.

**nat8546** · 17/05/2017, 15h58

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

**Julien N** · 18/05/2017, 09h51

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

**nat8546** · 18/05/2017, 10h54

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 !

[Transformée de Fourier][Numpy][Scipy] [Python 3.X]

Python

Discussions similaires

Partager

Partager