Python科学计算——复杂信号FFT

Posted 2023-03-19

参考技术A FFT (Fast Fourier Transform, 快速傅里叶变换) 是离散傅里叶变换的快速算法，也是数字信号处理技术中经常会提到的一个概念。用快速傅里叶变换能将时域的数字信号转换为频域信号，转换为频域信号后我们可以很方便地分析出信号的频率成分。

当我们把双频信号FFT示例中的 fft_size 的值改为 2**12 时，这时，基频为 16Hz，不能被 1kHz整除，所以 1kHz 处发生了频谱泄露，而它能被 4kHz 整除，所以 4kHz 可以很好地被采样。

由于波形的前后不是连续的，出现波形跳变，而跳变处有着非常广泛的频谱，因此FFT的结果中出现了频谱泄漏。

为了减小FFT所截取的数据段前后的跳变，可以对数据先乘以一个窗函数，使得其前后数据能平滑过渡。常用的hanning窗函数的定义如下：

50Hz 正弦波与hann窗函数乘积之后的重复波形如下：

我们对频谱泄漏示例中的1kHz 和 4kHz 信号进行了 hann 窗函数处理，可以看出能量更加集中在 1kHz 和 4kHz，在一定程度上抑制了频谱泄漏。

以 1kHz 三角波为例，我们知道三角波信号中含有丰富的频率信息，它的傅里叶级数展开为：

当数字信号的频率随时间变化时，我们称之为扫频信号。以频率随时间线性变化的扫频信号为例，其数学形式如下：

其频率随时间线性变化，当我们在 [0,1] 的时间窗口对其进行采样时，其频率范围为 0~5kHz。当时间是连续时，扫频信号的频率也是连续的。但是在实际的处理中，是离散的点采样，因此时间是不连续的，这就使扫频信号的快速傅里叶变换问题退化为多点频信号快速傅里叶变换问题。其快速傅里叶变换得到的频谱图如下所示：

以 50Hz 正弦信号相位调制到 1kHz 的信号为例，其信号形式如下：

它的时域波形，频率响应和相位响应如下图所示：

以扫频信号为例，当我们要探究FFT中的能量守恒时，我们要回归到信号最初的形式：

使用 scipy.fft 进行Fourier Transform：Python 信号处理

摘要：Fourier transform 是一个强大的概念，用于各种领域，从纯数学到音频工程甚至金融。

本文分享自华为云社区《使用 scipy.fft 进行Fourier Transform：Python 信号处理》，作者： Yuchuan。

scipy.fft模块

傅立叶变换是许多应用中的重要工具，尤其是在科学计算和数据科学中。因此，SciPy 长期以来一直提供它的实现及其相关转换。最初，SciPy 提供了该scipy.fftpack模块，但后来他们更新了他们的实现并将其移到了scipy.fft模块中。

SciPy 充满了功能。有关该库的更一般介绍，请查看Scientific Python：使用 SciPy 进行优化。

安装 SciPy 和 Matplotlib

在开始之前，您需要安装 SciPy 和Matplotlib。您可以通过以下两种方式之一执行此操作：

1. 使用 Anaconda 安装：下载并安装Anaconda Individual Edition。它带有 SciPy 和 Matplotlib，因此一旦您按照安装程序中的步骤操作，您就大功告成了！
2. 安装方式pip：如果您已经pip安装，那么您可以使用以下命令安装库：

$ python -m pip install -U scipy matplotlib

您可以通过在终端中键入python并运行以下代码来验证安装是否有效：

>>>
>>> import scipy, matplotlib
>>> print(scipy.__file__)
/usr/local/lib/python3.6/dist-packages/scipy/__init__.py
>>> print(matplotlib.__file__)
/usr/local/lib/python3.6/dist-packages/matplotlib/__init__.py

此代码导入 SciPy 和 Matplotlib 并打印模块的位置。您的计算机可能会显示不同的路径，但只要它打印路径，安装就成功了。

SciPy 现已安装！现在是时候看看scipy.fft和之间的区别了scipy.fftpack。

scipy.fft 对比 scipy.fftpack

在查看 SciPy 文档时，您可能会遇到两个看起来非常相似的模块：

1. scipy.fft
2. scipy.fftpack

该scipy.fft模块较新，应该优先于scipy.fftpack. 您可以在SciPy 1.4.0的发行说明中阅读有关更改的更多信息，但这里有一个快速摘要：

• scipy.fft 有一个改进的 API。
• scipy.fft允许使用多个 worker，这可以在某些情况下提供速度提升。
• scipy.fftpack被认为是遗留的，SciPy 建议scipy.fft改用。

除非您有充分的理由使用scipy.fftpack，否则您应该坚持使用scipy.fft.

scipy.fft 对比 numpy.fft

SciPy 的快速傅立叶变换 (FFT)实现包含更多功能，并且比 NumPy 的实现更可能修复错误。如果有选择，您应该使用 SciPy 实现。

NumPy 维护了 FFT 实现以实现向后兼容性，尽管作者认为像傅立叶变换这样的功能最好放在 SciPy 中。有关更多详细信息，请参阅SciPy 常见问题解答。

Fourier Transform

Fourier 分析是研究如何将数学函数分解为一系列更简单的三角函数的领域。傅立叶变换是该领域的一种工具，用于将函数分解为其分量频率。

好吧，这个定义非常密集。就本教程而言，傅立叶变换是一种工具，可让您获取信号并查看其中每个频率的功率。看看这句话中的重要术语：

• 一个信号是随时间变化的信息。例如，音频、视频和电压轨迹都是信号的例子。
• 甲频率是某物重复的速度。例如，时钟以一赫兹 (Hz) 的频率滴答，或每秒重复一次。
• 在这种情况下，功率仅表示每个频率的强度。

下图是一些正弦波的频率和功率的直观演示：

所述的峰值的高频正弦波比那些更靠近在一起低频正弦波，因为它们重复得更频繁。所述低功率正弦波具有比其它两个正弦波较小的峰。

为了更具体地说明这一点，假设您对某人同时在钢琴上弹奏三个音符的录音使用了傅立叶变换。结果频谱将显示三个峰值，每个音符一个。如果该人弹奏一个音符比其他音符更轻柔，那么该音符的频率强度将低于其他两个。

这是钢琴示例在视觉上的样子：

钢琴上最高的音符比其他两个音符更安静，因此该音符的频谱具有较低的峰值。

为什么需要Fourier Transform?

傅立叶变换在许多应用中都很有用。例如，Shazam和其他音乐识别服务使用傅立叶变换来识别歌曲。JPEG 压缩使用傅立叶变换的变体来去除图像的高频分量。语音识别使用傅立叶变换和相关变换从原始音频中恢复口语。

通常，如果您需要查看信号中的频率，则需要进行傅立叶变换。如果在时域中处理信号很困难，那么使用傅立叶变换将其移动到频域中是值得尝试的。在下一节中，您将了解时域和频域之间的差异。

Time Domain vs Frequency Domain

在本教程的其余部分，您将看到术语time domain和frequency domain。这两个术语指的是查看信号的两种不同方式，要么是其分量频率，要么是随时间变化的信息。

在时域中，信号是幅度（y 轴）随时间（x 轴）变化的波。您很可能习惯于在时域中看到图形，例如这个：

这是一些音频的图像，它是一个时域信号。横轴表示时间，纵轴表示幅度。

在频域中，信号表示为一系列频率（x 轴），每个频率都有相关的功率（y 轴）。下图是

上述经过Fourier Transform 后的音频信号：

此处，之前的音频信号由其组成频率表示。底部的每个频率都有一个相关的功率，产生您看到的频谱。

有关频域的更多信息，请查看DeepAI 词汇表条目。

Fourier Transform的类型

傅立叶变换可以细分为不同类型的变换。最基本的细分是基于变换操作的数据类型：连续函数或离散函数。本教程将仅处理离散傅立叶变换 (DFT)。

即使在本教程中，您也会经常看到 DFT 和 FFT 这两个术语互换使用。然而，它们并不完全相同。的快速傅立叶变换（FFT）是用于计算离散傅立叶变换（DFT）的算法，而DFT是变换本身。

您将在scipy.fft库中看到的另一个区别是不同类型的输入之间的区别。fft()接受复数值输入，并rfft()接受实数值输入。跳到使用快速傅立叶变换 (FFT) 部分以了解复数和实数。

另外两个变换与 DFT 密切相关：离散余弦变换 (DCT)和离散正弦变换 (DST)。您将在离散余弦和正弦变换部分中了解这些内容。

实际示例：从音频中去除不需要的噪音

为了帮助您理解傅立叶变换以及您可以用它做什么，您将过滤一些音频。首先，您将创建一个带有高音嗡嗡声的音频信号，然后您将使用傅立叶变换去除嗡嗡声。

创建信号

正弦波有时被称为纯音，因为它们代表单一频率。您将使用正弦波来生成音频，因为它们将在生成的频谱中形成不同的峰值。

正弦波的另一个优点是它们可以使用 NumPy 直接生成。如果您之前没有使用过 NumPy，那么您可以查看什么是 NumPy？

这是一些生成正弦波的代码：

import numpy as np
from matplotlib import pyplot as plt

SAMPLE_RATE = 44100  # Hertz
DURATION = 5  # Seconds

def generate_sine_wave(freq, sample_rate, duration):
    x = np.linspace(0, duration, sample_rate * duration, endpoint=False)
    frequencies = x * freq
    # 2pi because np.sin takes radians
    y = np.sin((2 * np.pi) * frequencies)
    return x, y

# Generate a 2 hertz sine wave that lasts for 5 seconds
x, y = generate_sine_wave(2, SAMPLE_RATE, DURATION)
plt.plot(x, y)
plt.show()

你以后导入与NumPy和Matplotlib，可以定义两个常量：

1. SAMPLE_RATE确定信号每秒使用多少个数据点来表示正弦波。因此，如果信号的采样率为 10 Hz，并且是 5 秒的正弦波，那么它就会有10 * 5 = 50数据点。
2. DURATION 是生成样本的长度。

接下来，您定义一个函数来生成正弦波，因为您将在以后多次使用它。该函数采用频率 ，freq然后返回用于绘制波形的x和y值。

正弦波的 x 坐标在0和之间均匀分布DURATION，因此代码使用 NumPylinspace()来生成它们。它需要一个起始值、一个结束值和要生成的样本数。设置endpoint=False对于傅立叶变换正常工作很重要，因为它假设信号是周期性的。

np.sin()计算每个 x 坐标处的正弦函数值。结果乘以频率，使正弦波在该频率振荡，乘积乘以 2π 以将输入值转换为弧度。

注意：如果您之前没有学过太多三角学，或者您需要复习，那么请查看可汗学院的三角学课程。

定义函数后，您可以使用它生成一个持续 5 秒的 2 赫兹正弦波，并使用 Matplotlib 绘制它。您的正弦波图应如下所示：

x 轴以秒为单位表示时间，由于每一秒的时间有两个峰值，您可以看到正弦波每秒振荡两次。此正弦波的频率太低而无法听到，因此在下一部分中，您将生成一些更高频率的正弦波，您将了解如何混合它们。

混合音频信号

好消息是混合音频信号只包括两个步骤：

1. 将信号相加
2. 标准化结果

在将信号混合在一起之前，您需要生成它们：

_, nice_tone = generate_sine_wave(400, SAMPLE_RATE, DURATION)
_, noise_tone = generate_sine_wave(4000, SAMPLE_RATE, DURATION)
noise_tone = noise_tone * 0.3

mixed_tone = nice_tone + noise_tone

此代码示例中没有任何新内容。它生成分别分配给变量 nice_tone和 的中音和高音noise_tone。您将使用高音作为您不需要的噪音，因此它会乘以0.3降低其功率。然后代码将这些音调加在一起。请注意，您使用下划线 ( _) 来丢弃由x返回的值generate_sine_wave()。

下一步是归一化，或缩放信号以适应目标格式。由于您稍后将如何存储音频，您的目标格式是一个 16 位整数，范围从 -32768 到 32767：

normalized_tone = np.int16((mixed_tone / mixed_tone.max()) * 32767)

plt.plot(normalized_tone[:1000])
plt.show()

在这里，代码进行缩放mixed_tone以使其适合 16 位整数，然后使用 NumPy 的np.int16. 除mixed_tone以其最大值将其缩放到-1和之间1。当这个信号乘以 时32767，它在-32767和之间缩放32767，大致是 的范围np.int16。该代码仅绘制第一个1000样本，以便您可以更清楚地看到信号的结构。

你的情节应该是这样的：

信号看起来像失真的正弦波。您看到的正弦波是您生成的 400 Hz 音调，失真是 4000 Hz 音调。如果仔细观察，您会发现失真呈正弦波形状。

要收听音频，您需要将其存储为音频播放器可以读取的格式。最简单的方法是使用 SciPy 的wavfile.write方法将其存储在 WAV 文件中。16 位整数是 WAV 文件的标准数据类型，因此您需要将信号标准化为 16 位整数：

from scipy.io.wavfile import write

# Remember SAMPLE_RATE = 44100 Hz is our playback rate
write("mysinewave.wav", SAMPLE_RATE, normalized_tone)

此代码将写入mysinewave.wav您运行 Python 脚本的目录中的文件。然后，您可以使用任何音频播放器甚至Python收听此文件。您会听到较低的音调和较高的音调。这些是您混合的 400 Hz 和 4000 Hz 正弦波。

完成此步骤后，您的音频样本就准备好了。下一步是使用傅立叶变换去除高音！

使用快速Fourier Transform (FFT)

是时候在生成的音频上使用 FFT 了。FFT 是一种实现傅立叶变换的算法，可以计算时域中信号的频谱，例如您的音频：

from scipy.fft import fft, fftfreq

# Number of samples in normalized_tone
N = SAMPLE_RATE * DURATION

yf = fft(normalized_tone)
xf = fftfreq(N, 1 / SAMPLE_RATE)

plt.plot(xf, np.abs(yf))
plt.show()

此代码将计算生成的音频的Fourier Transform 并绘制它。在分解它之前，先看看它产生的情节：

您可以看到正频率中的两个峰值和负频率中这些峰值的镜像。正频率峰值位于 400 Hz 和 4000 Hz，这对应于您输入音频的频率。

Fourier transform 已经将您复杂的、微弱的信号转化为它所包含的频率。因为你只输入了两个频率，所以只有两个频率出来了。负正对称性是将实值输入放入 Fourier transform 的副作用，但稍后您会听到更多相关信息。

在前几行中，您导入scipy.fft稍后将使用的函数，并定义一个变量N，用于存储信号中的样本总数。

在这之后是最重要的部分，计算 Fourier transform ：

yf = fft(normalized_tone)
xf = fftfreq(N, 1 / SAMPLE_RATE)

代码调用了两个非常重要的函数：

1. fft() 计算变换本身。
2. fftfreq()计算 的输出中每个bin中心的频率fft()。没有这个，就无法在频谱上绘制 x 轴。

甲箱是已经被分组，就像在一个值的范围的直方图。有关bin 的更多信息，请参阅此信号处理堆栈交换问题。出于本教程的目的，您可以将它们视为单个值。

一旦您获得了 Fourier transform 的结果值及其相应的频率，您就可以绘制它们：

plt.plot(xf, np.abs(yf))
plt.show()

这段代码有趣的部分是您yf在绘制之前所做的处理。你呼吁np.abs()是yf因为它的价值是复杂的。

甲复数是一个数，其具有两个部分，即实部和虚部。定义这样的数字很有用的原因有很多，但您现在需要知道的是它们存在。

数学家通常以a + bi的形式书写复数，其中a是实部，b是虚部。b后面的i表示b是虚数。

注意：有时您会看到使用i编写的复数，有时您会看到使用j编写的复数，例如 2 + 3 i和 2 + 3 j。两者是一样的，但i被数学家用得更多，而j被工程师用得更多。

有关复数的更多信息，请查看可汗学院的课程或数学很有趣页面。

由于复数有两个部分，将它们与二维轴上的频率作图需要您从它们计算一个值。这就是np.abs()进来的地方。它计算复数的 √(a² + b²)，这是两个数字的整体大小，重要的是单个值。

注意：顺便说fft()一句，您可能已经注意到，准确地说，返回的最大频率刚好超过 20,000 赫兹，即 22050Hz。这个值正好是我们采样率的一半，称为奈奎斯特频率。

这是信号处理中的一个基本概念，意味着您的采样率必须至少是信号最高频率的两倍。

让它更快 rfft()

fft()输出的频谱绕y轴反射，因此负半部是正半部的镜子。这种对称性是由向变换输入实数（不是复数）引起的。

您可以使用这种对称性，通过只计算它的一半来使您的 Fourier transform 更快。scipy.fft以rfft().

很棒的一点rfft()是，它是fft(). 记住之前的 FFT 代码：

yf = fft(normalized_tone)
xf = fftfreq(N, 1 / SAMPLE_RATE)

换入rfft()，代码基本保持不变，只是有几个关键的变化：

from scipy.fft import rfft, rfftfreq

# Note the extra 'r' at the front
yf = rfft(normalized_tone)
xf = rfftfreq(N, 1 / SAMPLE_RATE)

plt.plot(xf, np.abs(yf))
plt.show()

由于rfft()只返回一半的输出fft()，它使用不同的函数来获取频率映射，rfftfreq()而不是fftfreq()。

rfft()仍然会产生复杂的输出，因此绘制其结果的代码保持不变。但是，该图应如下所示，因为负频率将消失：

您可以看到上图只是fft()产生的频谱的正侧。rfft()从不计算频谱的负半部分，这使得它比使用fft().

usingrfft()可以比 using 快两倍fft()，但某些输入长度比其他输入长度快。如果你知道你只会使用实数，那么这是一个值得了解的速度技巧。

现在您有了信号的频谱，您可以继续对其进行滤波。

过滤信号

傅立叶变换的一大优点是它是可逆的，因此您在频域中对信号所做的任何更改都将在您将其变换回时域时应用。您将利用这一点来过滤音频并去除高频。

警告：本节中演示的过滤技术不适用于现实世界的信号。有关原因的解释，请参阅避免过滤陷阱部分。

返回的值rfft()代表每个频率仓的功率。如果您将给定 bin 的功率设置为零，则该 bin 中的频率将不再出现在生成的时域信号中。

使用 的长度xf、最大频率以及频率区间均匀分布的事实，您可以计算出目标频率的索引：

# The maximum frequency is half the sample rate
points_per_freq = len(xf) / (SAMPLE_RATE / 2)

# Our target frequency is 4000 Hz
target_idx = int(points_per_freq * 4000)

然后，您可以在目标频率周围的索引处设置yf为0以摆脱它：

yf[target_idx - 1 : target_idx + 2] = 0

plt.plot(xf, np.abs(yf))
plt.show()

您的代码应生成以下图：

既然只有一个峰，看来是奏效了！接下来，您将应用Fourier Transform返回时域。

应用逆 FFT

应用逆 FFT 类似于应用 FFT：

from scipy.fft import irfft

new_sig = irfft(yf)

plt.plot(new_sig[:1000])
plt.show()

由于您正在使用rfft()，您需要使用irfft()来应用逆。但是，如果您使用了fft()，则反函数将是ifft()。你的情节现在应该是这样的：

如您所见，您现在有一个以 400 Hz 振荡的正弦波，并且您已成功消除了 4000 Hz 噪声。

再一次，您需要在将信号写入文件之前对其进行标准化。你可以像上次那样做：

norm_new_sig = np.int16(new_sig * (32767 / new_sig.max()))

write("clean.wav", SAMPLE_RATE, norm_new_sig)

当您收听此文件时，您会听到烦人的噪音消失了！

避免过滤陷阱

上面的示例更多用于教育目的，而不是实际使用。在真实世界的信号（例如一首音乐）上复制该过程可能会引入比消除更多的嗡嗡声。

在现实世界中，您应该使用scipy.signal包中的滤波器设计函数来过滤信号。过滤是一个涉及大量数学的复杂主题。有关详细介绍，请查看科学家和工程师数字信号处理指南。

离散余弦和正弦变换

scipy.fft如果不了解离散余弦变换 (DCT)和离散正弦变换 (DST)，则有关该模块的教程将是不完整的。这两个变换与 Fourier transform 密切相关，但完全对实数进行运算。这意味着他们将一个实值函数作为输入，并产生另一个实值函数作为输出。

SciPy 将这些转换实现为dct()和dst()。的i*和*n变体是逆和Ñ的功能维版本，分别。

DCT 和 DST 有点像共同构成 Fourier transform 的两半。这并不完全正确，因为数学要复杂得多，但它是一个有用的心智模型。

因此，如果 DCT 和 DST 就像 Fourier transform 的一半，那么它们为什么有用？

一方面，它们比完整的 Fourier transform 更快，因为它们有效地完成了一半的工作。它们甚至可以比rfft(). 最重要的是，它们完全以实数工作，因此您永远不必担心复数。

在学习如何在它们之间进行选择之前，您需要了解偶函数和奇函数。偶函数关于 y 轴对称，而奇函数关于原点对称。要直观地想象这一点，请查看以下图表：

您可以看到偶数函数关于 y 轴对称。奇函数关于y = -x对称，这被描述为关于原点对称。

当您计算傅立叶变换时，您假装正在计算它的函数是无限的。完整的傅立叶变换 (DFT) 假设输入函数无限重复。然而，DCT 和 DST 假设函数是通过对称扩展的。DCT 假设函数以偶对称扩展，而 DST 假设它以奇对称扩展。

下图说明了每个变换如何想象函数扩展到无穷大：

在上图中，DFT 按原样重复了该函数。DCT 垂直镜像函数以扩展它，而 DST 则水平镜像它。

请注意，DST 隐含的对称性会导致函数出现大幅跳跃。这些被称为不连续性，并在结果频谱中产生更多的高频分量。因此，除非您知道您的数据具有奇对称性，否则您应该使用 DCT 而不是 DST。

DCT 非常常用。还有更多的例子，但 JPEG、MP3 和 WebM 标准都使用 DCT。

结论

Fourier transform 是一个强大的概念，用于各种领域，从纯数学到音频工程甚至金融。您现在已经熟悉了离散 Fourier transform ，并且可以很好地使用该scipy.fft模块将其应用于过滤问题。

点击关注，第一时间了解华为云新鲜技术~