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