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()