[PYTHON] FFT & trend removal using Numpy

Introduction

In Python, I would like to introduce how to use Numpy to perform FFT (Fast Fourier Transform) on time series data and how to remove trends in time series data. FFT is a calculation method that processes DFT (Discrete Fourier Transform) at high speed. This article does not touch on theory and shows the minimum code to perform an FFT. The reference document is "Spectrum Analysis: Mikio Hino (Asakura Shoten)". From the basics of Fourier analysis to the theory of FFT, this one book is enough.

table of contents

1. Time series data
2. FFT execution
3. Trend removal
4. Frequency smoothing of Fourier component (smoothing) 5.Appendix

Main subject

1. Time series data

30-minute data with a sampling frequency of 10 Hz. The orange line is the moving average. You can see that there is a trend. Figure 1. Original time series data

We will FFT this data.

2. FFT execution

N =len(X)      #Data length
fs=10          #Sampling frequency
dt =1/fs       #Sampling interval
t = np.arange(0.0, N*dt, dt) #Time axis
freq = np.linspace(0, fs,N) #Frequency axis
fn=1/dt/2     #Nyquist frequency

FFT is a calculation method that gives the fastest calculation speed when the data length is a power of 2, but it can be calculated even if it is not (although the processing time will be slightly longer). However, when the data length is a prime number, it takes a relatively long time to process compared to the power of 2, so it seems better to pad 0 so that it becomes a power of 2. I have posted the code to confirm it in the Appendix, so please check it out.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0 #Cut after Nyquist frequency

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, it is given by the complex Fourier component, so the figure below plots the absolute value. The trend component is around 0Hz. In the next section, let's remove the trend. Figure 2. Perform FFT on original time series data

3. Trend removal

Looking at the moving average (orange line) in Fig. 1, you can see that the axis of the time series data deviates from the x axis and there is a trend in the data. Let's remove the trend by setting 0.03Hz or less to 0 by the following operation.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0
F[(freq<=0.03)]=0 #0.Removed below 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. Time series data after trend removal

You can see that it is a periodic function starting from the x-axis.

4. Frequency smoothing of Fourier component (smoothing)

FFT the data in Fig. 3 gives Fig. 4. Figure 4. FFT for time series data after trend removal

Since it is rattling, I will try smoothing by applying a smoothing window.

window=np.ones(5)/5 #Smoothing window
F3=np.convolve(F2,window,mode='same') #Convolution
F3=np.convolve(F3,window,mode='same') #Convolution
F3=np.convolve(F3,window,mode='same') #Convolution

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

Figure 5. Smoothed

5. Appendix Prepare three types of data lengths and compare the calculation times. ① 2 ^ 19 (power of 2) ② 2 ^ (19) -1 (prime number) ③ 2 ^ (19) -2 (neither a prime number nor a power of 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]")

Calculation result ①0.04197[sec] ②0.1679[sec] ③0.05796[sec]

If the data length is a prime number, it seems better to do 0 padding and make it a power of 2.