fft python 举例_快速傅里叶变换（fft）的矩阵理解和在numpy中的举例说明

本文简单总结一下快速傅里叶变换的矩阵理解角度和在numpy中的语法和使用举例。

1.傅里叶变换的矩阵表示

我们在学习数字信号处理时遇到的离散时间傅里叶变换的公式都是以求和的形式出现的，即：

其中

，公式(1)对一个序列

做

点傅里叶变换。我们将这个公式写成矩阵乘以向量的形式，令

，

，即：

公式(2)用矩阵变量表示为：

其中

称为傅里叶矩阵，它是一个具有特殊结构的范德蒙德矩阵，其(

,

)的元素为

。由定义可知，傅里叶具有以下性质：

(1)

是对称矩阵：

；

(2)容易验证

；

(3)由于

是一个具有特殊结构的范德蒙德矩阵，因此它可逆，且

。

那么由公式(3)可得傅里叶逆变换公式为：

写成求和公式就是：

可见，求得了傅里叶矩阵

，其逆矩阵就是它的共轭再除以

，也就是傅里叶逆变换很容易求出来。

2.快速傅里叶变换

实际中，快速傅里叶变换算法是利用反向指标向量进行迭代得到一个指标反转的傅里叶矩阵，这里不做说明，直接给出公式：

其中：

表示克罗内克积。

则公式(3)变成：

例如，对于4点快速傅里叶变换，公式(8)的结果为：

8点快速傅里叶变换：

对于傅里叶反变换:

公式(8)中

的逆矩阵为：

公式(12)中：

其中：

先不用管上述公式什么得到的，我个人感觉矩阵形式的快速傅里叶变换比数字信号书里讲的蝶形运算好理解很多。

3.numpy中的快速傅里叶变换举例

考虑一个正弦函数，其中幅度

，频率

，初相

，即

采样频率取

(过采样了)。

(1)信号的产生

产生一个正弦信号，其公式是：

其中

,

。

Python3代码如下：

import numpy as np

import matplotlib.pyplot as plt

plt.style.use('seaborn')

A = 0.2

fc = 10

phase = np.pi / 6.0

fs = 32 * fc

t = np.arange(0, 2, 1.0 / fs) # 2 seconds of sampling time

x = A * np.sin(2 * np.pi * fc * t + phase)

plt.plot(t ,x)信号的产生(一个2秒的正弦信号)

(2)在numpy中计算离散时间傅里叶变换

N = 256 # N 点傅里叶变换

X = 1.0 / N * np.fft.fftshift(np.fft.fft(x, N))

函数的意思是将零频率分量移至频谱中心。下面我们比较一下如果不移动的话，频谱图是什么样子的。

# 没有将零频率分量移至频谱中心后的频谱

X_no_shift = 1.0 / N * np.fft.fft(x, N)

将零频率分量移至频谱中心后的频谱图：

plt.plot(freqs, np.abs(X));

plt.title('shifting');

没有将零频率分量移至频谱中心后的频谱图：

plt.plot(freqs, np.abs(X_no_shift));

plt.title('no shifting');

(3)相位谱

相位计算公式是将频谱数据中的虚部除以它的实部，再取反正切即可。

phases = np.arctan2(np.imag(X), np.real(X)) * 180 / np.pi # 虚部除以实部得到相位谱，返回值是弧度

plt.plot(freqs, phases)

相位谱如下：相位谱

可以看到，相位谱全是噪声，为什么会这样呢。我们看下算出来的频域数据是多少，这里给出前5个数据如下：

[-1.12757026e-17+0.00000000e+00j, -4.34079651e-18+1.52642204e-18j,

-6.87570451e-18+7.00872400e-18j, -1.21284986e-17-7.76954426e-18j,

-2.82817394e-18-6.98260581e-18j]

可以看到，这些数据特别小，由于计算机计算的精度问题，在进行反正切时导致数据异常如溢出等。

解决方法是：设置一容差阈值，小于该阈值直接置零。即：

X_backup = np.copy(X)

threshold = np.max(np.abs(X)) / 1000.0

X_backup[X_backup < threshold] = 0.0

phases = np.arctan2(np.imag(X_backup), np.real(X_backup)) * 180 / np.pi # 虚部除以实部得到相位谱，返回值是弧度

plt.stem(freqs, phases, use_line_collection=True)

最后画出来的相位谱是：处理后的相位谱

(4)从频域采样重构时域信号

直接作傅里叶反变换即可：

x_recon = np.fft.ifft(np.fft.ifftshift(X), N)

t = np.arange(0, len(x_recon)) / fs

plt.plot(t, x_recon)从频域信息重构时域信号

4.总结

以上从矩阵的角度简略地讲述了傅里叶变换和快速傅里叶变换，提供一个看待傅里叶变换的角度，如有错误，恳请各位知友批评指正，蟹蟹~~

参考文献

[1] 张贤达-《矩阵分析》

[2] 下列网站：https://www.gaussianwaves.com/2015/11/interpreting-fft-results-obtaining-magnitude-and-phase-information​www.gaussianwaves.com