在運用之前我們需要知道他是什么?是怎么來的?怎么去應用。
傅立葉變換是一種分析信號的方法,它可分析信號的組成成分,也可用這些成分合成信號。許多波形可作為信號的成分,比如正弦波、方波、鋸齒波等,傅立葉變換用正弦波作為信號的組成成分,在時域他們是相互重疊在一起的,我們需要運用傅里葉變換把他們分開並在頻域顯示出來。
連續傅里葉變換(Fourier Transform)如下:
連續傅里葉變換的反變換為:
滿足傅里葉變換的條件是f(t)在整個定義域是絕對可積的(不發散),只有這樣積分才有效。
快速傅里葉變換歸屬於離散傅里葉變換,理論部分大家可以自行百度。
那下面我們直接上代碼:
1 import numpy as np 2 import time as tm 3 import math 4 #******************************************** 5 #這個函數的作用是輸入轉換點數,然后求出多少次方 6 #******************************************** 7 def __fft_inpute(x): 8 for i in range(11): 9 if(2**i == x): 10 return i 11 return -1 12 #************************ 13 #這個函數是把輸入倒序排列 14 #************************ 15 def __fft_daoxu(x): 16 if(x == -1): 17 return -1 18 daoxu = [] 19 for i in range(2**x): 20 d = 0 21 for j in range(x): 22 if i & 1 == 1: 23 i >>= 1 24 d <<= 1 25 d += 1 26 else: 27 i >>= 1 28 d <<= 1 29 daoxu.append(d) 30 return daoxu 31 #**************************** 32 #快速傅里葉運算 基於時域的分解 33 #參數: 34 # data:是輸入數據,以列表的形式傳入 35 # N:輸入鏈表data的長度 36 #**************************** 37 def fft(data,N): 38 _data =[] 39 _post_n = __fft_inpute(N)#判斷運算級數 40 if(_post_n == -1): 41 return -1 42 43 44 daoxu_post = __fft_daoxu(_post_n) 45 print("原序列的倒置:",daoxu_post) 46 for i in daoxu_post: 47 _data.append(data[i]) 48 49 50 #蝶形運算,按時間抽選(DIT)的基-2 fft算法(庫里-圖基算法) 51 for i in range(_post_n):#第一次循環把就是各級運算分開 52 step = 1 << (i+1)#確定運算幾大步, 53 54 Wn = 1.0 + 0.0j; #給定第一次旋轉因子的值 55 56 print("第%d次蝶形運算"%(i)) 57 58 59 for j in range(int(step/2)):#不相同的Wn的個數 60 for x1 in range(j,N,step):#把同級的相同的Wn全部運算完 61 x2 = x1 + int(step / 2) 62 print(x1,x2) 63 64 #下面是兩點傅里葉變換 65 tem_variable = 0 + 0j 66 tem_variable = (_data[x1].real + _data[x2].real*Wn.real - _data[x2].imag*Wn.imag) + (_data[x1].imag + _data[x2].real*Wn.imag + _data[x2].imag*Wn.real)*1j 67 _data[x2] = (_data[x1].real - _data[x2].real*Wn.real + _data[x2].imag*Wn.imag) + (_data[x1].imag - _data[x2].real*Wn.imag - _data[x2].imag*Wn.real)*1j 68 _data[x1] = tem_variable 69 70 71 Wn = (math.cos(2 * math.pi*(j +1) / (N >> (_post_n - (i + 1))))) - (math.sin(2 * math.pi*(j+1 )/ (N >> (_post_n - (i + 1)))))*1j#歐拉公式求旋轉因子 72 73 74 return _data 75 76 num =[i for i in range(8)] 77 print("參與傅里葉變換的序列:",num) 78 star = tm.perf_counter() 79 print("我的fft:",fft(num,len(num)),len(num)) 80 end1 = tm.perf_counter() 81 print("numpy的fft:",np.fft.fft(num)) 82 end2 = tm.perf_counter() 83 print("自帶庫fft的運行時間:",end2 - end1,"我的fft的運行時間:",end1 - star)
下面我講解下旋轉因子:
Wn = (math.cos(2 * math.pi*(j +1) / (N >> (_post_n - (i + 1))))) - (math.sin(2 * math.pi*(j+1 )/ (N >> (_post_n - (i + 1)))))*1j#歐拉公式求旋轉因子
旋轉因子為: ,所以Wnk經過歐拉公式為:
。
首先我們看一張圖這個旋轉因子是怎么變化的:
是不是恍然大悟,N的變化是從2開始,往右依次增加,增加的倍數是以前的兩倍。所以我們可以寫成如下形式變化:
1 (N >> (_post_n - (i + 1))
N的右移一次表現的效果就是除以一次2。
接下來看看我們的運行結果:
1 參與傅里葉變換的序列: [0, 1, 2, 3, 4, 5, 6, 7] 2 原序列的倒置: [0, 4, 2, 6, 1, 5, 3, 7] 3 第0次蝶形運算 4 0 1 5 2 3 6 4 5 7 6 7 8 第1次蝶形運算 9 0 2 10 4 6 11 1 3 12 5 7 13 第2次蝶形運算 14 0 4 15 1 5 16 2 6 17 3 7 18 我的fft: [(28+0j), (-3.9999999999999996+9.65685424949238j), (-4+4j), (-4+1.6568542494923797j), (-4+0j), (-4-1.6568542494923806j), (-3.9999999999999996-4j), (-3.9999999999999987-9.65685424949238j)] 8 19 numpy的fft: [28.+0.j -4.+9.65685425j -4.+4.j -4.+1.65685425j 20 -4.+0.j -4.-1.65685425j -4.-4.j -4.-9.65685425j] 21 自帶庫fft的運行時間: 0.00610923200000002 我的fft的運行時間: 0.003303109999999998
經過測試,長序列的傅里葉變換numpy自帶的fft要快得多,這里只是驗證程序的可行性。