[PYTHON] FFT & Trendentfernung mit Numpy

Einführung

Ich möchte vorstellen, wie man Zeitreihendaten mit Numpy in Python FFT (Fast Fourier Transform) verwendet und wie man den Trend von Zeitreihendaten entfernt. FFT ist eine Berechnungsmethode, die DFT (Discrete Fourier Transform) mit hoher Geschwindigkeit verarbeitet. Dieser Artikel geht nicht auf die Theorie ein und zeigt den Mindestcode für die Durchführung der FFT. Das Referenzdokument lautet "Spektrumanalyse: Mikio Hino (Asakura Shoten)". Von den Grundlagen der Fourier-Analyse bis zur Theorie der FFT reicht dieses eine Buch aus.

Inhaltsverzeichnis

1. Zeitreihendaten
2. FFT-Ausführung
3. Trendentfernung
4. Frequenzglättung der Fourier-Komponente (Glättung) 5.Appendix

Hauptthema

1. Zeitreihendaten

30-Minuten-Daten mit einer Abtastfrequenz von 10 Hz. Die orange Linie ist der gleitende Durchschnitt. Sie sehen, dass es einen Trend gibt. Abbildung 1. Ursprüngliche Zeitreihendaten

Wir werden diese Daten FFT.

2. FFT-Ausführung

N =len(X)      #Datenlänge
fs=10          #Abtastfrequenz
dt =1/fs       #Abtastintervall
t = np.arange(0.0, N*dt, dt) #Zeitachse
freq = np.linspace(0, fs,N) #Frequenzachse
fn=1/dt/2     #Nyquist-Frequenz

FFT ist eine Berechnungsmethode, die die schnellste Berechnungsgeschwindigkeit liefert, wenn die Datenlänge eine Potenz von 2 ist. Sie kann jedoch auch dann berechnet werden, wenn dies nicht der Fall ist (obwohl die Verarbeitungszeit etwas länger ist). Wenn die Datenlänge jedoch eine Primzahl ist, dauert die Verarbeitung im Vergleich zur Potenz von 2 relativ lange. Daher scheint es besser, 0 aufzufüllen, damit sie zur Potenz von 2 wird. Ich habe den Code zur Bestätigung im Anhang veröffentlicht. Bitte überprüfen Sie ihn.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0 #Schnitt nach Nyquist-Frequenz

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

In np.fft.fft ist es durch die komplexe Fourier-Komponente gegeben, daher zeigt die folgende Abbildung den absoluten Wert. Die Trendkomponente liegt bei 0 Hz. Lassen Sie uns im nächsten Abschnitt den Trend entfernen. Figur 2. FFT-Ausführung für die ursprünglichen Zeitreihendaten

3. Trendentfernung

Betrachtet man den gleitenden Durchschnitt (orange Linie) in Fig. 1, so ist ersichtlich, dass die Achse der Zeitreihendaten von der x-Achse abweicht und es einen Trend in den Daten gibt. Entfernen Sie den Trend, indem Sie mit der folgenden Operation 0,03 Hz oder weniger auf 0 setzen.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0
F[(freq<=0.03)]=0 #0.Unter 03HZ entfernt
X_1=np.real(np.fft.ifft(F))*N

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

Figur 3. Zeitreihendaten nach Trendentfernung

Sie können sehen, dass es sich um eine periodische Funktion handelt, die von der x-Achse ausgeht.

4. Frequenzglättung der Fourier-Komponente (Glättung)

Die FFT der Daten in Fig. 3 ergibt Fig. 4. Figur 4. FFT für Zeitreihendaten nach Trendentfernung

Da es klappert, werde ich versuchen, es durch Anwenden eines Glättungsfensters zu glätten.

window=np.ones(5)/5 #Fenster glätten
F3=np.convolve(F2,window,mode='same') #Falten
F3=np.convolve(F3,window,mode='same') #Falten
F3=np.convolve(F3,window,mode='same') #Falten

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

Abbildung 5. Geglättet

5. Appendix Bereiten Sie drei Arten von Datenlängen vor und vergleichen Sie die Berechnungszeiten. ① 2 ^ 19 (Potenz von 2) ② 2 ^ (19) -1 (Primzahl) ③ 2 ^ (19) -2 (weder Primzahl noch Potenz von 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]")

Berechnungsergebnis ①0.04197[sec] ②0.1679[sec] ③0.05796[sec]

Wenn die Datenlänge eine Primzahl ist, ist es besser, 0 aufzufüllen und eine Potenz von 2 zu machen.