摘要:傅立叶变换是一个强大的概念,被广泛应用于从纯数学到音频工程甚至金融等各个领域。scipy.fft模块傅立叶变换是许多应用程序中的重要工具,尤其是在科学计算和数据科学中。因此,SciPy早就提供了它的实现及其相关的转换。最初,SciPy提供了scipy.fftpack模块,但后来他们更新了它们的实现并将其移动到scipy.fft模块中。SciPy充满了特性。有关该库的更一般性介绍,请查看ScientificPython:使用SciPy进行优化。安装SciPy和Matplotlib在开始之前,您需要安装SciPy和Matplotlib。您可以通过以下两种方式之一执行此操作:使用Anaconda安装:下载并安装AnacondaIndividualEdition。它带有SciPy和Matplotlib,所以一旦您按照安装程序中的步骤进行操作,就大功告成了!通过pip安装:如果安装了pip,则可以使用以下命令安装库:$python-mpipinstall-Uscipymatplotlib您可以通过在终端中键入python并运行来验证安装是否有效:>>>>>>导入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文档时,您可能会遇到两个看起来非常相似的模块:scipy.fftscipy.fftpackscipy.fft模块较新,应该优先于scipy.fftpack。您可以在阅读更多有关SciPy1.4.0发行说明中的??更改的信息中找到它,但这里有一个快速摘要:scipy.fft具有改进的API。scipy.fft允许使用多个worker,这在某些情况下可以提高速度。scipy.fftpack被认为是遗留的,SciPy建议改用scipy.fft。除非你有充分的理由使用scipy.fftpack,否则你应该坚持使用scipy.fft.scipy.fftvs.numpy.fftSciPy的快速傅立叶变换(FFT)实现包含更多功能,并且比NumPy的实现更有可能修复错误.如果有选择,您应该使用SciPy实现。NumPy维护FFT实现以实现向后兼容性,尽管作者认为像傅里叶变换这样的函数最好放在SciPy中。有关详细信息,请参阅SciPy常见问题解答。傅里叶变换傅里叶分析是研究如何将一个数学函数分解为一系列更简单的三角函数的领域。傅立叶变换是该领域中用于将函数分解为其分量频率的工具。好吧,这个定义非常密集。出于本教程的目的,傅里叶变换是一种工具,可让您获取信号并查看其中每个频率的功率。看看这句话中的重要术语:信号是随时间变化的信息。例如,音频、视频和电压轨迹都是信号的例子。频率是重复某事的速率。例如,时钟以一赫兹(Hz)的速率滴答作响,或每秒重复一次。在这种情况下,功率仅表示每个频率的强度。下图是一些正弦波的频率和功率的直观演示:高频正弦波的峰值比低频正弦波的峰值靠得更近,因为它们重复的频率更高。低功率正弦波的峰值小于其他两个正弦波。为了使这个更具体,假设您将傅立叶变换应用于某人在钢琴上同时弹奏三个音符的录音。生成的频谱将显示三个峰值,每个音符一个。如果此人弹奏一个音符比其他音符更轻柔,则该音符的频率强度将低于其他两个音符。这是钢琴示例的视觉效果:钢琴上的最高音符比其他两个音符更安静,因此该音符的频谱具有较低的峰值。为什么需要傅里叶变换?傅立叶变换在许多应用中都很有用。例如,Shazam和其他音乐识别服务使用傅里叶变换来识别歌曲。JPEG压缩使用傅立叶变换的变体来去除图像的高频分量。语音识别使用傅立叶变换和相关变换从原始音频中恢复口语。通常,如果您需要查看信号中的频率,则需要进行傅立叶变换。如果在时域中处理信号很困难,则值得尝试使用傅立叶变换将其移至频域。在下一节中,您将了解时域和频域之间的区别。时域与频域在本教程的其余部分,您将看到时域和频域这两个术语。这两个术语指的是看待信号的两种不同方式,要么是它的分量频率,要么是它随时间变化的信息。在时域中,信号是幅度(y轴)随时间(x轴)变化的波。您可能习惯于在时域中查看图表,例如:这是一些音频的图像,它是时域信号。横轴表示时间,纵轴表示振幅。在频域中,信号表示为一系列频率(x轴),每个频率都具有相关联的功率(y轴)。下图是上述经过傅里叶变换后的音频信号:这里,前面的音频信号用它的分量频率来表示。底部的每个频率都有一个相关的功率,产生您看到的频谱。有关频域的更多信息,请查看DeepAI词汇表条目。傅里叶变换的类型傅里叶变换可以细分为不同类型的变换。最基本的细分是基于转换操作的数据类型:连续或离散。本教程将仅处理离散傅立叶变换(DFT)。即使在本教程中,您也会经常看到术语DFT和FFT互换使用。然而,它们并不相同。快速傅立叶变换(FFT)是一种用于计算离散傅立叶变换(DFT)的算法,它本身就是变换。您将在scipy.fft库中看到的另一个区别是不同类型输入之间的区别。fft()接受复数值输入,而rfft()接受实数值输入。跳至使用快速傅立叶变换(FFT)部分以了解复数和实数。另外两个变换与DFT密切相关:离散余弦变换(DCT)和离散正弦变换(DST)。您将在离散余弦和正弦变换部分了解这些内容。实例:从音频中去除不需要的噪声为了帮助您理解傅立叶变换以及您可以用它做什么,您将过滤一些音频。首先,您将创建一个带有高音嗡嗡声的音频信号,然后您将使用傅立叶变换来消除嗡嗡声。创建信号正弦波有时被称为纯音,因为它们代表单一频率。您将使用正弦波来生成音频,因为它们会在生成的频谱中形成不同的峰值。正弦波的另一个优点是它们可以直接使用NumPy生成。如果您以前没有使用过NumPy,可以查看什么是NumPy?下面是一些生成正弦波的代码:importnumpyasnpfrommatplotlibimportpyplotaspltSAMPLE\_RATE=44100#HertzDURATION=5#Secondsdefgenerate\_sine\_wave(freq,sample\_rate,duration):x=np.linspace(0,duration,sample\_rate\*duration,endpoint=False)frequencys=x\*freq#2pi因为np.sin采用弧度y=np.sin((2\*np.pi)\*frequencies)returnx,y\#生成持续5秒的2赫兹正弦波x,y=generate\_sine\_wave(2,SAMPLE\_RATE,DURATION)plt.plot(x,y)plt.show()youlater使用NumPy和Matplotlib导入,可以定义两个常量:SAMPLE\_RATE确定信号每秒使用多少数据点来表示正弦波。因此,如果信号以10Hz采样并且是5秒的正弦波,它将具有10*5=50个数据点。DURATION是生成样本的长度。接下来,您定义一个函数来生成正弦波,因为您稍后会多次使用它。该函数采用频率freq并返回用于绘制波形的x和y值。正弦波的x坐标在0和DURATION之间均匀分布,因此代码使用NumPylinspace()来生成它们。它需要一个起始值、一个结束值和要生成的样本数。设置endpoint=False对于傅立叶变换的正常工作很重要,因为它假定信号是周期性的。np.sin()计算正弦函数在每个x坐标处的值。结果乘以正弦波振荡的频率,乘以2π将输入值转换为弧度。注意:如果您之前没有学过多少三角学,或者您需要复习,请查看可汗学院的三角学课程。定义函数后,您可以使用它生成持续5秒的2Hz正弦波,并使用Matplotlib绘制它。您的正弦波图应如下所示:x轴表示以秒为单位的时间,并且由于每一秒的时间有两个峰值,您可以看到正弦波每秒振荡两次。此正弦波的频率太低而无法听到,因此在下一节中您将生成一些更高频率的正弦波,并且您将了解如何混合它们。混合音频信号好消息是混合音频信号只包含两个步骤:添加信号以标准化结果在将信号混合在一起之前,您需要生成它们:\_,nice\_tone=generate\_sine\_wave(400,SAMPLE\_RATE,DURATION)\_,noise\_tone=generate\_sine\_wave(4000,SAMPLE\_RATE,DURATION)noise\_tone=noise\_tone\*0.3mixed\_tone=nice\_tone+noise\_tone这个代码示例中没有任何新内容。它生成中音和高音噪声\_tone,分别分配给变量nice\_tone和。您将使用高音作为您不想要的噪音,因此它乘以0.3以降低其功率。然后代码将这些音调加在一起。请注意,您使用下划线(\_)丢弃由generate\_sine\_wave()返回的x值。下一步是标准化或缩放信号以适合目标格式。由于您稍后将如何存储音频,您的目标格式是一个16位整数,范围从-32768到32767:normalized\_tone=np.int16((mixed\_tone/mixed\_tone.max())\*32767)plt.plot(normalized\_tone\[:1000\])plt.show()此处,代码缩放混合\_tone以适应16位整数,然后使用NumPy的np.int16。它的最大值将它缩放到-1到1之间。当这个信号乘以32767时,它缩放到-32767到32767之间,大致是np.int16的范围。该代码仅绘制前1000个样本,因此您可以更清楚地看到信号的结构。你的情节应该是这样的:信号看起来像一个失真的正弦波。您看到的正弦波是您生成的400Hz音调,失真是4000Hz音调。如果仔细观察,您会发现失真呈正弦波状。要收听音频,您需要以音频播放器可以读取的格式存储它。最简单的方法是使用SciPy的wavfile.write方法将其存储在WAV文件中。16位整数是WAV文件的标准数据类型,因此您需要将信号标准化为16位整数:fromscipy.io.wavfileimportwrite\#RememberSAMPLE\_RATE=44100Hzisourplaybackratewrite("mysinewave.wav",SAMPLE\_RATE,normalized\_tone)此代码将写入运行Python脚本的目录中的文件mysinewave.wav。然后,您可以使用任何音频播放器甚至Python收听此文件。您会听到较低的音调和较高的音调。这些是您混合的400Hz和4000Hz正弦波。完成此步骤后,您的音频样本就准备好了。下一步是使用傅里叶变换去除高音!使用快速傅里叶变换(FFT)是时候对生成的音频应用FFT了。FFT是一种实现傅立叶变换以计算时域信号频谱的算法,例如您的音频:fromscipy.fftimportfft,fftfreq\#Numberofsamplesinnormalized\_toneN=SAMPLE\_RATE\*DURATIONyf=fft(normalized\_tone)xf=fftfreq(N,1/SAMPLE\_RATE)plt.plot(xf,np.abs(yf))plt.show()此代码将计算结果音频的傅立叶变换并绘制它。在分解它之前,先看一下它生成的图:您可以看到正频率中的两个峰值和负频率中这些峰值的镜像。正频率峰值位于400Hz和4000Hz,对应于输入音频的频率。傅里叶变换已将您复杂的微弱信号转换为它包含的频率。因为你只输入了两个频率,所以只有两个频率出来了。负-正对称性是将实值输入放入傅里叶变换的副作用,但您稍后会听到更多相关信息。在前几行中,您导入了稍后将使用的scipy.fft函数,并定义了一个变量N来存储信号中的样本总数。之后是最重要的部分,计算傅立叶变换:yf=fft(normalized\_tone)xf=fftfreq(N,1/SAMPLE\_RATE)代码调用两个非常重要的函数:fft()计算变换本身。fftfreq()计算fft()输出中每个bin中心的频率。没有这个,就无法在光谱上绘制x轴。bin是已经分组的直方图,就像一系列值。有关bin的更多信息,请参阅此信号处理堆栈交换问题。出于本教程的目的,您可以将它们视为单个值。一旦你有了傅里叶变换的结果值及其相应的频率,你就可以绘制它们:plt.plot(xf,np.abs(yf))plt.show()这段代码有趣的部分是你绘制yf先前的处理。您对np.abs()的调用是yf因为它的值很复杂。复数是有两部分的数,实部和虚部。出于多种原因,像这样定义数字很有用,但您现在需要知道的是它们存在。数学家通常将复数写成a+bi的形式,其中a是实部,b是虚部。b后面的i表示b是虚数。注意:有时你会看到用i写的复数,有时你会看到用j写的复数,比如2+3i和2+3j。两者是一样的,只是数学家用的多,工程师用的多。有关复数的更多信息,请查看可汗学院课程或数学很有趣页面。由于复数有两个部分,在二维轴上绘制它们与频率的关系需要您从中计算出一个值。这就是np.abs()的用武之地。它计算复数的√(a2+b2),这是两个数字的总体大小,重要的是各个值。注意:顺便说一下fft(),您可能已经注意到返回的最大频率刚刚超过20,000赫兹,准确地说是22050Hz。这个值恰好是我们采样率的一半,称为奈奎斯特频率。这是信号处理中的一个基本概念,意味着您的采样率必须至少是信号最高频率的两倍。为了使rfft()更快,fft()输出的频谱围绕y轴反射,因此负半部分是正半部分的镜像。这种对称性是由输入实数(不是复数)到变换引起的。您可以使用这种对称性通过仅计算它的一半来使您的傅里叶变换更快。scipy.fft以rfft()开头。rfft()的优点在于它是fft()。记住之前的FFT代码:yf=fft(normalized\_tone)xf=fftfreq(N,1/SAMPLE\_RATE)切换到rfft(),代码基本保持不变,但有一些关键变化:来自scipy.fftimportrfft,rfftfreq\#注意前面的额外'r'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()的两倍。使用rfft()的速度可以是使用fft()的两倍,但对于某些输入长度比其他输入长度更快。如果您知道您只会使用实数,那么这是一个值得了解的速度技巧。现在您有了信号的频谱,您可以继续过滤它。滤波后信号的傅立叶变换的一大优点是它是可逆的,因此当您将信号转换回时域时,您在频域中对信号所做的任何更改都会应用。您将使用它来过滤音频并去除高频。警告:本节中演示的过滤技术不适用于真实世界的信号。有关原因的解释,请参阅避免过滤陷阱部分。rfft()返回的值表示每个频率仓的功率。如果将给定bin的功率设置为零,则该bin中的频率将不再出现在生成的时域信号中。利用长度xf,最大频率,以及频率间隔均匀分布的事实,可以计算出目标频率的指标:\#最大频率为采样率点数的一半\_per\_freq=len(xf)/(SAMPLE\_RATE/2)\#Ourtargetfrequencyis4000Hztarget\_idx=int(points\_per\_freq\*4000)然后,您可以在目标频率附近的索引处将yf设置为0以摆脱它:yf\[target\_idx-1:target\_idx+2\]=0plt.plot(xf,np.abs(yf))plt.show()你的代码应该产生下面的图:因为只有一个峰,它似乎正在工作!接下来,您将傅里叶变换应用回时域。应用逆FFT应用逆FFT类似于应用FFT:fromscipy.fftimportirfftnew\_sig=irfft(yf)plt.plot(new\_sig\[:1000\])plt.show()因为你正在使用rfft(),则需要使用irfft()来应用逆函数。但是,如果您使用了fft(),则反函数就是ifft()。您的图现在应该如下所示:如您所见,您现在有一个以400Hz振荡的正弦波,并且您已成功去除4000Hz噪声。同样,您需要在将信号写入文件之前对其进行归一化。你可以像上次那样做:norm\_new\_sig=np.int16(new\_sig\*(32767/new\_sig.max()))write("clean.wav",SAMPLE\_RATE,norm\_new\_sig)当你听这个文件时,你会听到烦人的噪音消失了!避免过滤器陷阱上面的例子更多的是出于教育目的,而不是实际使用。在现实世界的信号(例如一段音乐)上复制该过程可能会引入比消除的更多的嗡嗡声。在现实世界中,您应该使用scipy.signal包中的滤波器设计函数来过滤信号。过滤是一个涉及大量数学的复杂主题。有关详细介绍,请查看TheScientist'sandEngineer'sGuidetoDigitalSignalProcessing。离散余弦和正弦变换scipy.fft如果不了解离散余弦变换(DCT)和离散正弦变换(DST),则本模块的任何教程都是不完整的。这两个变换与傅里叶变换密切相关,但完全对实数进行操作。这意味着他们将一个实值函数作为输入并产生另一个实值函数作为输出。SciPy将这些转换实现为dct()和dst()。i*和*n变体分别是逆函数和n的函数维度版本。DCT和DST有点像傅里叶变换的两半。这不太正确,因为数学要复杂得多,但它是一个有用的心智模型。那么,如果DCT和DST就像傅里叶变换的一半,它们为什么有用呢?一方面,它们比完整的傅里叶变换更快,因为它们有效地完成了一半的工作。它们甚至可以与rfft()进行比较。最重要的是,它们完全适用于实数,因此您永远不必担心复数。在学习如何在它们之间进行选择之前,您需要了解偶函数和奇函数。偶函数关于y轴对称,而奇函数关于原点对称。为了形象化这一点,请看下图:您可以看到偶数函数关于y轴对称。奇函数关于y=-x对称,这被描述为关于原点对称。当您计算傅里叶变换时,您假装计算它的函数是无限的。全傅里叶变换(DFT)假设输入函数无限重复。但是,DCT和DST假设函数是对称扩展的。DCT假设函数以偶对称展开,而DST假设它以奇对称展开。下图说明了每个变换如何想象函数扩展到无穷大:在上图中,DFT按原样重复函数。DCT垂直镜像功能以扩展它,而DST水平镜像它。请注意,DST隐含的对称性可能会导致函数出现较大的跳跃。这些被称为不连续性,并在结果频谱中产生更多的高频分量。因此,除非您知道您的数据具有奇对称性,否则您应该使用DCT而不是DST。DCT非常常用。还有更多例子,但JPEG、MP3和WebM标准都使用DCT。结论傅立叶变换是一个强大的概念,可用于从纯数学到音频工程甚至金融等多个领域。您现在已经熟悉了离散傅立叶变换,并且可以使用scipy.fft模块很好地将其应用于过滤问题。如果上一篇文章对你有帮助,记得给作者点个赞
