信号处理基础：信号的时域和频域分析_（10）.离散傅里叶变换
离散傅里叶变换
引言
在信号处理领域中，离散傅里叶变换（DFT）扮演着非常关键的角色。它实现了时域信号向频域的转换，在众多应用场景中发挥着不可替代的作用。其应用范围极为广泛，在图像处理、音频处理、通信系统以及控制系统等领域均有显著的应用。本节将深入探讨离散傅里叶变换的基本原理及其数学理论基础，并展示如何在Python中实现这一方法。
离散傅里叶变换的数学基础
定义
离散傅里叶变换（DFT）是一种将有限时长的离散时间信号转换为有限频长的离散频率成分分析方法。该技术适用于分析周期性数字信号的频谱特性，在信息处理和通信工程领域具有重要应用价值。对于长度为 N 的序列 x[n] 来说，其数学表达式即为：
X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} k n}
其中，k 是频域的索引，j 是虚数单位，e^{-j \frac{2\pi}{N} k n} 是复指数函数。
逆离散傅里叶变换
逆离散傅里叶变换（Inverse Discrete Fourier Transform, IDFT）将频域信号还原至时域信号。对于一个长度为 N 的频域信号 X[k]，其 IDFT 定义为：
x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} k n}
性质
- 线性和：\text{DFT}是一种线性的数学运算。
- 周期特性和：\text{DFT}的结果展现出明显的周期特性（周期为N）。
- 对称特性和：对于实数值信号x[n]而言，在频域中其\text{DFT}结果满足X[k] = X^*[N-k]这一性质（即共轭对称）。
- 频移影响特性：当信号在时间轴上发生平移时，在频率轴上会引起相位偏移现象。
- 时移影响特性：如果在频率轴上对信号进行平移，则会导致其时间响应发生变化（即产生卷积效应））。
离散傅里叶变换的物理意义
DFT的本质是通过线性组合的方式将时域信号分解为多个正弦波叠加的结果。每个频域样本X[k]则对应于特定频率\frac{k}{N}f_s处的振幅与相位信息，并包含了这些关键参数的数据。
采样定理
该定理表明，在避免失真的情况下实现对连续时间信号的准确还原时所必需遵循的基本原则是将采样频率设定为其最高频率的最小两倍即f_s \geq 2 f_{\text{max}}这一理论对于DFT中的分析至关重要因为它决定了频域中能够区分各组成部分的能力以及最终重建信号的质量
离散傅里叶变换的实现
Python 实现
在 Python 编程语言中可以通过 numpy 和 scipy 库来完成离散傅里叶变换（DFT）和逆离散傅里叶变换（IDFT）。例如 下面展示了一个具体的示例 说明如何利用这些库进行 DFT 和 IDFT 的计算过程。
import numpy as np
import matplotlib.pyplot as plt
# 生成一个时域信号
fs = 1000 # 采样频率
t = np.linspace(0, 1, fs, endpoint=False) # 时间向量
f1 = 5 # 5 Hz 的正弦波
f2 = 15 # 15 Hz 的正弦波
x = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
# 计算 DFT
X = np.fft.fft(x)
# 计算频域的频率
frequencies = np.fft.fftfreq(fs, 1/fs)
# 绘制时域信号
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(t, x)
plt.title('时域信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
# 绘制频域信号
plt.subplot(1, 2, 2)
plt.plot(frequencies, np.abs(X))
plt.title('频域信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 20) # 限制频率范围为 0 到 20 Hz
plt.grid(True)
plt.show()代码解释
生成时域信号 ：
-
表示采样频率fs被设定为1000Hz。
-
变量t代表时间序列，在该区间内生成一系列离散的时间点。
-
参数f_1和f_2分别对应于频率值为5赫兹和15赫兹的正弦波的基本参数。
-
变量x表示由这两个不同频率的正弦波按照一定比例组成的线性组合。
计算 DFT ：
* `np.fft.fft(x)` 计算信号 `x` 的 DFT，结果存储在 `X` 中。计算频域的频率 ：
* `np.fft.fftfreq(fs, 1/fs)` 计算频域的频率点，结果存储在 `frequencies` 中。绘制时域和频域信号 ：
- 使用
matplotlib生成时域信号和频率域信号的大小。 plt.xlim(0, 20)设置频谱图的频率区间在零到二十赫兹之间，从而使五赫兹和十五赫兹的组成部分更加容易辨识。
离散傅里叶变换的应用
信号滤波
DFT 在信号处理领域发挥着重要作用。它能够通过频域中的特定频率成分调整能力来生成多种类型的滤波器；以下是一个经典的低通滤波器案例：
import numpy as np
import matplotlib.pyplot as plt
# 生成一个时域信号
fs = 1000 # 采样频率
t = np.linspace(0, 1, fs, endpoint=False) # 时间向量
f1 = 5 # 5 Hz 的正弦波
f2 = 15 # 15 Hz 的正弦波
x = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
# 计算 DFT
X = np.fft.fft(x)
# 计算频域的频率
frequencies = np.fft.fftfreq(fs, 1/fs)
# 设计低通滤波器
cutoff = 10 # 截止频率
X_filtered = X.copy()
X_filtered[np.abs(frequencies) > cutoff] = 0 # 高于截止频率的成分设为 0
# 计算 IDFT
x_filtered = np.fft.ifft(X_filtered)
# 绘制结果
plt.figure(figsize=(12, 6))
# 时域信号
plt.subplot(2, 2, 1)
plt.plot(t, x)
plt.title('原始时域信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
# 频域信号
plt.subplot(2, 2, 2)
plt.plot(frequencies, np.abs(X))
plt.title('原始频域信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 20)
plt.grid(True)
# 滤波后的频域信号
plt.subplot(2, 2, 3)
plt.plot(frequencies, np.abs(X_filtered))
plt.title('滤波后的频域信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 20)
plt.grid(True)
# 滤波后的时域信号
plt.subplot(2, 2, 4)
plt.plot(t, np.real(x_filtered))
plt.title('滤波后的时域信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
plt.tight_layout()
plt.show()代码解释
生成时域信号 ：
* 与之前的例子相同，生成一个包含 5 Hz 和 15 Hz 正弦波的信号。计算 DFT ：
* 使用 `np.fft.fft(x)` 计算信号 `x` 的 DFT。设计低通滤波器 ：
-
将截止频率设定为 10 Hz。
-
通过计算
np.abs(frequencies)是否大于截止频率cutoff来创建一个布尔掩码。 -
并将其对应的 DFT 值置零。
计算 IDFT ：
通过 np.fft.ifft(X_filtered) 求取滤波后的频域信号的逆离傅里叶变换，并将计算所得的结果赋值给变量 x_filtered。
绘制结果 ：
依次绘制原始时间域信号、原始频率域信号、经过滤波处理后的频率域信号以及经过滤波处理后的时域信号
信号压缩
DFT 也可应用于信号压缩。通过由频域中保留主要频率成分并舍弃次要频率成分的方法, 可以实现数据体积缩减. 例如, 下面是一个简化的例子:
import numpy as np
import matplotlib.pyplot as plt
# 生成一个时域信号
fs = 1000 # 采样频率
t = np.linspace(0, 1, fs, endpoint=False) # 时间向量
f1 = 5 # 5 Hz 的正弦波
f2 = 15 # 15 Hz 的正弦波
x = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
# 计算 DFT
X = np.fft.fft(x)
# 计算频域的频率
frequencies = np.fft.fftfreq(fs, 1/fs)
# 保留前 10% 的频率成分
compression_ratio = 0.1
X_compressed = X.copy()
X_compressed[int(compression_ratio * len(X)):] = 0
X_compressed[:int(compression_ratio * len(X))] = 0
# 计算 IDFT
x_compressed = np.fft.ifft(X_compressed)
# 绘制结果
plt.figure(figsize=(12, 6))
# 时域信号
plt.subplot(2, 2, 1)
plt.plot(t, x)
plt.title('原始时域信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
# 频域信号
plt.subplot(2, 2, 2)
plt.plot(frequencies, np.abs(X))
plt.title('原始频域信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 20)
plt.grid(True)
# 压缩后的频域信号
plt.subplot(2, 2, 3)
plt.plot(frequencies, np.abs(X_compressed))
plt.title('压缩后的频域信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 20)
plt.grid(True)
# 压缩后的时域信号
plt.subplot(2, 2, 4)
plt.plot(t, np.real(x_compressed))
plt.title('压缩后的时域信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
plt.tight_layout()
plt.show()代码解释
生成时域信号 ：
* 与之前的例子相同，生成一个包含 5 Hz 和 15 Hz 正弦波的信号。计算 DFT ：
* 使用 `np.fft.fft(x)` 计算信号 `x` 的 DFT。压缩信号 ：
-
将压缩比例设定为10%。
-
通过计算公式
int(compression_ratio * len(X))来确定需要保留的频率成分索引。 -
将不属于上述计算结果范围内的频率成分设值为零。
计算 IDFT ：
注
绘制结果 ：
依次呈现时间域信号、频率域信号及其压缩后的变化过程
离散傅里叶变换的快速算法
快速傅里叶变换（FFT）
该算法被称为 Fast Fourier Transform（FFT），是高效执行 Discrete Fourier Transform (DFT) 的方法。它的计算复杂度达到了 O(N \log N) 级别，在直接计算 DFT 时的计算复杂度则高达 O(N^2)。本方法通过 divide-and-conquer 策略将大规模问题划分为更小的问题来显著提升运算效率。
Python 实现
Python 提供了 numpy 和 scipy 库中 fft 函数的可调用接口以执行 FFT 任务。以下是演示如何利用这些库完成 FFT 转换过程的一个简单示例。
import numpy as np
import matplotlib.pyplot as plt
# 生成一个时域信号
fs = 1000 # 采样频率
t = np.linspace(0, 1, fs, endpoint=False) # 时间向量
f1 = 5 # 5 Hz 的正弦波
f2 = 15 # 15 Hz 的正弦波
x = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
# 计算 FFT
X = np.fft.fft(x)
# 计算频域的频率
frequencies = np.fft.fftfreq(fs, 1/fs)
# 绘制时域和频域信号
plt.figure(figsize=(12, 6))
# 时域信号
plt.subplot(1, 2, 1)
plt.plot(t, x)
plt.title('时域信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
# 频域信号
plt.subplot(1, 2, 2)
plt.plot(frequencies, np.abs(X))
plt.title('频域信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 20)
plt.grid(True)
plt.show()代码解释
生成时域信号 ：
* 与之前的例子相同，生成一个包含 5 Hz 和 15 Hz 正弦波的信号。计算 FFT ：
* 使用 `np.fft.fft(x)` 计算信号 `x` 的 FFT。计算频域的频率 ：
* 使用 `np.fft.fftfreq(fs, 1/fs)` 计算频域的频率点。绘制时域和频域信号 ：
* 分别绘制时域信号和频域信号的幅度。离散傅里叶变换的局限性和改进
局限性
- 频率分辨能力：离散傅里叶变换（DFT）的频谱分辨能力主要由采样点数量N决定，其频率分辨度\Delta f可表示为\Delta f = \frac{f_s}{N}。当N较小时，在频谱分析中将导致较低的频率区分能力。
- 泄漏现象：当信号周期与所选时间窗长度不一致时，在频域中将产生能量泄露现象（即所谓的"泄漏效应"），从而影响频谱准确性。
- 运算复杂度：采用直接计算法进行DFT运算时其时间复杂度为O(N^2)这一特性使得该方法在处理大规模数据集时存在较高的运算效率问题。
改进方法
- 提升采样密度 ：提升N值能够显著地改善频率分辨率。
- 选择合适的窗函数 ：合理选择窗函数有助于有效降低频谱泄露现象；目前常用的主要有矩形窗（rectangular window）、汉宁窗（Hanning window）、哈明窗（Hamming window）等。
- 采用多层FFT算法 ：采用多层快速傅里叶变换算法可明显地提升计算速度和资源利用率。
窗口函数的例子
以下是一个使用汉宁窗函数减少泄漏效应的例子：
import numpy as np
import matplotlib.pyplot as plt
# 生成一个时域信号
fs = 1000 # 采样频率
t = np.linspace(0, 1, fs, endpoint=False) # 时间向量
f1 = 5 # 5 Hz 的正弦波
f2 = 15 # 15 Hz 的正弦波
x = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
# 应用汉宁窗函数
window = np.hanning(fs)
x_windowed = x * window
# 计算 DFT
X_windowed = np.fft.fft(x_windowed)
# 计算频域的频率
frequencies = np.fft.fftfreq(fs, 1/fs)
# 绘制时域和频域信号
plt.figure(figsize=(12, 6))
# 时域信号
plt.subplot(1, 2, 1)
plt.plot(t, x, label='原始信号')
plt.plot(t, x_windowed, label='应用汉宁窗函数后的信号')
plt.title('时域信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
plt.legend()
# 频域信号
plt.subplot(1, 2, 2)
plt.plot(frequencies, np.abs(np.fft.fft(x)), label='无窗口')
plt.plot(frequencies, np.abs(X_windowed), label='应用汉宁窗函数')
plt.title('频域信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 20)
plt.grid(True)
plt.legend()
plt.show()代码解释
生成时域信号 ：
-
如同前例的做法一样，在前面的例子中已经展示了如何创建一个由不同频率正弦波组成的信号。这里同样采用这种方法来生成一个包含5Hz和15Hz正弦波的新信号。
-
在代码实现中，
- 首先设定变量fs表示采样频率并将其赋值为1000Hz；
- 然后创建时间变量t使其覆盖从0到1秒的时间点；
- 接着分别指定f_1和f_2分别为5Hz和15Hz的正弦波的基本参数；
- 最后将这两个不同频率的正弦波进行叠加形成新的信号变量x。
应用汉宁窗函数 ：
*通过调用 np.hanning(fs) 来生成汉宁窗函数 。 *对时域信号 x 进行乘法运算与汉宁窗函数以获得结果 x_windowed *。
计算 DFT ：
* 使用 `np.fft.fft(x_windowed)` 计算应用汉宁窗函数后的信号的 DFT，结果存储在 `X_windowed` 中。计算频域的频率 ：
* 使用 `np.fft.fftfreq(fs, 1/fs)` 计算频域的频率点。绘制时域和频域信号 ：
- 在时域图上描绘原始信号及其施加汉宁窗处理后的结果。
- 在频域图上展示原始信号及其经过汉宁窗处理后对应的DFT变换结果。
- 使用
plt.xlim(0, 20)设置频域图的频率范围至0至20Hz，并以此聚焦于5Hz和15Hz分量的分析。 - 通过
plt.legend()指令添加图例标识符，在图形中明确区分不同来源的DFT变换结果。
离散傅里叶变换的高级应用
频谱分析
频谱分析是DFT的一个关键性应用，在对信号进行频域表示时能够识别出其频率成分、幅度和相位信息；这不仅有助于更好地理解信号特征并对其进一步处理具有重要意义
频谱分析的例子
以下是一个使用 DFT 进行频谱分析的例子：
import numpy as np
import matplotlib.pyplot as plt
# 生成一个时域信号
fs = 1000 # 采样频率
t = np.linspace(0, 1, fs, endpoint=False) # 时间向量
f1 = 5 # 5 Hz 的正弦波
f2 = 15 # 15 Hz 的正弦波
x = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
# 计算 DFT
X = np.fft.fft(x)
# 计算频域的频率
frequencies = np.fft.fftfreq(fs, 1/fs)
# 计算幅度谱和相位谱
magnitude_spectrum = np.abs(X)
phase_spectrum = np.angle(X)
# 绘制频谱
plt.figure(figsize=(12, 6))
# 幅度谱
plt.subplot(1, 2, 1)
plt.plot(frequencies, magnitude_spectrum)
plt.title('幅度谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 20)
plt.grid(True)
# 相位谱
plt.subplot(1, 2, 2)
plt.plot(frequencies, phase_spectrum)
plt.title('相位谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('相位 (弧度)')
plt.xlim(0, 20)
plt.grid(True)
plt.show()代码解释
生成时域信号 ：
* 与之前的例子相同，生成一个包含 5 Hz 和 15 Hz 正弦波的信号。计算 DFT ：
* 使用 `np.fft.fft(x)` 计算信号 `x` 的 DFT，结果存储在 `X` 中。计算频域的频率 ：
* 使用 `np.fft.fftfreq(fs, 1/fs)` 计算频域的频率点。计算幅度谱和相位谱 ：
* 使用 `np.abs(X)` 计算幅度谱。
* 使用 `np.angle(X)` 计算相位谱。绘制频谱 ：
- 分别绘制出幅度谱与相位谱。
plt.xlim(0,20)设置频域图的频率范围为0到20Hz，并有助于更好地识别其5Hz和15Hz的部分。
信号重建
利用IDFT技术, 能够将频域信号转换回时域信号, 并实现信号重建. 该方法不仅有助于验证DFT的结果, 还能支持有效的信号处理工作.
信号重建的例子
以下是一个使用 IDFT 进行信号重建的例子：
import numpy as np
import matplotlib.pyplot as plt
# 生成一个时域信号
fs = 1000 # 采样频率
t = np.linspace(0, 1, fs, endpoint=False) # 时间向量
f1 = 5 # 5 Hz 的正弦波
f2 = 15 # 15 Hz 的正弦波
x = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
# 计算 DFT
X = np.fft.fft(x)
# 计算 IDFT
x_reconstructed = np.fft.ifft(X)
# 绘制结果
plt.figure(figsize=(12, 6))
# 原始时域信号
plt.subplot(1, 2, 1)
plt.plot(t, x, label='原始信号')
plt.plot(t, np.real(x_reconstructed), label='重建信号', linestyle='--')
plt.title('时域信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
plt.legend()
# 频域信号
plt.subplot(1, 2, 2)
plt.plot(frequencies, np.abs(X), label='频域信号')
plt.title('频域信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 20)
plt.grid(True)
plt.legend()
plt.show()代码解释
生成时域信号 ：
* 与之前的例子相同，生成一个包含 5 Hz 和 15 Hz 正弦波的信号。计算 DFT ：
* 使用 `np.fft.fft(x)` 计算信号 `x` 的 DFT，结果存储在 `X` 中。计算 IDFT ：
* 使用 `np.fft.ifft(X)` 计算 DFT 的逆变换，结果存储在 `x_reconstructed` 中。绘制结果 ：
- 生成时域图中的原始信号与重建信号曲线。
- 通过DFT算法在频域图上呈现分析结果。
- 使用Python的matplotlib库中的legend函数，在图像中添加图例标识以区分各条曲线。
总结
在信号处理领域中，在离散傅里叶变换（DFT）的基础上进行研究与应用是一项十分关键的工作；它将时域信号转换为其频域表示；从而使得对信号的分析和处理变得更加高效；借助快速傅里叶变换（FFT），我们可以高效地计算出DFT的结果；这一方法特别适合用于处理大规模的数据集；应用适当的窗函数能够有效减少频谱泄漏现象；从而进一步提升了频谱分析的准确性；在信号滤波、数据压缩以及重构造等方面有着广泛的应用；作为该领域的重要基石之一；理解其原理对于深入研究 signal processing 技术具有重要意义