- 实数傅利叶变换
import numpy as np
def demo_fft():
fig = plt.figure()
f = 10
f_s = 100 # Sampling rate
t = np.linspace(0, 1, f_s, endpoint=False)
y = np.sin(f * 2 * np.pi * t)
ax1 = fig.add_subplot(2, 1, 1)
ax1.plot(t, y)
X = np.fft.fft(y)
freqs = np.fft.fftfreq(len(y)) * len(y)
ax2 =fig.add_subplot(2, 1, 2)
spectrum = np.abs(X)/f_s*2
ax2.stem(freqs, spectrum)
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Magnitude')
ax2.set_xlim(-f_s / 2, f_s / 2)
plt.show()
demo_fft()
Figure_1.png
正弦波的频率为10，为什么FFT的结果是正负10呢？
傅利叶变换是应用在复频域的，原始信号是在复数平面上旋转的圆，所以单独一个频率为10的信号，除了位于X轴，它在任何一个时刻的值都是复数。上面的输入信号中只有实数部分，也就是信号合成后只在X轴上振动。这就要求两个频率相同的信号沿相反方向运动，这样虚部值互相抵消，只留实部。
- 虚数傅利叶变换
根据欧拉公式
一个在复平面内频率为10的信号可以表示为：
def demo_fft():
fig = plt.figure()
f = 10
f_s = 100 # Sampling rate
t = np.linspace(0, 1, f_s, endpoint=False)
y = np.exp(1j * f * 2 * np.pi * t)
ax1 = fig.add_subplot(2, 1, 1)
ax1.plot(t, y)
X = np.fft.fft(y)
freqs = np.fft.fftfreq(len(y)) * len(y)
ax2 =fig.add_subplot(2, 1, 2)
spectrum = np.abs(X)/f_s*2
ax2.stem(freqs, spectrum)
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Magnitude')
ax2.set_xlim(-f_s / 2, f_s / 2)
plt.show()
demo_fft()
Figure_1.png
注意上面计算振幅时没有乘以2，其实不应该乘以2，第一个例子是忽略了负频率所以要求乘以2来补足。
- 从采样数据中获取频率
def get_data():
f = 10
f_s = 100 # Sampling rate
t = np.linspace(0, 1, f_s, endpoint=False)
y = np.sin(f * 2 * np.pi * t)
return y
def demo_fft():
fig = plt.figure()
y = get_data()
ax1 = fig.add_subplot(2, 1, 1)
N = len(y)
ax1.plot(y)
X = np.fft.fft(y)
freqs = np.fft.fftfreq(N) * N
plus_fre = slice(int(N/2))
ax2 = fig.add_subplot(2, 1, 2)
spectrum = np.abs(X)/N*2
ax2.stem(freqs[plus_fre], spectrum[plus_fre])
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Magnitude')
plt.show()
demo_fft()
4.截取信号一段进行分析时
def get_data():
f = 10
f_s = 100 # Sampling rate
t = np.linspace(0, 1, f_s, endpoint=False)
y = np.sin(f * 2 * np.pi * t)
return y
def demo_fft():
fig = plt.figure()
n_list = [67, 79, 88, 97, 100]
for i in range(len(n_list)):
sub = slice(0,n_list[i])
y = get_data()[sub]
ax = fig.add_subplot(len(n_list), 1, i+1)
N = len(y)
X = np.fft.fft(y)
# the time of sample,the whole 100 data is 1 second
s_t = N/100*1
freqs = np.fft.fftfreq(N) * N/s_t
plus_fre = slice(int(N/2))
spectrum = np.abs(X)/N*2
ax.stem(freqs[plus_fre], spectrum[plus_fre])
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude')
ax.set_xlim([-2,50])
ax.set_title(f"N={N}", loc='right', pad=-20)
plt.show()
Figure_1.png
由此可以看出采样点数不同时，有些结果频率不是那么纯，这是因为截取后的采样信号已经不是全部由完整周期构成的了。
5.已知采样时间，便可得出信号频率
之前的例子，采样时间是1秒
def get_data(sample_time,f_s):
""""""
f = 10
t = np.linspace(0, sample_time, sample_time*f_s, endpoint=False)
y = np.sin(f * 2 * np.pi * t)
return y
def demo_fft():
sample_time = 2
fig = plt.figure()
n_list = [67, 79, 88, 97, 100]
y_data = get_data(sample_time, 200)
total = len(y_data)
n_list = [int(i/100*total) for i in n_list]
for i in range(len(n_list)):
sub = slice(0,n_list[i])
y = y_data[sub]
ax = fig.add_subplot(len(n_list), 1, i+1)
N = len(y)
X = np.fft.fft(y)
# the time of sample,the whole data are sample_time second
s_t = N/total*sample_time
freqs = np.fft.fftfreq(N) * N/s_t
plus_fre = slice(int(N/2))
spectrum = np.abs(X)/N*2
ax.stem(freqs[plus_fre], spectrum[plus_fre])
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude')
ax.set_xlim([-2,50])
ax.set_title(f"N={N}", loc='right', pad=-20)
plt.show()
demo_fft()
这样改变采样时间，采样频率不会改变结果，但采样频率要足够高。
网友评论