您好,登錄后才能下訂單哦!
這篇文章給大家分享的是有關(guān)使用python如何實(shí)現(xiàn)DFT的內(nèi)容。小編覺(jué)得挺實(shí)用的,因此分享給大家做個(gè)參考,一起跟隨小編過(guò)來(lái)看看吧。
DFT
DFT(Discrete Fourier Transform),離散傅里葉變化,可以將離散信號(hào)變換到頻域,它的公式非常簡(jiǎn)單:
離散頻率下標(biāo)為k時(shí)的頻率大小
離散時(shí)域信號(hào)序列
信號(hào)序列的長(zhǎng)度,也就是采樣的個(gè)數(shù)
如果你剛接觸DFT,并且之前沒(méi)有信號(hào)處理的相關(guān)經(jīng)驗(yàn),那么第一次看到這個(gè)公式,你可能有一些疑惑,為什么這個(gè)公式就能進(jìn)行時(shí)域與頻域之間的轉(zhuǎn)換呢?
這里,我不打算去解釋它,因?yàn)槲宜接邢蓿f(shuō)的不清楚。相反,在這里我想介紹,作為一個(gè)程序員,如何如實(shí)現(xiàn)DFT
從矩陣的角度看DFT
DFT的公式,雖然簡(jiǎn)單,但是理解起來(lái)比較麻煩,我發(fā)現(xiàn)如果用矩陣相乘的角度來(lái)理解上面的公式,就會(huì)非常簡(jiǎn)單,直接上矩陣:
OK,通過(guò)上面的表示,我們很容易將DFT理解成為一種矩陣相乘的操作,這對(duì)于我們編碼是很容易的。
Talk is cheap, show me the code
根據(jù)上面的理解,我們只需要構(gòu)建出S SS矩陣,然后做矩陣相乘,就等得到DFT的結(jié)果
在這之前,我們先介紹如何生成正弦信號(hào),以及如何用scipy中的fft模塊進(jìn)行DFT操作,以驗(yàn)證我們的結(jié)果是否正確
正弦信號(hào)
A: 幅度
f: 信號(hào)頻率
n: 時(shí)間下標(biāo)
T: 采樣間隔, 等于 1/fs,fs為采樣頻率
? \phi?: 相位
下面介紹如何生成正弦信號(hào)
import numpy as np import matplotlib.pyplot as plt %matplotlib inline
def generate_sinusoid(N, A, f0, fs, phi): ''' N(int) : number of samples A(float) : amplitude f0(float): frequency in Hz fs(float): sample rate phi(float): initial phase return x (numpy array): sinusoid signal which lenght is M ''' T = 1/fs n = np.arange(N) # [0,1,..., N-1] x = A * np.cos( 2*f0*np.pi*n*T + phi ) return x N = 511 A = 0.8 f0 = 440 fs = 44100 phi = 0 x = generate_sinusoid(N, A, f0, fs, phi) plt.plot(x) plt.show()
# 另一種生成正弦信號(hào)的方法,生成時(shí)長(zhǎng)為t的序列 def generate_sinusoid_2(t, A, f0, fs, phi): ''' t (float) : 生成序列的時(shí)長(zhǎng) A (float) : amplitude f0 (float) : frequency fs (float) : sample rate phi(float) : initial phase returns x (numpy array): sinusoid signal sequence ''' T = 1.0/fs N = t / T return generate_sinusoid(N, A, f0, fs, phi) A = 1.0 f0 = 440 fs = 44100 phi = 0 t = 0.02 x = generate_sinusoid_2(t, A, f0, fs, phi) n = np.arange(0, 0.02, 1/fs) plt.plot(n, x)
Scipy FFT
介紹如何Scipy的FFT模塊計(jì)算DFT
注意,理論上輸入信號(hào)的長(zhǎng)度必須是才能做FFT,而scipy中FFT卻沒(méi)有這樣的限制
這是因?yàn)楫?dāng)長(zhǎng)度不等于時(shí),scipy fft默認(rèn)做DFT
from scipy.fftpack import fft # generate sinusoid N = 511 A = 0.8 f0 = 440 fs = 44100 phi = 1.0 x = generate_sinusoid(N, A, f0, fs, phi) # fft is X = fft(x) mX = np.abs(X) # magnitude pX = np.angle(X) # phase # plot the magnitude and phase plt.subplot(2,1,1) plt.plot(mX) plt.subplot(2,1,2) plt.plot(pX) plt.show()
自己實(shí)現(xiàn)DFT
自己實(shí)現(xiàn)DFT的關(guān)鍵就是構(gòu)造出S,有兩種方式:
def generate_complex_sinusoid(k, N): ''' k (int): frequency index N (int): length of complex sinusoid in samples returns c_sin (numpy array): the generated complex sinusoid (length N) ''' n = np.arange(N) c_sin = np.exp(1j * 2 * np.pi * k * n / N) return np.conjugate(c_sin) def generate_complex_sinusoid_matrix(N): ''' N (int): length of complex sinusoid in samples returns c_sin_matrix (numpy array): the generated complex sinusoid (length N) ''' n = np.arange(N) n = np.expand_dims(n, axis=1) # 擴(kuò)充維度,將1D向量,轉(zhuǎn)為2D矩陣,方便后面的矩陣相乘 k = n m = n.T * k / N # [N,1] * [1, N] = [N,N] S = np.exp(1j * 2 * np.pi * m) # 計(jì)算矩陣 S return np.conjugate(S)
# 生成信號(hào),用于測(cè)試 N = 511 A = 0.8 f0 = 440 fs = 44100 phi = 1.0 x = generate_sinusoid(N, A, f0, fs, phi) # 第一種方式計(jì)算DFT X_1 = np.array([]) for k in range(N): s = generate_complex_sinusoid(k, N) X_1 = np.append(X_1, np.sum(x * s)) mX = np.abs(X_1) pX = np.angle(X_1) # plot the magnitude and phase plt.subplot(2,1,1) plt.plot(mX) plt.subplot(2,1,2) plt.plot(pX) plt.show() # 結(jié)果和scipy的結(jié)果基本相同
# 第二種方法計(jì)算DFT S = generate_complex_sinusoid_matrix(N) X_2 = np.dot(S, x) mX = np.abs(X_2) pX = np.angle(X_2) # plot the magnitude and phase plt.subplot(2,1,1) plt.plot(mX) plt.subplot(2,1,2) plt.plot(pX) plt.show()
感謝各位的閱讀!關(guān)于“使用python如何實(shí)現(xiàn)DFT”這篇文章就分享到這里了,希望以上內(nèi)容可以對(duì)大家有一定的幫助,讓大家可以學(xué)到更多知識(shí),如果覺(jué)得文章不錯(cuò),可以把它分享出去讓更多的人看到吧!
免責(zé)聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點(diǎn)不代表本網(wǎng)站立場(chǎng),如果涉及侵權(quán)請(qǐng)聯(lián)系站長(zhǎng)郵箱:is@yisu.com進(jìn)行舉報(bào),并提供相關(guān)證據(jù),一經(jīng)查實(shí),將立刻刪除涉嫌侵權(quán)內(nèi)容。