概述

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

测试信号

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

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()