使用python如何实现DFT

发布时间:2021-03-23 12:53:16 作者:小新
来源:亿速云 阅读:1554

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

DFT

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

使用python如何实现DFT

使用python如何实现DFT离散频率下标为k时的频率大小

使用python如何实现DFT 离散时域信号序列

使用python如何实现DFT 信号序列的长度,也就是采样的个数

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

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

从矩阵的角度看DFT

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

使用python如何实现DFT

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

Talk is cheap, show me the code

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

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

正弦信号

使用python如何实现DFT

A: 幅度

f: 信号频率

n: 时间下标

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

ϕ \phiϕ: 相位

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

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()

使用python如何实现DFT

# 另一种生成正弦信号的方法,生成时长为t的序列
def generate_sinusoid_2(t, A, f0, fs, phi):
 '''
 t (float) : 生成序列的时长
 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)

使用python如何实现DFT

Scipy FFT

介绍如何Scipy的FFT模块计算DFT

注意,理论上输入信号的长度必须是使用python如何实现DFT才能做FFT,而scipy中FFT却没有这样的限制

这是因为当长度不等于使用python如何实现DFT时,scipy fft默认做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()

使用python如何实现DFT

自己实现DFT

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

使用python如何实现DFT

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)  # 扩充维度,将1D向量,转为2D矩阵,方便后面的矩阵相乘
 
 k = n
 
 m = n.T * k / N     # [N,1] * [1, N] = [N,N]
 
 S = np.exp(1j * 2 * np.pi * m)  # 计算矩阵 S
 
 return np.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([])
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()

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

使用python如何实现DFT

# 第二种方法计算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()

使用python如何实现DFT

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

推荐阅读:
  1. 怎么使用Python实现队列
  2. 使用python怎么实现列表切片

免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。

python dft

上一篇:使用Python实现分数序列求和的案例

下一篇:如何使用Python实现企业微信机器人每天定时发消息示例

相关阅读

您好,登录后才能下订单哦!

密码登录
登录注册
其他方式登录
点击 登录注册 即表示同意《亿速云用户服务条款》