离散傅立叶变换的Python实现
DFT原理、公式、Python代码实现
😀 基本概念
离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是指傅里叶变换在时域和频域上都呈现离散的形式，将时域信号的采样变换为在离散时间傅里叶变换（DTFT）频域的采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号做DFT，也应当对其经过周期延拓成为周期信号再进行变换。在实际应用中，通常采用快速傅里叶变换来高效计算DFT。
🫠 为何需要离散傅立叶变换
傅立叶变换本身具有的三个特点：
- 时间的积分长度是无穷的；
- 频率空间是无穷的；
- 函数f(t)是连续的，其本身也包含了无穷多的点。
正是因为傅立叶变换中这些“无穷”的特点，导致了其不能在计算机上实现，所以就出现了离散傅立叶变换。
现实世界中获得的数据，只能是有限的时间段，且我们只能针对其中有限个点进行采样。那么我们采样得到的数据能让我们对函数原本的形状了解到什么样的程度呢？
这个时候就出现了采样定理：设时间连续信号f(t)，其最高截至频率为f_m，如果用时间间隔为T≤1/(2f_m)对信号f(t)进行抽样，则连续信号f(t)可以被抽样信号唯一的表示。采样定理告诉我们，如果采样频率为f_s，则我们对原函数最高了解到[0, f_s/2]范围内的频率信息。
采样定理，又称香农采样定理，奈奎斯特采样定理，只要采样频率大于或等于有效信号最高频率的两倍，采样值就可以包含原始信号的所有信息，被采样的信号就可以不失真地还原成原始信号。
为了获取连续函数的离散值，我们的抽取时间间隔取T_s。其实在做信号分析前，我们对信号是一无所知的。根据“采样定理”，当我们选取了抽样时间间隔，其实已经确定了原信号能分析的频率范围[0, 1/(T_s*2)](采样频率f_s=1/(T_s))。
采样定理其实我也没搞太明白，但是不影响后面的理解。
😁 DFT的Python实现
Tips / 提示 接下来使用一个Python编程实例，来了解DFT究竟可以干什么。 使用Numpy创建三个不同频率、不同振幅的正弦函数y_0, y_1, y_2，然后将其相加合并成为一个函数y_3。从时域角度看，y_0, y_1, y_2已经完全“融入”y_3中，如果想从y_3中“挑出来”y_0, y_1, y_2中的一个，显然不可能。但是当我们对该函数进行DFT，从频域的角度我们会发现y_0, y_1, y_2函数中的振幅、频率信息仍然保存在函数y_3中。
引包：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft分别创建4个正弦函数：
$$ y_0=3\sin \left( 2\pi \cdot 1\cdot t \right) \\ y_1=\sin \left( 2\pi \cdot 4.3\cdot t \right) \\ y_2=0.5\sin \left( 2\pi \cdot 7\cdot t \right) \\ y_3=y_0+y_1+y_2 $$
# 样本采样率(每秒钟采几个样)
sr = 100
# 样本采样时间间隔
T_s = 1. / sr # T_s=0.01
# 样本采样的时间长度为5s，样本采样的个数为500，所以说最小频率=1/(N*T_s)=1/(500*0.01)=0.2Hz
t = np.arange(0, 5, T_s)
y_0 = 3 * np.sin(2 * np.pi * 1 * t) # ω=2π/T=2πf,sin(ωt)
y_1 = np.sin(2 * np.pi * 4.3 * t)
y_2 = 0.5 * np.sin(2 * np.pi * 7 * t)
y_3 = y_0 + y_1 + y_2绘制出4个函数的图像：
fig, ax = plt.subplots(2, 2, figsize=(6, 6), sharex='col')
fig.subplots_adjust(hspace=0.2, wspace=0.3)
ax[0, 0].plot(t, y_0)
ax[0, 1].plot(t, y_1)
ax[1, 0].plot(t, y_2)
ax[1, 1].plot(t, y_3)
ax[1, 0].set_xlabel('Time')
ax[1, 1].set_xlabel('Time')
ax[0, 0].set_ylabel('Amplitude')
ax[1, 0].set_ylabel('Amplitude')
ax[0, 0].set_title('$y_0=3\sin(2 \pi \cdot t)$')
ax[0, 1].set_title('$y_1=\sin(2 \pi \cdot 4.3t)$')
ax[1, 0].set_title('$y_2=0.5\sin(2 \pi \cdot 7t)$')
ax[1, 1].set_title('$y_3=y_0+y_1+y_2$')
plt.show()
从图中时域角度可以看出，y_0, y_1, y_2共同“融合”组成了y_3。从时域角度，如果想从y_3中分离出y_0, y_1, y_2其中的一个，显然是不可能的。
下面我们对y_3进行傅立叶变换，换一个角度，从频域的角度来看看会有什么不一样的。
# 傅里叶变换结果，返回长度=1/2奈奎斯频率/最小频率=1/2*100/0.2=250，250*2=500
y_3_fft = fft(y_3)
N = len(t)
# 采样数据的idx
n = np.arange(N)
# 总的采样时间
T = N / sr # 5
# 频率区间：奈奎斯频率/2；n/T*sr = n/(N*T_s)=n*w；刚好奈奎斯频率限制和最小刻度值一起给出了频率空间，直接查看freq就懂了
freq = n / T
# 实际幅度
y_3_fft_norm = y_3_fft / N * 2 # 为何要÷N×2？绘图可视化DFT结果：
fig, ax = plt.subplots(1, 2, figsize=(10, 4), dpi=150)
ax[0].stem(freq, np.abs(y_3_fft), 'b', markerfmt=' ', basefmt='-b')
# ax[0].set_xlim(0, 10)
ax[0].set_xlabel('Freq(Hz)')
ax[0].set_ylabel('FFT Amplitude |X(freq)|')
ax[1].stem(freq, np.abs(y_3_fft_norm), 'b', markerfmt=' ', basefmt='-b')
ax[1].set_xlim(0, 10)
ax[1].set_xlabel('Freq(Hz)')
ax[1].set_ylabel('FFT Amplitude |X(freq)|')
ax[1].annotate(
r'$y_0=3\sin(2 \pi \cdot t)$',
xy=(1, 2.8),
xytext=(+30, -30),
textcoords='offset points',
arrowprops=dict(arrowstyle='->', connectionstyle="arc3,rad=.2")
)
ax[1].annotate(
r'$y_1=\sin(2 \pi \cdot 4.3t)$',
xy=(4.2, .4),
xytext=(+10, 50),
textcoords='offset points',
arrowprops=dict(arrowstyle='->', connectionstyle="arc3,rad=-.5")
)
ax[1].annotate(
r'$y_2=0.5\sin(2 \pi \cdot 7t)$',
xy=(7, .4),
xytext=(+10, 30),
textcoords='offset points',
arrowprops=dict(arrowstyle='->', connectionstyle="arc3,rad=-.5")
)
plt.show()
上图中，左边是scipy.fftpack.fft直接返回的结果，可以看出DFT的输出结果是关于采样率的一半对称的（上面我们设置采样率sr = 100，表示每秒钟采100个样）。这一半的采样率被称为奈奎斯特频率（Nyquist frequency）或折叠频率，它是以电子工程师哈里·奈奎斯特的名字命名的。他和克劳德•香农（Claude Shannon）共同定义了Nyquist-Shannon采样定理——如果一个信号只包含低于采样频率一半的频率成分，那么以一定速率采样的信号就可以完全重建，因此，从DFT得到的最高频率输出是采样频率的一半。
还有一个问题是左图中虽然有明显的三个振幅，但是这三个振幅对应的值却与原来函数y_0, y_1, y_2不对应，这是因为离散傅立叶内部公式实现上的原因导致，细节不用纠结，记住这一步就行了。除以N是因为scipy包中封装的离散傅立叶变换公式为了和傅立叶变换公式保持一致，所以内部没有除以N；乘以2是因为由于复数的引入，同一个振幅被分配至两个共轭复数上。
这也就是为什么我们需要将函数返回的振幅值y_3_fft进行y_3_fft_norm = y_3_fft / N * 2后，才可以得到真正的振幅值。
🥰 DFT应用——信号降噪
创建一个原始信号，并为原始信号加上一个随机数噪音：
# 采样频率
sr = 1000
# 采样时间间隔
ts = 1. / sr
# 样本采样点
t = np.arange(0, 2, ts) # 在2s的时间内，每隔(1./sr)进行一次采样
# 原始信号
f_clean = 5
freq = 10
f_clean += 2 * np.sin(2 * np.pi * freq * t + 3) # ω=2π/T=2πf,sin(ωt)=sin(2πf·t)
freq = 30
f_clean += 5 * np.sin(2 * np.pi * freq * t)
# 信号噪音
f_noise = f_clean + 3 * np.random.randn(len(t))又上面的代码可以看出，原始信号由三项组成：常数项5、正弦项2\sin(2\pi \cdot 10 \cdot t+3)、正弦项5\sin(2\pi \cdot 30 \cdot t)以及一个随机噪音项组成。
绘制信号：
# 绘制信号
fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(t, f_noise, linewidth=.5, color='c', label='Noisy')
ax.plot(t, f_clean, linewidth=.5, color='r', label='Clean')
ax.set_xlabel('Sampling Time')
ax.set_ylabel('Amplitude')
ax.legend()
plt.show()
对比原始信号与加上噪音后的信号，可以发现加上噪音后，数据整体大概的趋势不变，但是局部信息受到较大扰动。
进行傅立叶变换：
# 傅里叶变换
"""
1、采样频率sr = 1000
2、采样时间间隔 = 1/sr = 0.001
3、采样点个数 = 2000
4、最小频率 = 1/(N*T_s)=1/(2000/1000) = 0.5Hz
"""
# 傅立叶变换结果，返回长度=1/2奈奎斯频率/最小频率=1/2*1000/0.2=250，再加上负频率，250*2=500
X = fft(f_noise)
N = len(X)
# 采样数据的idx
n = np.arange(N)
# 总的采样时间
T = N / sr
# 频率区间：奈奎斯频率/2；n/T=n/N*sr=n/(N*T_s)=n*w；刚好奈奎斯频率限制和最小刻度值一起给出了频率空间，直接查看freq就懂了
freq = n / T绘图可视化观察DFT结果：
fig, ax = plt.subplots(2, 1, figsize=(10, 8))
ax[0].stem(freq, np.abs(X), 'b', markerfmt=' ', basefmt='-b')
# ax[0].set_xlim(0, 10)
ax[0].set_xlabel('Freq(Hz)')
ax[0].set_ylabel('FFT Amplitude |X(freq)|')
# 实际幅度
X_norm = X / N * 2 # 为何获取真实幅度的时候要✖(2/N)？
ax[1].stem(freq, np.abs(X_norm), 'b', markerfmt=' ', basefmt='-b')
ax[1].set_xlim(-1, 40)
ax[1].set_xlabel('Freq(Hz)')
ax[1].set_ylabel('FFT Amplitude |X(freq)|')
ax[1].annotate(
r'5(Constant terms)',
xy=(0, 5),
xytext=(+10, 50),
textcoords='offset points',
arrowprops=dict(arrowstyle='->', connectionstyle="arc3,rad=-.5")
)
ax[1].annotate(
r'$2\sin(2 \pi \cdot 10 t+3)$',
xy=(10, 2),
xytext=(+10, 30),
textcoords='offset points',
arrowprops=dict(arrowstyle='->', connectionstyle="arc3,rad=-.5")
)
ax[1].annotate(
r'$5\sin(2 \pi \cdot 30 t)$',
xy=(30, 4),
xytext=(+10, 30),
textcoords='offset points',
arrowprops=dict(arrowstyle='->', connectionstyle="arc3,rad=-.5")
)
plt.show()
对DFT结果进行噪音去除，过滤除去低振幅项：
freq_clean = pd.Series(X).apply(lambda x: x if np.abs(x) > 1000 else 0).to_numpy()可视化去除后的效果：
fig, ax = plt.subplots(figsize=(10, 4))
ax.stem(freq, np.abs(freq_clean), 'b', markerfmt=' ', basefmt='-b')
ax.set_xlim(-1, 100)
ax.set_xlabel('Freq(Hz)')
ax.set_ylabel('FFT Amplitude |X(freq)|')
ax.annotate(
r'5(Constant terms)',
xy=(0, 5 * N / 2),
xytext=(+10, 50),
textcoords='offset points',
arrowprops=dict(arrowstyle='->', connectionstyle="arc3,rad=-.5")
)
ax.annotate(
r'$2\sin(2 \pi \cdot 10 t+3)$',
xy=(10, 2 * N / 2),
xytext=(+10, 30),
textcoords='offset points',
arrowprops=dict(arrowstyle='->', connectionstyle="arc3,rad=-.5")
)
ax.annotate(
r'$5\sin(2 \pi \cdot 30 t)$',
xy=(30, 4 * N / 2),
xytext=(+10, 30),
textcoords='offset points',
arrowprops=dict(arrowstyle='->', connectionstyle="arc3,rad=-.5")
)
plt.show()
将过滤后的信号进行傅立叶逆变换：
# 逆傅里叶变换
ix = ifft(freq_clean)可视化观察过滤后的结果：
# 绘制信号
fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(t, ix, linewidth=.5, color='c', label='IFFT')
ax.plot(t, f_clean, linewidth=.5, color='r', linestyle='-.', label='Raw signal', alpha=0.7)
ax.set_xlabel('Sampling Time')
ax.set_ylabel('Amplitude')
ax.legend()
plt.show()
😜 纯代码实现DFT
上面这个公式就是DFT实现的算法，这么一看觉得好简单~😁
import numpy as np
from scipy.fft import fft
def DFT_slow(x):
x = np.asarray(x, dtype=float) # ensure the data type
N = x.shape[0] # get the x array length
n = np.arange(N) # 1d array
k = n.reshape((N, 1)) # 2d array, 10 x 1, aka column array
M = np.exp(-2j * np.pi * k * n / N)
return np.dot(M, x) # [a,b] . [c,d] = ac + bd, it is a sum
x = np.random.random(1024)
np.allclose(DFT_slow(f_noise), fft(f_noise))🙃 傅立叶变换中的哲学
- 你眼中看似落叶纷飞变化无常的世界，实际只是躺在上帝怀中一份早已谱好的乐章。
- The life you live on may have a completely different new face that can only be seen with a kind of transform.
- The big circle is our country, our era. We as an individual is the small tiny inner circle. Without the big circle that drives everything, we are really can’t do too much.
- If you are rich or very successful sometimes, somewhere, it is probably not all because of your own merit, but largely because of your country, people around you, or the good company you work with. Without those big circles that drive you forward, you may not able to achieve what you have now.
腾讯云开发者
