淺談快速離散傅里葉變換的實現


在運用之前我們需要知道他是什么?是怎么來的?怎么去應用。

傅立葉變換是一種分析信號的方法,它可分析信號的組成成分,也可用這些成分合成信號。許多波形可作為信號的成分,比如正弦波、方波、鋸齒波等,傅立葉變換用正弦波作為信號的組成成分,在時域他們是相互重疊在一起的,我們需要運用傅里葉變換把他們分開並在頻域顯示出來。

連續傅里葉變換(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要快得多,這里只是驗證程序的可行性。


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM