时间序列分析之:傅里叶变换找周期
時間序列分析
萬萬沒想到吧,信號處理的技術,能用在數據分析中。誰叫我是學通信出生的呢?
承接上一篇:函數分解
本節承接上文找函數的周期。
文章目錄
- 時間序列分析
- 傅里葉變換
- 一、傅里葉變換(FFT)是什么?
- 二、使用步驟
- 1.新建FFT函數
- 2.測試函數
- 總結
傅里葉變換
通信專業的我,看到找周期時,不由自主想起了傅里葉變換。
傅里葉變換就是把時域上的信號,變換到頻域上,用很多個正弦波來合成時域信號。所以,我們找信號幅度最大的那個正弦波的頻率,作為函數的周期。
傅里葉變換最詳細的介紹見這個文:詳細得令人發指
一、傅里葉變換(FFT)是什么?
世間萬物都在隨著時間不停的改變,并且永遠不會靜止下來。比如一個周期變化的函數,但如果我告訴你,用另一種方法來觀察世界,你會發現世界是永恒不變的,這個靜止的世界就叫做頻域。
比如一年四季每天都在變化,但是變化的頻率卻永遠都是12個月。所以在頻率的視角看,四季更迭就是固定的12個月。
我們知道經濟有周期,那么周期是多少呢?傅里葉變換讓我們站在頻率的角度看,就更加直觀。
二、使用步驟
1.新建FFT函數
功能:把函數進行傅里葉變換,變換到頻域,以期獲得函數的周期。
然后我們只用找到局部增幅最大的那個正弦波作為周期。
比如下圖的上半部分是我們觀察到的實際信號,下半部分進行了傅里葉變換,得到了他的頻率為10,相對于上面一直變化的曲線,下面頻率的圖,只有在10的位置有一個波動點,其余都是0,所以是相對靜止的。我們可以說,上半部分的圖形的頻率為10,根據總計的1000個輸入點,因此這個圖形的周期是100,也就是100個采樣點就是一個完整周期。:
新建的fft代碼如下(示例):
# 功能:把函數進行傅里葉變換,變換到頻域,以期獲得函數的周期 # 輸入:時間序列,獲取頻率點數值n(可選),頻率對應幅度的下限值fmin(可選) # 輸入序列的X軸需要歸一化為1 # 輸出: n個序列的下標以及對應的幅度值 # 創建時間: 2021-1-26import pandas as pd import numpy as np import math import numpy as np from scipy.fftpack import fft, ifft import matplotlib.pyplot as plt import seaborn import scipy.signal as signaldef fftTransfer(timeseries, n=10, fmin=0.2):yf = abs(fft(timeseries)) # 取絕對值yfnormlize = yf / len(timeseries) # 歸一化處理yfhalf = yfnormlize[range(int(len(timeseries) / 2))] # 由于對稱性,只取一半區間yfhalf = yfhalf * 2 # y 歸一化xf = np.arange(len(timeseries)) # 頻率xhalf = xf[range(int(len(timeseries) / 2))] # 取一半區間plt.subplot(211)x = np.arange(len(timeseries)) # x軸plt.plot(x, timeseries)plt.title('Original wave')plt.subplot(212)plt.plot(xhalf, yfhalf, 'r')plt.title('FFT of Mixed wave(half side frequency range)', fontsize=10, color='#7A378B') # 注意這里的顏色可以查詢顏色代碼表fwbest = yfhalf[signal.argrelextrema(yfhalf, np.greater)]xwbest = signal.argrelextrema(yfhalf, np.greater)plt.plot(xwbest[0][:n], fwbest[:n], 'o', c='yellow')plt.show(block=False)plt.show()xorder = np.argsort(-fwbest) # 對獲取到的極值進行降序排序,也就是頻率越接近,越排前print('xorder = ', xorder)print(type(xorder))xworder = list()xworder.append(xwbest[x] for x in xorder) # 返回頻率從大到小的極值順序fworder = list()fworder.append(fwbest[x] for x in xorder) # 返回幅度if len(fwbest) <= n:fwbest = fwbest[fwbest >= fmin].copy()return len(timeseries)/xwbest[0][:len(fwbest)], fwbest #轉化為周期輸出else:fwbest = fwbest[fwbest >= fmin].copy()print(len(fwbest))print(xwbest)return len(timeseries)/xwbest[0][:len(fwbest)], fwbest # 只返回前n個數 #轉化為周期輸出2.測試函數
簡單的測試就是新建一個模擬信號,然后fft變換,核心代碼如下(生成sin(20πx)的信號,對應的頻率為10):
xtime = np.arange(0,1000,1) xnorm = xtime/len(xtime) queshi = xnorm +1 plt.plot(xtime,queshi) zouqi = [sin(x*20*math.pi) for x in xnorm] plt.plot(xtime,zouqi) signal = zouqi y= signal y = pd.Series(y).astype('float') df_price=y x , y=fftTransfer(df_price,n=5,fmin = 0.015) #快速傅里葉變換 print('x = ',x) print('y = ',y)但是,實際的生活中的信號就像愛情一樣復雜,不是簡單的一個周期,而是多個周期的疊加。
比如我們產生一個頻率逐漸加快的函數,
他的周期是多少呢?
首先fft變換到頻域:
選取增幅最大的5個點,作為周期。
所以可以看出5個周期不同的函數合成。
測試代碼如下:
這個周期就是我們想要的,在函數分解的時候,需要輸入的周期。
總結
通過fft變換確實可以提取到一些周期信息,配合函數分解來使用。
后面我們用創業板指數來試探在現實中是否有意義。
總結
以上是生活随笔為你收集整理的时间序列分析之:傅里叶变换找周期的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 时间序列分析之:函数分解decompos
- 下一篇: 时间序列分析源资料汇总