背景
晚上看到一篇文章，介紹了【對20世紀科學與工程的發展和實踐影響最大的10大演算法】，看到了 FFT 演算法，於是寫點。
Fast Fourier Transform (FFT) 是一種用於計算離散傅立葉變換 (DFT) 的演算法。
傅立葉變換是一種將訊號從時間域轉換到頻率域的數學變換，而FFT透過最佳化計算過程使這一變換更加高效。
快速傅立葉變換
考慮以下場景，如何高效的解決？
Q：假設你有一段音樂，它由多個頻率成分組成。
A：FFT可以幫助你找到這些成分的頻率和強度。
比如說，你可以透過FFT發現這段音樂中主要的頻率成分是100Hz、200Hz和300Hz，且它們的強度分別是50、30和20。
什麽是傅立葉變換？
傅立葉變換將復雜的訊號分解成一系列簡單的正弦波（即不同頻率的正弦和余弦波）。
這類似於將一首復雜的交響樂分解成單獨的樂器聲部，使我們能夠分析每種頻率成分的強度。
什麽是快速傅立葉變換？
FFT，英文全稱，Fast Fourier Transform。
傅立葉變換的直接計算（即DFT）非常耗時，特別是當訊號長度較長時。
FFT是一種最佳化演算法，使傅立葉變換的計算變得更快。
它利用訊號的對稱性和周期性，將計算復雜度從 O ( N 2) 降低到 O ( N log N ) 。
電力系統的套用場景
一個完整的視覺化範例
這是一個使用快速傅立葉變換（FFT）進行音訊雜訊消除的範例，用於模擬和實作音訊訊號的降噪系統。
1、匯入包
import matplotlib.pyplot as pltimport numpy as np
2、建立一個包含多個頻率成分的訊號，並為訊號添加了雜訊
# 設定圖形大小和字型大小plt.rcParams['figure.figsize']=[16,12]plt.rcParams.update({'font.size':18})# 建立時間軸dt =0.001t = np.arange(0,1, dt)# 建立包含10種不同頻率的訊號f1 = np.sin(2* np.pi *20* t)f2 = np.sin(2* np.pi *50* t)f3 = np.sin(2* np.pi *70* t)f4 = np.sin(2* np.pi *90* t)f5 = np.sin(2* np.pi *120* t)f6 = np.sin(2* np.pi *150* t)f7 = np.sin(2* np.pi *180* t)f8 = np.sin(2* np.pi *210* t)f9 = np.sin(2* np.pi *240* t)f10 = np.sin(2* np.pi *270* t)# 合並多個頻率f_clean = f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9 + f10# 為訊號添加雜訊f_noisy = f_clean +4* np.random.randn(len(t))
3、生成幹凈訊號的音訊形式
# 匯入音訊播放功能from IPython.display import Audio# 擴充套件訊號f_clean_extended = np.tile(f_clean, 100)# 轉換訊號數據型別f_clean_audio_data = f_clean_extended.astype(np.float32)# 播放音訊display(Audio(f_clean_audio_data , autoplay=True, rate=44100))
4、將包含白雜訊的訊號轉換為音訊形式
# 將訊號轉換為音訊數據f_noisy_extended = np.tile(f_noisy, 100)f_noisy_audio_data = f_noisy_extended.astype(np.float32)display(Audio(f_noisy_audio_data, autoplay=True, rate=44100))
5、繪制幹凈訊號和包含白雜訊的訊號的時域波形圖
# 繪制幹凈和雜訊訊號的時域波形圖plt.plot(t , f_noisy, color='c', linewidth=1.5, label='Noisy')plt.plot(t, f_clean, color='k', linewidth=2, label='Clean')# 設定繪圖限制和標簽plt.xlim(t[0], t[-1]) # 設定x軸的範圍為從t的第一個值到最後一個值plt.legend() #顯示圖例plt.show()
6、進行頻譜分析，並繪制包含雜訊訊號的時域和頻域圖
n = len(t)fhat = np.fft.fft(f_noisy,n)PSD = fhat * np.conj(fhat)/ n #power spectral densityfreq =(1/(dt*n))* np.arange(n)L = np.arange(1,np.floor(n/2), dtype ='int')fig,axs = plt.subplots (2,1)plt.sca(axs[0])plt.plot(t,f_noisy,color='c', label='Noisy')plt.plot(t,f_clean,color='k', label='Clean')plt.xlim(t[0],t[-1])plt.legend()plt.sca(axs[1])plt.plot(freq[L],PSD[L],color='c', label='Noisy')plt.xlim(freq[L[0]], freq[L[-1]])plt.legend()plt.show()
7、利用功率譜密度（PSD）清除訊號中的雜訊成分
# 使用功率譜密度（PSD）清除雜訊indices = PSD > 150 # 找出所有功率大的頻率PSDclean = PSD * indices # 將其他頻率清零fhat = indices * fhat # 將Y中小的Fourier系數清零ffilt = np.fft.ifft(fhat) # 反FFT
# 將訊號轉換為音訊數據ffilt_extended = np.tile(ffilt , 100)ffilt_audio_data = ffilt_extended.astype(np.float32)display(Audio(ffilt_audio_data, autoplay=True, rate=44100))
8、繪制多個子圖，展示不同訊號處理步驟的結果
# 繪制多個子圖fig,axs = plt.subplots(3,1)plt.sca(axs[0])plt.plot(t, f_noisy, color='c', label='Noisy')# 繪制包含雜訊的訊號時域圖plt.plot(t, f_clean, color='k', label='Clean')# 繪制幹凈訊號時域圖plt.xlim(t[0], t[-1])# 設定x軸的範圍為從時間陣列t的第一個值到最後一個值plt.legend()# 顯示圖例plt.sca(axs[1])plt.plot(t, ffilt, color='k', label='Filtered')# 繪制濾波後的訊號時域圖plt.xlim(t[0], t[-1])# 設定x軸的範圍為從時間陣列t的第一個值到最後一個值plt.legend()# 顯示圖例plt.sca(axs[2])plt.plot(freq[L], PSD[L], color='c', label='Noisy')# 繪制包含雜訊訊號的頻域圖plt.plot(freq[L],PSDclean[L], color='k', label='Filtered')# 繪制經過濾波後訊號的頻域圖plt.xlim(freq[L[0]], freq[L[-1]])# 設定x軸的範圍為從頻率陣列L的第一個值到最後一個值plt.legend()# 顯示圖例plt.show() # 顯示繪制的圖形
也可以參看之前文章： 傅立葉變換及Python實作
