FFT(快速傅里葉變換）示例

#FFT變換是針對一組數值進行運算的，這組數的長度N必須是2的整數次冪，例如64, 128, 256等等； 數值可以是實數也可以是複數，通常我們的時域信號都是實數，因此下面都以實數為例。我們可以把這一組實數想像成對某個連續信號按照一定取樣週期進行取樣而得來，如果對這組N個實數值進行FFT變換，將得到一個有N個複數的數組，我們稱此複數數組為頻域信號，此複數數組符合如下規律：#其結果數組有以下特點：

#下標為0和N/2的兩個複數的虛數部分為0，

#下標為i和N-i的兩個複數共軛，也就是其虛數部分數值相同、符號相反

#首先下標為0的實數表示了時域信號中的直流成分的多少

#下標為i的複數a+b*j表示時域信號中週期為N/i個取樣值的正弦波和餘弦波的成分的多少，

#其中a表示cos波形的成分，b表示sin波形的成分

import numpy as np

import matplotlib

import matplotlib.pyplot as plt

pi = np.pi

time_len = 2.0 #時長

N = 2000 #數據點數,須為偶數，FFT的要求

fs = N/ time_len #[Hz] 取樣頻率

f = np.arange(0,N//2) * fs / N #實際頻率

t = np.arange(0,time_len, time_len/N)

y = 5*np.sin(2*pi*20*t)+8*np.sin(2*pi*40*t) +14.14*np.sin(2*pi*100*t) +14.14*np.cos(2*pi*100*t)+ 16

yf = np.fft.fft(y)# 振幅

bias = (yf[0] / N).real #直流分量,虛部為0

if bias<1e-15: bias = 0

yf_scaled = 2* np.abs(yf) / N

yf_scaled[0] = bias #直流分量(0 Hz處)修正

yf_scaled_half = yf_scaled[:N//2]

#print(np.where(yf_scaled>0.001)])

#np.angle

matplotlib.rcParams["font.sans-serif"] = ["SimHei"]

matplotlib.rcParams["axes.unicode_minus"] = False

plt.subplot(211)

plt.plot(t, y)

plt.xlabel("time 時間")

plt.ylabel("信號值")

plt.grid()

plt.title("時域信號")

plt.subplot(212)

plt.plot(f, yf_scaled_half, "r-")

plt.xlabel("Frequency 頻率[Hz]")

#plt.xlim(0,120)

plt.ylabel("Amplitude 幅值")

plt.grid()

plt.title("FFT 頻域信號")

#plt.suptitle("FFT 示例")

plt.tight_layout()

plt.show()