[PYTHON] Suppression de la FFT et des tendances avec Numpy

introduction

Je voudrais vous présenter comment FFT (Fast Fourier Transform) des données de séries chronologiques à l'aide de Numpy en Python et comment supprimer la tendance des données de séries chronologiques. La FFT est une méthode de calcul qui traite la DFT (Discrete Fourier Transform) à grande vitesse. Cet article n'aborde pas la théorie et montre le code minimum pour effectuer la FFT. Le document de référence est "Analyse du spectre: Mikio Hino (Asakura Shoten)". Des bases de l'analyse de Fourier à la théorie de la FFT, ce livre suffit.

table des matières

1. Données de séries chronologiques
2. Exécution FFT
3. Suppression de la tendance
4. Lissage de fréquence de la composante de Fourier (lissage) 5.Appendix

Sujet principal

1. Données de séries chronologiques

Données de 30 minutes avec une fréquence d'échantillonnage de 10 Hz. La ligne orange est la moyenne mobile. Vous pouvez voir qu'il y a une tendance. Figure 1. Données de série chronologique d'origine

Nous allons FFT ces données.

2. Exécution FFT

N =len(X)      #Longueur des données
fs=10          #Fréquence d'échantillonnage
dt =1/fs       #Intervalle d'échantillonnage
t = np.arange(0.0, N*dt, dt) #Axe du temps
freq = np.linspace(0, fs,N) #Axe de fréquence
fn=1/dt/2     #Fréquence de Nyquist

La FFT est une méthode de calcul qui donne la vitesse de calcul la plus rapide lorsque la longueur des données est une puissance de 2, mais elle peut être calculée même si ce n'est pas le cas (bien que le temps de traitement soit légèrement plus long). Cependant, lorsque la longueur des données est un nombre premier, le traitement prend un temps relativement long par rapport à la puissance de 2, il semble donc préférable de compléter 0 pour qu'il devienne la puissance de 2. J'ai posté le code pour le confirmer en annexe, veuillez donc le vérifier.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0 #Couper après la fréquence de Nyquist

plt.plot(freq,np.abs(F))#
plt.xlabel("[Hz]")
plt.ylabel("Amp")
plt.xlim(-0.01,0.5)
plt.grid()
plt.show()

Dans np.fft.fft, il est donné par la composante complexe de Fourier, donc la figure ci-dessous trace la valeur absolue. La composante de tendance est d'environ 0 Hz. Dans la section suivante, supprimons la tendance. Figure 2. Exécution FFT sur les données de la série chronologique d'origine

3. Suppression de la tendance

En regardant la moyenne mobile (ligne orange) sur la figure 1, on peut voir que l'axe des données de la série temporelle s'écarte de l'axe des x et il y a une tendance dans les données. Supprimons la tendance en définissant 0,03 Hz ou moins à 0 par l'opération suivante.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0
F[(freq<=0.03)]=0 #0.Supprimé en dessous de 03HZ
X_1=np.real(np.fft.ifft(F))*N

plt.xlabel("Time [s]")
plt.ylabel("Signal")
plt.xlim(-50,1850)
plt.grid()
plt.show()

Figure 3. Données de séries chronologiques après suppression de la tendance

Vous pouvez voir qu'il s'agit d'une fonction périodique à partir de l'axe des x.

4. Lissage de fréquence de la composante de Fourier (lissage)

La FFT des données de la figure 3 donne la figure 4. Graphique 4. FFT pour les données de séries chronologiques après suppression de la tendance

Comme il claque, je vais essayer de le lisser en appliquant une fenêtre de lissage.

window=np.ones(5)/5 #Fenêtre de lissage
F3=np.convolve(F2,window,mode='same') #Plier
F3=np.convolve(F3,window,mode='same') #Plier
F3=np.convolve(F3,window,mode='same') #Plier

plt.plot(freq,np.abs(F3))
plt.xlabel("[Hz]")
plt.ylabel("Amp")
plt.xlim(-0.01,0.5)
plt.grid()
plt.show()

Figure 5. Lissé

5. Appendix Préparez trois types de longueurs de données et comparez les temps de calcul. ① 2 ^ 19 (puissance de 2) ② 2 ^ (19) -1 (premier) ③ 2 ^ (19) -2 (ni premier ni puissance de 2)

import time

if __name__ == '__main__':
    start = time.time()
    x2 = np.random.uniform(size=2**19-2)#2**19 , 2**19-1
    print(np.fft.fft(x2))
    elapsed_time = time.time() - start
    print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

Résultat du calcul ①0.04197[sec] ②0.1679[sec] ③0.05796[sec]

Si la longueur des données est un nombre premier, il semble préférable de faire un remplissage à 0 et d'en faire une puissance de 2.