SciPy 傅立叶变换

对时域信号计算傅立叶变换（Fourier Transformation）以检查其在频域中的行为。傅里叶变换在信号和噪声处理、图像处理、音频信号处理等学科中应用很广。SciPy 提供了 fftpack 模块，用户可以使用该模块计算快速傅里叶变换。

下面我们通过示例来看一下如何使用FFTpack 模块计算傅立叶变换

下面是一个正弦函数的例子，然后我们使用 fftpack 模块计算傅立叶变换。

快速傅立叶变换

让我们详细了解什么是快速傅立叶变换。

一维离散傅立叶变换

一个长度为N的序列 x[n]，我们对其长度值N使用 fft() 函数计算 FFT y[k]。并且使用 ifft() 函数计算逆变换。

import numpy as np
# 从 fftpack 包中导入 fft 和 ifft
from scipy.fftpack import fft, ifft

# 随机生成5数字元素的数组
x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])

# 然后使用fft函数计算傅立叶变换
y = fft(x)
print(y)

运行示例

运行结果如下

[ 4.50000000+0.j           2.08155948-1.65109876j   -1.83155948+1.60822041j
 -1.83155948-1.60822041j   2.08155948+1.65109876j ]

接下来让我们使用 ifft 进行逆变换

yinv = ifft(y)

print(yinv)

运行示例

运行结果如下

[ 1.0+0.j  2.0+0.j  1.0+0.j -1.0+0.j  1.5+0.j]

scipy.fftpack模块允许计算快速傅立叶变换。例如，我们有一个输入信号是三个正弦波的叠加

import numpy as np
time_step = 0.02
period = 5.
time_vec = np.arange(0, 20, time_step)
sig =7*np.sin(2*np.pi*200*time_vec) + 5*np.sin(2*np.pi*400*time_vec)+3*np.sin(2*np.pi*600*time_vec)

上述示例创建一个时间步长为0.02秒的信号，然后我们打印一下 sig 的大小

print(sig)

运行结果如下

1000

此时我们不知道信号频率；我们只知道信号sig的采样时间步长。信号应该来自实函数，因此傅立叶变换将是对称的。然后我们使用scipy.fftpack.fftfreq()函数计算采样频率，使用scipy.fftpack.fft()计算快速傅立叶变换。

from scipy import fftpack
sample_freq = fftpack.fftfreq(sig.size, d = time_step)
sig_fft = fftpack.fft(sig)
print(sig_fft)

运行结果如下

[  2.21575538e-09 +0.00000000e+00j   4.03133698e-10 +7.37393674e-10j
   4.52465238e-10 +4.70407906e-10j   3.07544640e-10 +5.58941973e-10j
   1.34328382e-10 +6.56307244e-10j  -2.42793467e-10 +2.23010863e-10j
   7.49744033e-12 +4.72055738e-10j  -1.75221507e-10 +1.39737727e-11j
  -2.08093400e-10 +1.82583141e-10j   5.85773449e-10 -4.80043066e-10j
   8.42668121e-11 +9.46190469e-10j  -5.01835896e-10 -3.49354610e-11j
  -2.92218896e-11 +1.46160379e-11j   2.28621881e-10 -2.27646026e-10j
   1.20131398e-11 +6.11882396e-11j   5.01082439e-11 +1.54999502e-10j
  -1.97902675e-11 +3.06585824e-10j   9.75383035e-12 +2.21351889e-10j
   3.81423410e-11 -3.31741298e-10j   5.56370204e-11 -3.71389590e-11j
....]

完整示例如下，我们使用 matplotlib 将信号图形打印出来

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl

mpl.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
mpl.rcParams['axes.unicode_minus'] = False  # 显示负号

time_step = 0.02
period = 5.
x = np.arange(0, 20, time_step)
sig = 7 * np.sin(2 * np.pi * 200 * x) + 5 * np.sin(2 * np.pi * 400 * x) + 3 * np.sin(2 * np.pi * 600 * x)

sample_freq = fftpack.fftfreq(sig.size, d=time_step)
sig_fft = fftpack.fft(sig)

plt.rcParams["font.sans-serif"] = "SimHei"
plt.figure()
plt.plot(sample_freq, sig_fft)
plt.title('sig_fft')

plt.show()

结果如下

离散余弦变换

离散余弦变换（DCT）表示的是余弦函数在不同的频率振荡的总和产生的数据点的有限序列。SciPy 提供dct函数计算离散余弦变换，使用idct进行逆变换。

import numpy as np
from scipy.fftpack import dct
print(dct(np.array([4., 3., 5., 10., 5., 3.])))

运行示例

运行结果如下

[ 60.          -3.48476592 -13.85640646  11.3137085    6.          -6.31319305]

逆离散余弦变换从其离散余弦变换 (DCT) 系数重建序列。idct 函数是 dct 函数的反函数。

import numpy as np
from scipy.fftpack import dct,idct
print(idct(np.array([4., 3., 5., 10., 5., 3.])))

运行示例

运行结果如下

[ 39.15085889 -20.14213562  -6.45392043   7.13341236   8.14213562
  -3.83035081]

下面我们使用 matplotlib 将上述示例的图形画出来，这样比只是看一堆数字更直观一些，我们这里说明一下，上述示例的数据并没有什么实际的意义，仅仅为如何使用相应的函数而设。

完整代码如下

from scipy.fftpack import dct,idct
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl

mpl.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
mpl.rcParams['axes.unicode_minus'] = False  # 显示负号

x = np.array([4., 3., 5., 10., 5., 3.])

y = dct(x)

plt.figure()
plt.plot(x, y)
plt.title('dct')

plt.show()

运行结果如下