摘要：Fourier transform 是一個強大的概念，用於各種領域，從純數學到音訊工程甚至金融。
本文分享自華為雲社群《使用 scipy.fft 進行Fourier Transform：Python 訊號處理》，作者： Yuchuan。
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 實現。
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，可以定義兩個常量：
- 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()進來的地方。它計算複數的 √(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模組將其應用於過濾問題。