傅里叶变换基础
采样定理：采样频率大于信号中最高频率的2倍时，采样后的数字信号可以完整地恢复出原始信号。一般实际应用中保证采样频率为信号最高频率的2.56~4倍。
假设信号频率为F，采样频率为Fs，采样点数为N。为了方便进行FFT运算，通常N取2的整数次幂。
N个采样点，经过fft变换后的结果为N个复数，每个复数对应一个频率(第n<=N/2个点对应的频率为(n-1)/N*Fs)，该复数的模值表示该频率的振幅特征。该振幅特征和原始信号的振幅之间的关系是：如果原始信号的振幅为A，则fft结果的每个点(除了第一个直流分量点)的模值就是A的N/2倍；而第一个点的模值是直流分量振幅的N倍。
注意：这N个复数点去掉第一个点后剩下的N-1个点是关于其中心共轭对称的，因此实际只需要取前一半点的频谱即可，因为共轭对称的两个点的模值(振幅)相同。
In: np.fft.fft(np.array([1,2,3]))
Out: array([ 6. +0.j , -1.5+0.8660254j, -1.5-0.8660254j])
In: np.fft.fft(np.array([1,2,3,4]))
Out: array([10.+0.j, -2.+2.j, -2.+0.j, -2.-2.j])
numpy实现
我们先模拟一个一维时序信号y，它由2V的直流分量(0Hz)，和振幅为3V，频率为50Hz的交流信号，以及振幅为1.5V，频率为75Hz的交流信号组成：
y = 2 + 3*np.cos(2*np.pi*50*t) + 1.5*np.cos(2*np.pi*75*t)
然后我们采用256Hz的采样频率，总共采样256个点。
傅里叶变换结果如下图所示，no normalization对应的是原始的fft结果，normalization是将fft结果的振幅特征转换为原始信号的振幅，可以看到振幅为2V，3V，1.5V的信号分别被解析了出来。
image.png
import numpy as np
import matplotlib.pyplot as plt
Fs = 256 # 采样频率 要大于信号频率的两倍可恢复信号
t = np.arange(0, 1, 1.0/Fs) # 1s采样Fs个点
F1 = 50 # 信号1的频率
F2 = 75 # 信号2的频率
y = 2 + 3*np.cos(2*np.pi*F1*t) + 1.5*np.cos(2*np.pi*F2*t)
N = len(t) # 采样点数
freq = np.arange(N) / N * Fs
Y1 = np.fft.fft(y) # 复数
print('Y[0]', Y1[0])
print('Y[128]', abs(Y1)[120:137])
Y = Y1 / (N/2) # 换算成实际的振幅
Y[0] = Y[0] / 2
freq_half = freq[range(int(N/2))]
Y_half = Y[range(int(N/2))]
fig, ax = plt.subplots(4, 1, figsize=(12, 12))
ax[0].plot(t, y)
ax[0].set_xlabel('Time (s)')
ax[0].set_ylabel('Amplitude')
ax[1].plot(freq, abs(Y1), 'r', label='no normalization')
ax[1].set_xlabel('Freq (Hz)')
ax[1].set_ylabel('Amplitude')
ax[1].legend()
ax[2].plot(freq, abs(Y), 'r', label='normalization')
ax[2].set_xlabel('Freq (Hz)')
ax[2].set_ylabel('Amplitude')
ax[2].set_yticks(np.arange(0, 3))
ax[2].legend()
ax[3].plot(freq_half, abs(Y_half), 'b', label='normalization')
ax[3].set_xlabel('Freq (Hz)')
ax[3].set_ylabel('Amplitude')
ax[3].set_yticks(np.arange(0, 3))
ax[3].legend()
# plt.show()
plt.savefig('a.png')
plt.close()
网友评论