Discussion :

[nat8546]

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]

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]

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]

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]

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]

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