信號(hào)生成及DFT的python實(shí)現(xiàn)方式

更新時(shí)間：2020年02月25日 13:02:28 作者：芥末的無(wú)奈

今天小編就為大家分享一篇信號(hào)生成及DFT的python實(shí)現(xiàn)方式，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(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)槲宜接邢?，說(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()