Python?scipy利用快速傅里葉變換實現(xiàn)濾波

更新時間：2024年01月22日 09:02:25 作者：微小冷

這篇文章主要為大家詳細介紹了Python?scipy如何利用快速傅里葉變換實現(xiàn)濾波,文中的示例代碼講解詳細,感興趣的小伙伴可以跟隨小編一起學習一下

fft模塊簡介

scipy官網(wǎng)宣稱，fftpack模塊將不再更新，或許不久之后將被廢棄，也就是說fft將是唯一的傅里葉變換模塊。

Fourier變換極其逆變換在數(shù)學上的定義如下

下表整理出一部分與Fourier變換相關(guān)的函數(shù)，其中FFT為快速Fourier變換(Fast Fourier Transform)；DCT為離散余弦變換(Discrete Cosine Transform)；DST為離散正弦變換(discrete sine transform)，另外，函數(shù)的前綴和后綴有如下含義

i表示逆變換；
2, n分別表示2維和n維

	正變換	逆變換
通用	fft, fft2, fftn	ifft, ifft2, ifftn
實數(shù)域	rfft, rfft2, rfftn	irfft, irfft2, irfftn
厄米對稱	hfft, hfft2, hfftn	ihfft, ihfft2, ihfftn
離散余弦變換	dct, dctn	idct, idctn
離散正弦變換	dst, dstn	idst, idstn
漢克爾變換	fht	ifht
移動零頻	fftshift	ifftshift
DFT采樣頻率	fftfreq	ifftfreq

fft函數(shù)示例

在數(shù)值計算中，一切輸入輸出均為離散值，所以實際上用到的是離散Fourier變換，即DFT，其功能是將離散的采樣變換為一個離散的頻譜分布。

下面將手動創(chuàng)建一組采樣點，并添加一點噪聲，然后通過FFT獲取其頻域信息。

import numpy as np
from scipy import fft

PI = np.pi*2
fs = 60     #采樣頻率
T = 100     #采樣周期數(shù)
N = fs*T    #采樣點
t = np.linspace(0, T, N)
noise = 2 * np.random.randn(*t.shape)
s = 2 * np.sin(PI * t) + 3 * np.sin(22 * PI * t) + noise
F = fft.fft(s)
f = fft.fftfreq(N, 1.0/fs)

其中，t為時間序列，s為模擬的采樣點，F是Fourier變換之后的結(jié)果。但由于fft默認是在復數(shù)域上的，故而可以查看其實部、虛部、模和輻角的值。

下面對采樣點以及Fourier變換的結(jié)果進行繪制

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(2,2,1)
ax.plot(t, s)
ax.set_title("t vs s")
f_abs = np.abs(F)
ax = fig.add_subplot(2,2,2)
ax.plot(f, f_abs)
ax.set_title("fs vs |F|")

xlims = [[0,2], [21,23]]
for i, xlim in enumerate(xlims):
    ax = fig.add_subplot(2,2,3+i)
    ax.stem(f, f_abs)
    ax.set_title("fs vs |F|")
    ax.set_xlim(xlim)

plt.show()

結(jié)果為

即f=1和f=22處被篩選了出來。

濾波

有了這個，就可以在頻域上對數(shù)據(jù)進行濾波，其思路是，對f_abs中的值進行閾值分割，例如，只篩選出低頻部分，然后看一下濾波效果

fig = plt.figure(1)
f_filt = F * (np.abs(f) < 2)
s_filt = fft.ifft(f_filt)
ax = fig.add_subplot()
ax.plot(t, s, lw=0.2)
ax.plot(t, s_filt.real, lw=2)
ax.set_title("threshold=2")
ax.set_xlim([0,10])
plt.show()

效果如下