快速傅里叶变换（FFT）的简单应用——提取信号中的主频率

概述

我们知道，现实生活中能接触到的信号大都可以展开成傅里叶级数的形式，即表现为无数个余弦波（谐波）的叠加。其中，有的谐波振幅比较大，在信号中占主导地位。利用快速傅里叶变换，我们可以很方便地提取出这些起主导作用的谐波的频率（即题中说的主频率），为进一步还原出这些主要信号做准备。

测试信号

首先，生成一个含有大量噪声的信号：

import numpy
t = numpy.arange(-10, 10, 1 / 1e4) # 时间序列
f1, f2 = 0.1, 2e2 
noise = 10 * numpy.random.randn(len(t)) # 用于干扰的噪音信号，把10改成0就可以获得无噪声的原始信号
# 将频率为f1和f2的两个信号以及噪声信号相叠加，生成我们的测试信号
y = 1.3 * numpy.sin(2 * numpy.pi * f1 * t) + 2.5 * numpy.cos(2 * numpy.pi * f2 * t) + noise
import matplotlib.pyplot as plt
plt.figure()
plt.plot(t, y)
plt.xlabel("t")
plt.ylabel("y")
plt.show()

这个信号在时域上大概长这个鬼样子：

可见由于噪声很大，我们几乎看不出主要的信号长啥样。实际上，在没有噪声干扰的时候，它长这样：

利用FFT提取主频率

from scipy.fft import fft, fftfreq
def analyse_freq_and_amp(x: numpy.ndarray, y: numpy.ndarray):
    """
    分析不同各频率谐波的幅值
    :param x: （几乎）等间隔的升序序列
    :param y:
    :return: 频率序列，幅值序列
    """
    n = len(x)
    sample_freq = (n - 1) / (x[-1] - x[0])  # 信号的采样频率
    freqs = fftfreq(n, 1. / sample_freq)[:n // 2]
    amplitudes = 2. / n * numpy.abs(fft(y)[:n // 2])
    return freqs, amplitudes

调用：

freqs, amps = analyse_freq_and_amp(t, y)
plt.figure()
plt.plot(freqs, amps)
plt.xlabel("Freq./Hz")
plt.ylabel("Amplitude")
plt.show()

在生成的如下所示的频域图上，我们可以看出，这个信号有两个频率的谐波具有很大的幅值：

放大可以看出，这两个频率一个在0.1左右，幅值约1.3；另一个在200左右，幅值约2.5，这正是我们测试信号的主频率和对应的振幅。

如果要还原出完整的信号，除了这两个主频外，我们还需要知道两个信号的初始相位。这一部分的内容将在后面的文章介绍。

完整代码

完整代码如下：

# -*- coding: utf-8 -*-
# @Time    : 2021/12/17 14:56
# @Author  : Z.J. Zhang
# @Email   : zijingzhang@mail.ustc.edu.cn
# @File    : ffttool.py
# @Software: PyCharm
import matplotlib.pyplot as plt
import numpy
from scipy.fft import fft, fftfreq


def analyse_freq_and_amp(x: numpy.ndarray, y: numpy.ndarray):
    """
    分析不同各频率谐波的幅值
    :param x: （几乎）等间隔的升序序列
    :param y:
    :return: 频率序列，幅值序列
    """
    n = len(x)
    sample_freq = (n - 1) / (x[-1] - x[0])  # 信号的采样频率
    freqs = fftfreq(n, 1. / sample_freq)[:n // 2]
    amplitudes = 2. / n * numpy.abs(fft(y)[:n // 2])
    return freqs, amplitudes


if __name__ == "__main__":
    t = numpy.arange(-10, 10, 1 / 1e4)
    f1, f2 = .1, 2e2
    noise = 10 * numpy.random.randn(len(t)) # 用于干扰的噪音信号
    # 将频率为f1和f2的两个信号以及噪声信号相叠加，生成我们的测试信号
    y = 1.3 * numpy.sin(2 * numpy.pi * f1 * t) + 2.5 * numpy.cos(2 * numpy.pi * f2 * t) + noise
    freqs, amps = analyse_freq_and_amp(t, y)
    plt.figure()
    plt.plot(t, y)
    plt.xlabel("t")
    plt.ylabel("y")
    plt.figure()
    plt.plot(freqs, amps)
    plt.xlabel("Freq./Hz")
    plt.ylabel("Amplitude")
    plt.show()