16. Scipy Tutorial- 快速傅立叶逆变换ifft

逆向快速傅里叶变换(IFFT)的计算原理是将频域(注意频域是复数)数据进行取共轭复数(虚部取反),然后再进行FFT变换,这样便将频域信号转换到时域。因为FFT变换的结果是复数,所以从频域进行FFT变换过来的结果也是复数,而此时只需取复数的实部,便是原时域信号。

signal_samples函数实现产生一个信号,由两个正弦波叠加而成,两个正弦波一个频率为$1Hz$另一个正弦波的频率为$20Hz$,下面的程序通过处理掉高频的信号而留下$1Hz$信号的数据,实现低通滤波,然后使用ifft将频率数据转为时域数据,恢复时域的波形。程序和前一章区别不大,一个是signal_samples函数是用两个正弦波组成的,另外一处不同是滤波,语句代码是F_filtered = F * (abs(f) < 10),新增了ifft即逆向快速傅里叶变换产生新的时域信号数据f_t_filtered

import numpy as np
import matplotlib.pyplot as plt
def signal_samples(t):
    return np.sin(2 * np.pi * 1 * t)  + np.sin(2 * np.pi * 20 *t)

B = 30.0
f_s = 2 * B
delta_f = 0.01
N = int(f_s / delta_f)
T = N / f_s
t = np.linspace(0, T, N)
f_t = signal_samples(t)
fig, axes = plt.subplots(1, 2, figsize=(8, 3), sharey=True)
axes[0].plot(t, f_t)
axes[0].set_xlabel("time (s)")
axes[0].set_ylabel("signal")
axes[1].plot(t, f_t)
axes[1].set_xlim(0, 5)
axes[1].set_xlabel("time (s)")
plt.show()

from scipy import fftpack
F = fftpack.fft(f_t)
f = fftpack.fftfreq(N, 1.0/f_s)
F_filtered = F * (abs(f) < 10)
f_t_filtered = fftpack.ifft(F_filtered)
mask = np.where(f >= 0)

fig, axes = plt.subplots(3, 1, figsize=(8, 6))
axes[0].plot(f[mask], np.log(abs(F[mask])), label="real")
axes[0].plot(B, 0, 'r*', markersize=10)
axes[0].set_ylabel("$\log(|F|)$", fontsize=14)
axes[1].plot(f[mask], abs(F[mask])/N, label="real")

axes[1].set_ylabel("$|F|$", fontsize=14)
axes[2].plot(t, f_t, label='original')
axes[2].plot(t, f_t_filtered.real, color="red", lw = 3, label='filtered')
axes[2].set_xlim(1, 11)
axes[2].set_xlabel("time (s)", fontsize=14)
axes[2].set_ylabel("$|F|$", fontsize=14)
plt.show()

程序执行结果:

上图画出了频率为$1Hz$、$20Hz$的两个正弦波叠加的信号的采样结果。

上图中间子图可以看出在频谱图里有2个峰值,出现在$1Hz$、$20Hz$和信号产生函数 signal_samples里的正弦波的频率相对应,最下面的子图绘制(红线)的是过滤掉频率为$20Hz$的信号频率后经ifft恢复的信号的时域数据,频率为$1Hz$和函数 signal_samples里频率为$1Hz$的正弦波相呼应。