摘要:Fourier transform 是一个弱小的概念,用于各种畛域,从纯数学到音频工程甚至金融。

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 domainfrequency 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模块将其利用于过滤问题。

最初

文章对你有帮忙的话,记得帮作者点点赞
接下来还会继续跟新无关Python的文章,点点关注不迷路