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

目录

所述Fourier Transform是用于分析的有力工具信号和一切从音频处理到图像压缩时使用。SciPy 在其scipy.fft模块中提供了一个成熟的实现，在本教程中，您将学习如何使用它。

该scipy.fft模块一开始可能看起来很吓人，因为有很多功能，通常名称相似，并且文档使用了大量技术术语而没有解释。好消息是您只需要了解一些核心概念即可开始使用该模块。

如果您对数学不满意，请不要担心！您将通过具体示例对算法有所了解，如果您想深入了解方程式，还将提供指向更多资源的链接。有关Fourier Transform如何工作的直观介绍。

在本教程中，您将学习：

• 如何以及何时使用Fourier Transform
• 如何scipy.fft为您的用例选择正确的函数
• 如何查看和修改频谱的信号的
• 哪些不同的转换可用scipy.fft

如果您希望在阅读完后保留本教程的摘要，请下载下面的备忘单。它解释了scipy.fft模块中的所有功能以及可用的不同类型转换的细分：

该scipy.fft模块

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

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

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

Fourier Transform

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

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

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

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

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

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

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

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

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

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

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

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

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样本，以便您可以更清楚地看到信号的结构。

你的情节应该是这样的：

from scipy.io.wavfile import write

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

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

您可以看到正频率中的两个峰值和负频率中这些峰值的镜像。正频率峰值位于 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是虚数。

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

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

让它更快 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()仍然会产生复杂的输出，因此绘制其结果的代码保持不变。但是，该图应如下所示，因为负频率将消失：

# 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[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()。你的情节现在应该是这样的：

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

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

您可以看到偶数函数关于 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模块将其应用于过滤问题。

在本教程中，您学习了：

• 如何以及何时使用 Fourier transform 
• 如何scipy.fft为您的用例选择正确的函数
• 时域和频域有什么区别
• 如何查看和修改频谱的信号的
• 如何使用rfft()可以加速你的代码

在上一节中，您还了解了离散余弦变换和离散正弦变换。你看到了调用哪些函数来使用它们，并且你学会了何时使用一个而不是另一个。

如果您想要本教程的摘要，则可以下载下面的备忘单。它解释了scipy.fft模块中的所有功能以及可用的不同类型转换的细分：