numpy实现快速傅里叶变换FFT

参考资料:
https://zhuanlan.zhihu.com/p/31584464
https://blog.csdn.net/weixin_39591031/article/details/110392352
https://blog.csdn.net/weixin_40146921/article/details/122006305

NumPy模块提供了快速傅里叶变换（FFT）的功能。FFT是一种高效计算离散傅里叶变换（DFT）及其逆变换的算法。它能够快速计算时间信号在频域上的表达，在NumPy中，进行快速傅里叶变换可以使用numpy.fft模块。该模块提供了多种函数，用于进行一维和多维的FFT操作。【具体可参看numpy文档】

np.fft.fft(a, n=None,axis=-1,norm=None)

该函数通过快速傅里叶变换算法计算一维(n个元素的)离散傅里叶变换(DFT)

参数：

• a: 输入数组 ，可以有复数
• n: 一个整数，可选参数。它代表选定的轴最后输出的长度，如果n比len(a的那个轴的数组)小，则会进行裁剪。如果更大，则填充0。n默认是len(a的那个轴的数组)
• axis:一个整数，可选参数。表示选定输入数组a的哪个轴，默认选最后一个轴
• norm:规范化方法，可选参数。{“backward”, “ortho”, “forward”}，默认"backward"，表明变换的后向会被放缩
  输出： 带复数的n个元素的数组 shape为(1,n)

如果x是一个1维数组，那么FFT就相当于：
y[k] = np.sum(x * np.exp(-2j * np.pi * k * np.arange(n)/n))

• $y_k={\sum }_{最后对x的每一个元素求和}x\cdot e^{\frac{-2\pi j\cdot k\cdot [0,1,...,n-1]}{n}}$

从上式可知频率 $f=k/n$
当k达到$n/2$时，$y_{n/2}$时已是Nyquist极限，这时进行转向，从k为最小的负数($-n/2$)开始求起。

比如一个x有8个元素，那么它对应的k的计算顺序如下：
[0, 1, 2, 3, -4, -3, -2, -1]。你也可以手动进行一个平移fftshift，使得从小到大。 [-4, -3, -2, -1, 0, 1, 2, 3]

在这个fft算法中，将一个实数序列x(n)映射到一个复数序列X(n)上，输出的复数序列代表了输入实数序列的频域表示。其中每个复数X(k)由实部和虚部组成，复数的虚部表示了输入信号中存在的一些非实部成分，例如一些高频噪声或者相位信息。
虽然FFT变换的结果是复数，但通常情况下我们只关注实部或幅度信息，因为这些信息对于很多信号处理任务来说是足够的。虚部通常被忽略或者用于一些特殊的应用，例如相位恢复或相位校准。

使用示例：

import numpy as np
import matplotlib
matplotlib.use('TkAgg')
matplotlib.rcParams['font.family'] = 'SimHei'
matplotlib.rcParams['font.size'] = 10
matplotlib.rcParams['axes.unicode_minus']=False
import matplotlib.pyplot as plt

# 画出原始信号
plt.figure(figsize=(8, 4))
plt.subplot(3, 1, 1)
plt.plot(np.linspace(0, 0.5, 200), np.cos(2 * np.pi * 5 * np.linspace(0, 0.5, 200)), label='原始信号')
plt.xlabel('时间')
plt.ylabel('幅度')
plt.title('原始信号时域表示')
plt.legend()

# 采样原始信号并进行FFT变换
N = 20  # 采样20个点
t = np.linspace(0, 0.5, N)  # 时间是0-0.5s，共有N个样点，则采样频率为 (N-1) / 0.5 Hz
signal = np.cos(2 * np.pi * 5 * t)  # 原始信号是频率为5，幅度为1的余弦波信号，假如说我们能采集到它这么N个点

# 进行傅里叶变换
fft_signal = np.fft.fft(signal)
# 计算频率轴
freq = np.fft.fftfreq(signal.shape[-1])
# freq = np.fft.fftfreq(signal.shape[-1], (t[-1] - t[0]) / (N - 1))
print(freq)  # f = k / n,  k = 0, 1, ... n/2, -n/2, -n/2+1, ...]
print(fft_signal)
# 绘制频域图

plt.subplot(3, 1, 2)
plt.plot(np.fft.fftshift(freq), np.fft.fftshift(np.abs(fft_signal)), label='频域表示')
plt.xlabel('频率')
plt.ylabel('幅度')
plt.title('频域图')
plt.legend()

signal_reconstructed = np.fft.ifft(fft_signal)  # 是复数，绘图的时候只取其实数部分
print(signal_reconstructed)
plt.subplot(3, 1, 3)
plt.plot(t, np.real(signal_reconstructed), label='重建的原始信号')
plt.legend()
plt.xlabel('时间')
plt.ylabel('幅度')
plt.title('重建的原始信号时域表示')
plt.tight_layout()
plt.show()

在频域图中，f为0.15和-0.15时(即k为3和-3时)，振幅最大为8.51 (3.86-7.58j)。

频率公式
在知道了采样频率Fs后，快速傅里叶变换(FFT)后的第x个（x从0开始）复数值对应的实际频率为
f(x) = x * (Fs / n)
对于实数信号，有对称性，只需要0 ~ N/2 这一半的频率，另一半与这一半对称

也就是说实际上算出来的频率是f = 0.15 * Fs= 0.15 * 38 = 5.7，这虽然不是实际频率5，但是非常接近。主要是我这里采样区间取得太小了，导致频率分辨率低。
我简单说明一下为什么：首先采样频率至少是5*2=10，我显然满足了。但我取的时间是0-0.5s，取N个点，这样我的采样频率是Fs = (N-1)/0.5=2N-1，而f(k) = k * (Fs / N) = 2k-k/N，这k总是取整的，所以根本不会出现f(k)=5的情况，最接近的就是f(3)=6-3/N了，关键就是0.5的区间取小了。

而对于频域中振幅与时域中振幅的关系如下：

则上述转换为时域后最大振幅应为8.51 * (20 / 0.5) = 0.851

也就是说频域的代码可以修改成这样:
plt.plot(np.fft.fftshift(freq * ((N-1)/0.5) ), np.fft.fftshift(np.abs(fft_signal)/ (N/2) ), label='频域表示') 【因为这里没有直流分量(表达式中无常数分量)，就直接忽略了】
或者求freq的时候使用参数采样时间间隔d: freq = np.fft.fftfreq(signal.shape[-1], d= 1.0/ Fs)

那么如何从FFT变换后的频域中提取主要的频率呢？
这里提供一个简单的提取最大频率的示例，使用了周期谱:

import numpy as np
import scipy
import scipy.io
from scipy.signal import butter
def select_max_freq(ori_signal, fs):  # 选出频域图中一小段一小段面积中最大的那个freq
    """
        using Fast Fourier transform (FFT).
    """
    ori_signal = np.expand_dims(ori_signal, 0)  # shape由(n,) 到 (1, n)   转换成了一个行向量
    # 可以这么理解将二维的频域图(frequency, amplitude)压缩到了一维数组中[freq1_amplitude, freq2_amplitude, ...]
    f, pxx = scipy.signal.periodogram(ori_signal, fs=fs,  detrend=False)
    # 从功率密度中挑出最大的那个对应的频率
    max_freq = np.take(f, np.argmax(pxx)) 
    return max_freq


max_freq = select_max_freq(signal, Fs)
print(max_freq)

最终代码如下：

import numpy as np
import matplotlib
import scipy
import scipy.io
from scipy.signal import butter

# matplotlib.use('TkAgg')
matplotlib.rcParams['font.family'] = 'SimHei'
matplotlib.rcParams['font.size'] = 10
matplotlib.rcParams['axes.unicode_minus'] = False
import matplotlib.pyplot as plt

# 画出原始信号
plt.figure(figsize=(8, 4))
plt.subplot(3, 1, 1)
plt.plot(np.linspace(0, 0.5, 200), np.cos(2 * np.pi * 5 * np.linspace(0, 0.5, 200)), label='原始信号')
plt.xlabel('时间')
plt.ylabel('幅度')
plt.title('原始信号时域表示')
plt.legend()

# 采样原始信号并进行FFT变换
N = 20  # 采样20个点
start_time = 0
stop_time = 0.5
t = np.linspace(start_time, stop_time, N)  # 时间是0-0.5s，共有N个样点，则采样频率为 (N-1) / 0.5 Hz
signal = np.cos(2 * np.pi * 5 * t)  # 原始信号是频率为5，幅度为1的余弦波信号，假如说我们能采集到它这么N个点

Fs = (N - 1) / (stop_time - start_time)

# 进行傅里叶变换
fft_signal = np.fft.fft(signal)
# 计算频率轴
# freq = np.fft.fftfreq(signal.shape[-1])
freq = np.fft.fftfreq(signal.shape[-1], 1. / Fs)
print(freq)
print(fft_signal)
# 绘制频域图

plt.subplot(3, 1, 2)
plt.plot(np.fft.fftshift(freq), np.fft.fftshift(np.abs(fft_signal)), label='频域表示')
plt.xlabel('频率')
plt.ylabel('幅度')
plt.title('频域图')
plt.legend()

signal_reconstructed = np.fft.ifft(fft_signal)  # 是复数，绘图的时候只取其实数部分
# 如果仍然要使用复数的虚数部分可以这么做
#signal_reconstructed = np.piecewise(signal_reconstructed,
#                                    [np.real(signal_reconstructed) < 0, np.real(signal_reconstructed) >= 0],
#                                    [lambda elem: -np.abs(elem), lambda elem: np.abs(elem)]).astype(np.float64)
print(signal_reconstructed)
plt.subplot(3, 1, 3)
plt.plot(t, np.real(signal_reconstructed), label='重建的原始信号')
plt.legend()
plt.xlabel('时间')
plt.ylabel('幅度')
plt.title('重建的原始信号时域表示')
plt.tight_layout()
plt.show()

def select_max_freq(ori_signal, fs):  # 选出频域图中一小段一小段面积中最大的那个freq
    """
        using Fast Fourier transform (FFT).
    """
    ori_signal = np.expand_dims(ori_signal, 0)  # shape由(n,) 到 (1, n)   转换成了一个行向量
    # 可以这么理解将二维的频域图(frequency, amplitude)压缩到了一维数组中[freq1_amplitude, freq2_amplitude, ...]
    f, pxx = scipy.signal.periodogram(ori_signal, fs=fs,  detrend=False)
    # 从功率密度中挑出最大的那个对应的频率
    max_freq = np.take(f, np.argmax(pxx)) 
    return max_freq


max_freq = select_max_freq(signal, Fs)
print(max_freq)   # 5.7