用Python實現DFT並繪製功率譜

知識背景：傅利葉變換可以分為連續傅利葉變化和離散傅利葉變換，分別是ft,fs,dtft,dtfs。其中dtft是我們常說的離散時間傅利葉變換，但這種變換並不一定能夠由計算機進行處理，因為對於非週期訊號來說其譜一般是連續譜，這樣就無法由計算機完成了。所以dft就出現了，我們知道dtft是以2π為週期的，我們一般只需要取其主值（可以看作取乙個完整的週期）進行分析，如果對dtft在0到2π內均勻盡行取樣的話得到的結果就一定是離散的，如果我們的取樣是遵循一定規則的那就可以用取樣後的譜完整的恢復原訊號。

完整的過程是：原訊號遵循取樣定理取樣得到離散時間訊號序列，假設其序列長度為n，如果我們的dft的點數為k，只要k>=n，那麼我們就可以完全重構其頻譜（也就是要恢復其dtft的譜），再由恢復的譜就可以恢復出原訊號。再到後來對dft的計算方式進行改進就有了fft，是按照dft具有週期性對稱性把多點的dft分為點數最少的dft來進行計算。

fft的演算法實現比較複雜，我們一般都是呼叫第三方庫中已經實現的通用的fft，但是dft這種不關心運算效率的演算法實現起來就很簡單，無非就是兩層迴圈。

外層迴圈就是對k的遍歷，內層迴圈就是對n的遍歷，相乘，求和。

比如我們要對

這樣乙個訊號做dft，f1=0.1 f2=0.3，實現過程如下：

from math import *
import numpy as np
import matplotlib.pyplot as plt
#懶得導包了，直接全抓進來
def signal
(n):
return
(sin
(0.2
*pi*n+pi/3)
+10*sin
(0.6
*pi*n+pi/4)
)#生成wn項
def wn_k
(k,n,n)
:return
complex
(cos(2
*pi*n*k/n)
,sin(-
2*pi*n*k/n)
)amplitude=
#準備乙個空列表
power_spectrum=
sums=
0#256點dft，x(0
)到x(
255)
for k in range(0
,256):
for n in range(1
,257):
#n的取值為從1到256
sums=sums+
signal
(n)*
wn_k
(k,n,
256)
amplitude.
(sums)
sums=
0print
(amplitude,
len(amplitude)
)for i in range(0
,256):
power_spectrum.
(amplitude[i]**
2)plt.
suptitle
("xn=sin(2*pi*f1*n+pi/3)+10*sin(2*pi*f2*n+pi/4)"
)plt.
subplot(2
,1,1
)plt.
plot
(np.
abs(amplitude)
)plt.
title
("amplitude_spectrum"
)plt.
subplot(2
,1,2
)plt.
plot
(np.
abs(power_spectrum)
)plt.
title
("power_spectrum"
)plt.
show
()

其中那一段for迴圈的結果就是把變換結果全部放入了乙個列表，列表裡儲存的是複數。之所以要把結果放入列表是因為matplotlib的繪圖函式的入口引數是列表型別。我們要繪製幅度譜的話就取列表裡每個復數值的模值，繪製功率譜的話就把每個值平方一下。為了驗證我們所求結果是否正確，我們可以呼叫numpy庫的fft來對照一下

庫函式版本：

import numpy as np
import matplotlib.pyplot as plt
from math import *
power_spectrum=
x=np.
linspace(1
,256
,256
)y=np.
sin(
0.2*np.pi*x+np.pi/3)
+10*np.
sin(
0.6*np.pi*x+np.pi/4)
#plt.subplot(3,1,1)
#plt.plot(y)
#plt.title("input xn")
trans=np.fft.
fft(y,n=
256)
plt.
suptitle
("xn=sin(2*pi*f1*n+pi/3)+10*sin(2*pi*f2*n+pi/4)"
)plt.
subplot(2
,1,1
)plt.
plot
(np.
abs(trans)
)plt.
title
("amplitude by numpy"
)for x in range(0
,256):
power_spectrum.
(trans[x]**
2)plt.
subplot(2
,1,2
)plt.
plot
(np.
abs(power_spectrum)
)plt.
title
("power_spectrum by numpy"
)plt.
show
()

對比結果如圖：結果一致

請注意：實際上是離散譜，但繪圖時自動把點與點之間用線段給連線起來了，改用stem繪圖就是離散圖了

用Python實現DFT並繪製功率譜

DFT 理解與python實現

用python繪製詞云

用python繪製玫瑰花

用Python實現DFT並繪製功率譜

DFT 理解與python實現

用python繪製詞云

用python繪製玫瑰花

相關推薦