SciPy - FFTpack



傅立葉變換是在時域訊號上計算的,用於檢查其在頻域中的行為。傅立葉變換應用於訊號和噪聲處理、影像處理、音訊訊號處理等學科。SciPy 提供了 fftpack 模組,允許使用者計算快速傅立葉變換。

以下是一個正弦函式的示例,它將用於使用 fftpack 模組計算傅立葉變換。

快速傅立葉變換

讓我們詳細瞭解什麼是快速傅立葉變換。

一維離散傅立葉變換

長度為 N 的序列 x[n] 的長度為 N 的 FFT y[k] 由 fft() 計算,逆變換由 ifft() 計算。讓我們考慮以下示例

#Importing the fft and inverse fft functions from fftpackage
from scipy.fftpack import fft

#create an array with random n numbers
x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])

#Applying the fft function
y = fft(x)
print y

上述程式將生成以下輸出。

[ 4.50000000+0.j           2.08155948-1.65109876j   -1.83155948+1.60822041j
 -1.83155948-1.60822041j   2.08155948+1.65109876j ]

讓我們看另一個例子

#FFT is already in the workspace, using the same workspace to for inverse transform

yinv = ifft(y)

print yinv

上述程式將生成以下輸出。

[ 1.0+0.j   2.0+0.j   1.0+0.j   -1.0+0.j   1.5+0.j ]

scipy.fftpack 模組允許計算快速傅立葉變換。例如,一個(噪聲)輸入訊號可能如下所示:

import numpy as np
time_step = 0.02
period = 5.
time_vec = np.arange(0, 20, time_step)
sig = np.sin(2 * np.pi / period * time_vec) + 0.5 *np.random.randn(time_vec.size)
print sig.size

我們正在建立一個時間步長為 0.02 秒的訊號。最後一條語句列印訊號 sig 的大小。輸出如下:

1000

我們不知道訊號頻率;我們只知道訊號 sig 的取樣時間步長。該訊號應該來自一個實函式,因此傅立葉變換將是對稱的。scipy.fftpack.fftfreq() 函式將生成取樣頻率,scipy.fftpack.fft() 將計算快速傅立葉變換。

讓我們透過一個例子來理解這一點。

from scipy import fftpack
sample_freq = fftpack.fftfreq(sig.size, d = time_step)
sig_fft = fftpack.fft(sig)
print sig_fft

上述程式將生成以下輸出。

array([ 
   25.45122234 +0.00000000e+00j,   6.29800973 +2.20269471e+00j,
   11.52137858 -2.00515732e+01j,   1.08111300 +1.35488579e+01j,
   …….])

離散餘弦變換

離散餘弦變換 (DCT) 用一系列不同頻率的餘弦函式之和來表示有限的數列資料點。SciPy 使用函式dct 提供 DCT,並使用函式idct 提供相應的 IDCT。讓我們考慮以下示例。

from scipy.fftpack import dct
print dct(np.array([4., 3., 5., 10., 5., 3.]))

上述程式將生成以下輸出。

array([ 60.,  -3.48476592,  -13.85640646,  11.3137085,  6.,  -6.31319305])

逆離散餘弦變換根據其離散餘弦變換 (DCT) 係數重建序列。idct 函式是 dct 函式的逆函式。讓我們透過以下示例來理解這一點。

from scipy.fftpack import dct
print idct(np.array([4., 3., 5., 10., 5., 3.]))

上述程式將生成以下輸出。

array([ 39.15085889, -20.14213562, -6.45392043, 7.13341236,
8.14213562, -3.83035081])
廣告

© . All rights reserved.