Python实现离散傅里叶变换(DFT)与正弦信号生成详解
本文详细解析离散傅里叶变换(DFT)的数学原理，并通过矩阵乘法视角简化理解。重点介绍在Python中如何生成正弦信号，利用scipy.fftpack进行FFT验证，并分别通过循环和矩阵运算两种方式手动实现DFT计算，提供完整代码示例与结果对比。
DFT基础概念
DFT
DFT（Discrete Fourier Transform，离散傅里叶变换）能够将离散时间信号转换至频域进行分析。其核心计算公式如下所示：
其中，
表示在离散频率下标为k时的信号幅度大小。
代表离散时域中的信号序列。
表示信号序列的长度，即采样的总点数。
对于初次接触DFT且缺乏信号处理背景的读者而言，直接面对该公式可能会产生疑惑：为何该公式能实现时域与频域之间的转换？
鉴于解释数学推导的复杂性，本文不深入探讨其理论证明，而是聚焦于程序员视角，展示如何在代码中实现DFT。
从矩阵角度理解DFT
从矩阵角度理解DFT
虽然DFT公式形式简洁，但直观理解具有一定难度。通过矩阵乘法的视角来解读该公式，可以极大地简化理解过程。具体矩阵表示如下：
通过上述矩阵表示，我们可以将DFT直观地理解为一种矩阵乘法操作，这为后续的编程实现提供了清晰的思路。
Python代码实现DFT
Talk is cheap, show me the code
基于上述理解，实现DFT的关键在于构建出矩阵S，随后执行矩阵乘法即可得到DFT的计算结果。
在正式实现之前，我们先介绍如何生成正弦信号，并使用scipy库中的fft模块进行DFT计算，以此验证我们手动实现结果的准确性。
正弦信号生成
正弦信号生成
A：表示信号的幅度。
f：表示信号的频率。
n：表示时间下标。
T：表示采样间隔，其值等于1/fs，其中fs为采样频率。
ϕ：表示信号的初始相位。
接下来介绍具体的正弦信号生成方法。
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()
# 另一种生成正弦信号的方法，生成时长为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)
Scipy FFT模块使用
Scipy FFT模块使用
本节介绍如何利用Scipy的FFT模块来计算DFT。
需要注意的是，理论上输入信号的长度必须为
才能执行标准的FFT运算，而scipy库中的FFT函数并没有这一严格限制。
当信号长度不等于时，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()
手动实现DFT
手动实现DFT
手动实现DFT的核心在于构造矩阵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) # 扩充维度，将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的结果基本相同
# 第二种方法计算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()
总结
总结
本文回顾了DFT的计算公式，并尝试从矩阵相乘的角度重新解读DFT的本质。
详细介绍了两种生成正弦信号的方法。
展示了两种不同的DFT计算实现方式。
完整代码示例见下文。
以上关于信号生成及DFT的Python实现方式的全部内容分享完毕，希望能为读者提供参考，并欢迎大家继续支持本站。
本站声明
本文内容仅供参考，不代表本站观点。转载或引用请注明出处。如发现内容问题，请通过联系我们页面反馈。
