python作傅里葉變換

傅里葉變換(fft)

　　法國科學家傅里葉提出，任何一條週期曲線，不管多麼跳躍或不規則，都能表示成一組光滑正弦曲線疊加之和。傅里葉變換便是把一條不規則的曲線拆解成一組光滑正弦曲線的過程。

　　傅里葉變換的目的是將時域（即時間域）上的信號轉變爲頻域（即頻率域）上的信號，隨着域的變換，對同一個事物的瞭解角度也就隨之改變，所以在時域中某些很差處理的地方，在頻域就能夠較爲簡單的處理。這就能夠大量減小處理信號存儲量。

例如：彈鋼琴

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

$$y = A_1sin(\omega_1x+\phi_1) +  A_2sin(\omega_2x+\phi_2) +  A_2sin(\omega_2x+\phi_2) + R$$

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

傅里葉變換相關函數

freqs = np.fft.fft.fftfreq(採樣數量, 採樣週期)　　經過採樣數與採樣週期求得傅里葉變換分解所得數組的頻率序列

np.fft.fft.fft(原序列)　　原函數值的序列通過快速傅里葉變換獲得一個複數數組，複數的模表明的是振幅，複數的輻角表明初相位

np.fft.fft.ifft(複數序列)　　複數數組 通過逆向傅里葉變換獲得合成的函數值數組

案例：針對合成波作快速傅里葉變換，獲得分解波數組的頻率、振幅、初相位數組，並繪製頻域圖像。

import numpy as np
import matplotlib.pyplot as plt
import numpy.fft as fft

x = np.linspace(-2 * np.pi, 2 * np.pi, 1000)

n = 1000
y = np.zeros(x.size)
for i in range(1, n + 1):
    y += 4 * np.pi / (2 * i - 1) * np.sin((2 * i - 1) * x)

complex_array = fft.fft(y)
print(complex_array.shape)  # (1000,) 
print(complex_array.dtype)  # complex128 
print(complex_array[0])  # (-2.1458390619955026e-12+0j)
y_new = fft.ifft(complex_array)

plt.subplot(311)
plt.grid(linestyle=':')
plt.plot(x, y, label='y')  # y是1000個相加後的正弦序列
plt.subplot(312)
plt.plot(x, y_new, label='y_new', color='orangered')  # y是ifft變換後的序列

# 獲得分解波的頻率序列
freqs = fft.fftfreq(x.size, x[1] - x[0])
# 複數的模爲信號的振幅（能量大小）
complex_array = fft.fft(y)
pows = np.abs(complex_array)

plt.subplot(313)
plt.title('Frequency Domain', fontsize=16)
plt.xlabel('Frequency', fontsize=12)
plt.ylabel('Power', fontsize=12)
plt.tick_params(labelsize=10)
plt.grid(linestyle=':')
plt.plot(freqs[freqs > 0], pows[freqs > 0], c='orangered', label='Frequency')
plt.legend()
plt.tight_layout()
plt.show()

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

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

經過FFT使含噪信號轉換爲含噪頻譜，去除低能噪聲，留下高能頻譜後再經過IFFT留下高能信號。

案例：基於傅里葉變換的頻域濾波爲音頻文件去除噪聲(noiseed.wav數據集地址)。

　　一、讀取音頻文件，獲取音頻文件基本信息：採樣個數，採樣週期，與每一個採樣的聲音信號值。繪製音頻時域的：時間/位移圖像

import numpy as np
import numpy.fft as nf
import scipy.io.wavfile as wf
import matplotlib.pyplot as plt

# 讀取音頻文件
sample_rate, noised_sigs = wf.read('./da_data/noised.wav')
print(sample_rate)  # sample_rate：採樣率44100
print(noised_sigs.shape)    # noised_sigs:存儲音頻中每一個採樣點的採樣位移(220500,)
times = np.arange(noised_sigs.size) / sample_rate

plt.figure('Filter')
plt.subplot(221)
plt.title('Time Domain', fontsize=16)
plt.ylabel('Signal', fontsize=12)
plt.tick_params(labelsize=10)
plt.grid(linestyle=':')
plt.plot(times[:178], noised_sigs[:178], c='orangered', label='Noised')
plt.legend()

 

　　二、基於傅里葉變換，獲取音頻頻域信息，繪製音頻頻域的：頻率/能量圖像

# 傅里葉變換後，繪製頻域圖像
freqs = nf.fftfreq(times.size, times[1] - times[0])
complex_array = nf.fft(noised_sigs)
pows = np.abs(complex_array)

plt.subplot(222)
plt.title('Frequency Domain', fontsize=16)
plt.ylabel('Power', fontsize=12)
plt.tick_params(labelsize=10)
plt.grid(linestyle=':')
# 指數增加座標畫圖
plt.semilogy(freqs[freqs > 0], pows[freqs > 0], c='limegreen', label='Noised')
plt.legend()

　　三、將低頻噪聲去除後繪製音頻頻域的：頻率/能量圖像

# 尋找能量最大的頻率值
fund_freq = freqs[pows.argmax()]
# where函數尋找那些須要抹掉的複數的索引
noised_indices = np.where(freqs != fund_freq)
# 複製一個複數數組的副本，避免污染原始數據
filter_complex_array = complex_array.copy()
filter_complex_array[noised_indices] = 0
filter_pows = np.abs(filter_complex_array)

plt.subplot(224)
plt.xlabel('Frequency', fontsize=12)
plt.ylabel('Power', fontsize=12)
plt.tick_params(labelsize=10)
plt.grid(linestyle=':')
plt.plot(freqs[freqs >= 0], filter_pows[freqs >= 0], c='dodgerblue', label='Filter')
plt.legend()

　　四、基於逆向傅里葉變換，生成新的音頻信號，繪製音頻時域的：時間/位移圖像

filter_sigs = nf.ifft(filter_complex_array).real
plt.subplot(223)
plt.xlabel('Time', fontsize=12)
plt.ylabel('Signal', fontsize=12)
plt.tick_params(labelsize=10)
plt.grid(linestyle=':')
plt.plot(times[:178], filter_sigs[:178], c='hotpink', label='Filter')
plt.legend()

　　五、從新生成音頻文件

# 生成音頻文件
wf.write('./da_data/filter.wav', sample_rate, filter_sigs)
plt.show()