使用 scipy.fft 進行Fourier Transform：Python 信號處理

目錄

scipy.fft 模塊

安裝 SciPy 和 Matplotlib

scipy.fft 與 scipy.fftpack

scipy.fft 與 numpy.fft

The Fourier Transform

為什么需要Fourier Transform?

Time Domain vs Frequency Domain

Types of Fourier Transforms

實際示例：從音頻中去除不需要的噪音

創建信號

混合音頻信號

使用快速Fourier Transform (FFT)

使用 rfft() 讓它更快

過濾信號

應用逆 FFT

Avoiding Filtering Pitfalls

The Discrete Cosine and Sine Transforms

結論

所述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 進行優化。

安裝 SciPy 和 Matplotlib

在開始之前，您需要安裝 SciPy 和Matplotlib。您可以通過以下兩種方式之一執行此操作：

使用 Anaconda 安裝：下載并安裝Anaconda Individual Edition。它帶有 SciPy 和 Matplotlib，因此一旦您按照安裝程序中的步驟操作，您就大功告成了！

安裝方式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 文檔時，您可能會遇到兩個看起來非常相似的模塊：

scipy.fft

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 實現。

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，可以定義兩個常量：

SAMPLE_RATE確定信號每秒使用多少個數據點來表示正弦波。因此，如果信號的采樣率為 10 Hz，并且是 5 秒的正弦波，那么它就會有10 * 5 = 50數據點。

DURATION?是生成樣本的長度。

接下來，您定義一個函數來生成正弦波，因為您將在以后多次使用它。該函數采用頻率 ，freq然后返回用于繪制波形的x和y值。

正弦波的 x 坐標在0和之間均勻分布DURATION，因此代碼使用 NumPylinspace()來生成它們。它需要一個起始值、一個結束值和要生成的樣本數。設置endpoint=False對于傅立葉變換正常工作很重要，因為它假設信號是周期性的。

np.sin()計算每個 x 坐標處的正弦函數值。結果乘以頻率，使正弦波在該頻率振蕩，乘積乘以 2π 以將輸入值轉換為弧度。

注意：如果您之前沒有學過太多三角學，或者您需要復習，那么請查看可汗學院的三角學課程。

定義函數后，您可以使用它生成一個持續 5 秒的 2 赫茲正弦波，并使用 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.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)

代碼調用了兩個非常重要的函數：

fft()?計算變換本身。

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()進來的地方。它計算復數的 √(a2 + b2)，這是兩個數字的整體大小，重要的是單個值。

注意：順便說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模塊將其應用于過濾問題。

在本教程中，您學習了：

如何以及何時使用?Fourier transform

如何scipy.fft為您的用例選擇正確的函數

時域和頻域有什么區別

如何查看和修改頻譜的信號的

如何使用rfft()可以加速你的代碼

在上一節中，您還了解了離散余弦變換和離散正弦變換。你看到了調用哪些函數來使用它們，并且你學會了何時使用一個而不是另一個。

如果您想要本教程的摘要，則可以下載下面的備忘單。它解釋了scipy.fft模塊中的所有功能以及可用的不同類型轉換的細分：

Python 機器學習

版權聲明：本文內容由網絡用戶投稿，版權歸原作者所有，本站不擁有其著作權，亦不承擔相應法律責任。如果您發現本站中有涉嫌抄襲或描述失實的內容，請聯系我們jiasou666@gmail.com 處理，核實后本網站將在24小時內刪除侵權內容。

版權聲明：本文內容由網絡用戶投稿，版權歸原作者所有，本站不擁有其著作權，亦不承擔相應法律責任。如果您發現本站中有涉嫌抄襲或描述失實的內容，請聯系我們jiasou666@gmail.com 處理，核實后本網站將在24小時內刪除侵權內容。