使用python如何实现DFT

这篇文章给大家分享的是有关使用python如何实现DFT的内容。小编觉得挺实用的，因此分享给大家做个参考，一起跟随小编过来看看吧。

DFT

DFT(Discrete Fourier Transform)，离散傅里叶变化，可以将离散信号变换到频域，它的公式非常简单:

离散频率下标为k时的频率大小

离散时域信号序列

信号序列的长度，也就是采样的个数

如果你刚接触DFT，并且之前没有信号处理的相关经验，那么第一次看到这个公式，你可能有一些疑惑，为什么这个公式就能进行时域与频域之间的转换呢？

这里，我不打算去解释它，因为我水平有限，说的不清楚。相反，在这里我想介绍，作为一个程序员，如何如实现DFT

从矩阵的角度看DFT

DFT的公式，虽然简单，但是理解起来比较麻烦，我发现如果用矩阵相乘的角度来理解上面的公式，就会非常简单，直接上矩阵：

OK，通过上面的表示，我们很容易将DFT理解成为一种矩阵相乘的操作，这对于我们编码是很容易的。

Talk is cheap, show me the code

根据上面的理解，我们只需要构建出S SS矩阵，然后做矩阵相乘，就等得到DFT的结果

在这之前，我们先介绍如何生成正弦信号，以及如何用scipy中的fft模块进行DFT操作，以验证我们的结果是否正确

正弦信号

A: 幅度

f: 信号频率

n: 时间下标

T: 采样间隔, 等于 1/fs，fs为采样频率

ϕ \phiϕ: 相位

下面介绍如何生成正弦信号

importnumpyasnp
importmatplotlib.pyplotasplt

%matplotlibinline

defgenerate_sinusoid(N,A,f0,fs,phi):
'''
N(int):numberofsamples
A(float):amplitude
f0(float):frequencyinHz
fs(float):samplerate
phi(float):initialphase

return
x(numpyarray):sinusoidsignalwhichlenghtisM
'''

T=1/fs
n=np.arange(N)#[0,1,...,N-1]
x=A*np.cos(2*f0*np.pi*n*T+phi)

returnx

N=511
A=0.8
f0=440
fs=44100
phi=0

x=generate_sinusoid(N,A,f0,fs,phi)

plt.plot(x)
plt.show()

#另一种生成正弦信号的方法，生成时长为t的序列
defgenerate_sinusoid_2(t,A,f0,fs,phi):
'''
t(float):生成序列的时长
A(float):amplitude
f0(float):frequency
fs(float):samplerate
phi(float):initialphase

returns
x(numpyarray):sinusoidsignalsequence
'''

T=1.0/fs
N=t/T

returngenerate_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模块计算DFT

注意，理论上输入信号的长度必须是才能做FFT，而scipy中FFT却没有这样的限制

这是因为当长度不等于时，scipy fft默认做DFT

fromscipy.fftpackimportfft

#generatesinusoid
N=511
A=0.8
f0=440
fs=44100
phi=1.0
x=generate_sinusoid(N,A,f0,fs,phi)

#fftis
X=fft(x)
mX=np.abs(X)#magnitude
pX=np.angle(X)#phase

#plotthemagnitudeandphase
plt.subplot(2,1,1)
plt.plot(mX)

plt.subplot(2,1,2)
plt.plot(pX)
plt.show()

自己实现DFT

自己实现DFT的关键就是构造出S，有两种方式：

defgenerate_complex_sinusoid(k,N):
'''
k(int):frequencyindex
N(int):lengthofcomplexsinusoidinsamples

returns
c_sin(numpyarray):thegeneratedcomplexsinusoid(lengthN)
'''

n=np.arange(N)

c_sin=np.exp(1j*2*np.pi*k*n/N)

returnnp.conjugate(c_sin)

defgenerate_complex_sinusoid_matrix(N):
'''
N(int):lengthofcomplexsinusoidinsamples

returns
c_sin_matrix(numpyarray):thegeneratedcomplexsinusoid(lengthN)
'''

n=np.arange(N)
n=np.expand_dims(n,axis=1)#扩充维度，将1D向量，转为2D矩阵，方便后面的矩阵相乘

k=n

m=n.T*k/N#[N,1]*[1,N]=[N,N]

S=np.exp(1j*2*np.pi*m)#计算矩阵S

returnnp.conjugate(S)

#生成信号，用于测试
N=511
A=0.8
f0=440
fs=44100
phi=1.0
x=generate_sinusoid(N,A,f0,fs,phi)

#第一种方式计算DFT
X_1=np.array([])
forkinrange(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)

#plotthemagnitudeandphase
plt.subplot(2,1,1)
plt.plot(mX)

plt.subplot(2,1,2)
plt.plot(pX)
plt.show()

#结果和scipy的结果基本相同

#第二种方法计算DFT
S=generate_complex_sinusoid_matrix(N)
X_2=np.dot(S,x)

mX=np.abs(X_2)
pX=np.angle(X_2)

#plotthemagnitudeandphase
plt.subplot(2,1,1)
plt.plot(mX)

plt.subplot(2,1,2)
plt.plot(pX)
plt.show()

感谢各位的阅读！关于“使用python如何实现DFT”这篇文章就分享到这里了，希望以上内容可以对大家有一定的帮助，让大家可以学到更多知识，如果觉得文章不错，可以把它分享出去让更多的人看到吧！