傅里叶分析本质上是一种将函数表示为周期成分之和的方法，并从这些成分中还原原始函数。
当函数及其傅里叶变换都被替换为离散化版本时，称为离散傅里叶变换（DFT）。
离散傅里叶变换能将输入分解成不同频率的组成部分。
- 变换的离散输入通常被称为信号（signal），存在于时域
- 输出则被称为频谱（spectrum） 或 变换（transform），存在于频域
实验步骤如下：
- 首先定义一个标准的方波波形。
- 利用快速傅里叶变化提取主要的频率特征。
- 利用主要的频率特征重建波形，与原始波形对比。
- 绘制三维图分解频率特征，直观感受频率域对原始波形的影响。
import numpy as np # 绘图工具 import matplotlib.pyplot as plt import matplotlib as mpl # 中文乱码处理（从window下找个中文ttf文件，即可） font_path = 'simhei.ttf' mpl.font_manager.fontManager.addfont(font_path) plt.rcParams['font.sans-serif'] = ['simhei'] plt.rcParams['axes.unicode_minus'] = False
标准方波
# 采样率和时长 sr = 16_000 T = 4 # 生成方波 x = np.linspace(0, T, T * sr, endpoint=False) y = 1 * np.sign( np.sin(2 * np.pi * 0.5 * x) ) # 绘制方波 plt.figure(figsize=(5, 3)) plt.plot(x, y) plt.grid(True) plt.ylabel('振幅') plt.xlabel('时间') plt.axhline(y=0, color='k', linestyle='-', linewidth=1) plt.show()
快速傅里叶变换
# 采样总数量 n = len(y) # 快速傅里叶变换 y_fft = np.fft.fft(y) # 计算频率范围 y_freq = np.fft.fftfreq(n, 1 / sr) # 对称的频率取一半 y_fft = y_fft[:n//2] y_freq = y_freq[:n//2] # 信号能量提取 magnitude = 2.0 * np.abs(y_fft) / n
# 绘制频域能量图 plt.figure(figsize=(5, 3)) plt.scatter(y_freq, magnitude, alpha=0.5) plt.title('幅度 vs 频率') plt.xlabel('频率') plt.ylabel('幅度') plt.grid(True) plt.show()
由上图可知：能量主要集中在，低频的那几个频率上。
主要的频率特征
# 这里以幅度的0.1为预置 threshold = 0.1 # 提取幅度大于0.1的频率及其幅度大小 mask = magnitude > threshold detected_freqs = y_freq[mask] detected_amps = magnitude[mask] # 输出主要频率及其能量 print(detected_freqs) print(detected_amps)
[0.5 1.5 2.5 3.5 4.5 5.5]
[1.27323954 0.42441318 0.25464791 0.18189136 0.14147105 0.11574904]
# 绘制能量及其频率关系图 plt.figure(figsize=(5, 3)) plt.plot(detected_freqs, detected_amps) plt.show()
重建波形
根据影响最大的几个频率幅度，重建波形；函数如下：
# 重建波形 re_wave = np.zeros_like(y, dtype=np.float64) for freq, amp in zip(detected_freqs, detected_amps): re_wave += amp * np.sin(2 * np.pi * freq * x) # 重建波形 与 原始波形对比 plt.figure(figsize=(5, 3)) plt.plot(x, y, label='原始信号') plt.plot(x, re_wave, label='重建信号') plt.grid(True) plt.legend() plt.show()
调小幅度阈值threshold，利用更多的频率特征及其幅度，重建信号更加逼近原始信号。
三维图分解频率特征
# 创建3d图形 fig = plt.figure(figsize=(20, 10)) ax = fig.add_subplot(111, projection='3d') # 原始波形 ax.plot(x, y, zs=-3, zdir='y', label='原始信号') # 近似整合波形 rebuild_wave = np.zeros_like(y, dtype=np.float64) # 分频波形 for freq, amp in zip(detected_freqs, detected_amps): if freq > 0: component = amp * np.sin(2 * np.pi * freq * x) rebuild_wave += component ax.plot( x, component, zs=freq, zdir='y', label=f'频率={freq},幅度={amp}' ) ax.plot(x, rebuild_wave, zs=-1, zdir='y', label='重建信号') # 基本配置 ax.set_xlabel('时间') ax.set_ylabel('频率') ax.set_zlabel('幅度') ax.set_title('三维波形分解图') # 增强配置 ax.view_init(elev=30., azim=-25.) # 显示图例 ax.legend() plt.tight_layout() plt.show()
如上图：
- 第一条是原始波形
- 第二条是重建波形
- 之后的波形是对应频率刻度上的特征波形
