numpy實現快速傅里葉變換

1、快速傅里葉變換的實現

什么是傅里葉定理？

　　法國科學家傅里葉提出，任何一條周期性曲線，無論多么跳躍或不規則，都能表示成一組光滑正弦曲線疊加之和。

什么是傅里葉變換？

　　傅里葉變換即是把一條周期性曲線拆解成一組光滑正弦曲線的過程。

　　傅里葉變換的目的是可將時域（即時間域）上的信號轉變為頻域（即頻率域）上的信號，隨着域的不同， 對同一個事物的認知角度也隨之改變，因此在時域中某些不好處理的地方，在頻域就可以較為簡單地處理。這就可以大量減少處理信號存儲量。

　　假設有一時間域函數：y = f(x)，根據傅里葉的理論它可以被分解為一系列正弦函數的疊加，他們的振幅A，頻率ω或初相位φ不同：

　　所以傅里葉變換可以把一個比較復雜的函數轉換為多個簡單函數的疊加，看問題的角度也從時間域轉到了頻率域，有些問題處理起來就會比較簡單。

傅里葉變換相關函數

　　導入快速傅里葉變換所需模塊

import numpy.fft as nf

　　通過采樣數域采樣周期求得傅里葉變換分解所得曲線的頻率序列

freqs = nf.fftfreq(采樣熟練，采樣周期)

　　通過原函數值得序列經過快速傅里葉變換得到一個復數數組，復數的模代表的是振幅，復數的輻角代表初相位

nf.fft(原函數序列值) -> 目標函數序列值（復數數組）

　　通過一個復數數組（復數的模代表振幅，復數的輻角代表初相位）經過逆傅里葉變換得到合成的函數值數組

nf.ifft(目標函數值序列（復數）) -> 原函數值序列

　　案例：

import numpy as np import matplotlib.pyplot as plt import numpy.fft as nf x = np.linspace(-2 * np.pi, 2 * np.pi, 1000) y = np.zeros(x.size) for i in range(1, 1000): y += 4 * np.pi / (2 * i - 1) * np.sin((2 * i - 1) * x) plt.figure('FFT', facecolor='lightgray') plt.subplot(121) plt.title('Time Domain', fontsize=16) plt.grid(linestyle=':') plt.plot(x, y, label=r'$y$') # 針對方波y做fft
comp_arr = nf.fft(y) y2 = nf.ifft(comp_arr).real plt.plot(x, y2, color='orangered', linewidth=5, alpha=0.5, label=r'$y$') # 繪制頻域圖形
plt.subplot(122) freqs = nf.fftfreq(y.size, x[1] - x[0]) pows = np.abs(comp_arr)  # 復數的模
plt.title('Frequency Domain', fontsize=16) plt.grid(linestyle=':') plt.plot(freqs[freqs > 0], pows[freqs > 0], color='orangered', label='frequency') plt.legend() plt.savefig('fft.png') plt.show()

　　運行結果：

2、基於傅里葉變換的頻域濾波

　　含噪信號是高能信號與低能噪聲疊加的信號，可以通過傅里葉變換的頻域濾波實現降噪。

　　通過FFT使含噪信號轉換為含噪頻譜，去除低能噪聲，留下高能頻譜后再通過IFFT留下高能信號。

　　案例：

import numpy as np import numpy.fft as nf import matplotlib.pyplot as plt import scipy.io.wavfile as sw # 獲取采樣率、采樣位移
sample_rate, sigs = sw.read('noised.wav') print(sample_rate, sigs.shape) sigs = sigs / (2 ** 15) times = np.arange(sigs.size) / sample_rate # 繪圖
plt.figure('Filter', facecolor='lightgray') plt.subplot(221) plt.title('Time Domain', fontsize=16) plt.ylabel('Signal', fontsize=14) plt.grid(linestyle=':') plt.plot(times[:178], sigs[:178], color='dodgerblue', label='Nosied Signal') plt.legend() # fft變換后繪制頻域圖
comp_arr = nf.fft(sigs) pows = np.abs(comp_arr) freqs = nf.fftfreq(sigs.size, times[1] - times[0]) plt.subplot(222) plt.title('Frequency Domain', fontsize=16) plt.ylabel('pow', fontsize=14) plt.grid(linestyle=':') plt.semilogy(freqs[freqs > 0], pows[freqs > 0], color='orangered', label='Freq') plt.legend() # 在頻譜中去除噪聲
maxpow_freq = freqs[pows.argmax()] noised_inds = np.where(freqs != maxpow_freq) # 噪聲的復數全改為0
comp_arr[noised_inds] = 0 # 繪制降噪之后的頻譜圖形
pows = np.abs(comp_arr) plt.subplot(224) plt.ylabel('pow', fontsize=14) plt.grid(linestyle=':') plt.plot(freqs[freqs > 0], pows[freqs > 0], color='orangered', label='Freq') plt.legend() # 逆向傅里葉變換，輸出時域函數圖形
filter_sigs = nf.ifft(comp_arr).real plt.subplot(223) plt.ylabel('Signal', fontsize=14) plt.grid(linestyle=':') plt.plot(times[:178], filter_sigs[:178], color='dodgerblue', label='Filter Signal') plt.legend() # 保存濾波后的信號
sw.write('filter.wav', sample_rate, (filter_sigs * 2 ** 15).astype(np.int16)) plt.tight_layout() plt.savefig('filter.png') plt.show()

　　運行結果：