cpp集成fftw [ Soymilk ]

cpp集成fftw

categories
>
Code
>
CPP

写在前面

FFTW (Fastest Fourier Transform in the West) 是业界公认最快的 FFT 库，支持实数/复数、多维变换、多线程并行、SIMD 优化，广泛用于科学计算、音频处理、图像分析、通信系统等领域。

本文采用 C++ 核心计算 + Python 可视化 分离架构：

C++：FFTW 实数 FFT (r2c/c2r)、低通滤波、数据导出 CSV
Python：Matplotlib 读取 CSV 生成高质量 PNG 图表（频谱、时域对比、滤波器响应、波形放大）
现代化 CMake (find_library) 集成 FFTW

适合读者：需要做频域分析、数字滤波、频谱可视化的 C++ 开发者。

1. 环境准备与编译

1.1 获取源码

wget http://www.fftw.org/fftw-3.3.10.tar.gz
tar xzf fftw-3.3.10.tar.gz
cd fftw-3.3.10

1.2 配置与编译

推荐开启共享库、单精度、SIMD、多线程：

./configure --enable-shared --enable-float --enable-sse2 --enable-avx --enable-threads
make -j$(nproc)
sudo make install
ldconfig

验证安装：

pkg-config --libs fftw3
# 输出: -lfftw3 -lm

Gentoo 用户：emerge sci-libs/fftw 自动处理 USE 标志（threads sse avx float）

2. 现代化 CMake 集成

使用 find_package 替代硬编码路径，支持系统安装与 FetchContent 无缝切换。

# CMakeLists.txt
cmake_minimum_required(VERSION 3.16)
project(FFTWDemo LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)

if(CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang")
    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -O3 -march=native -ffast-math")
endif()

# 查找 FFTW3 库
find_library(FFTW3_LIB fftw3)
find_library(M_LIB m)

if(NOT FFTW3_LIB)
    message(FATAL_ERROR "FFTW3 library not found")
endif()

add_executable(main main.cpp)
target_link_libraries(main PRIVATE ${FFTW3_LIB} ${M_LIB})
target_include_directories(main PRIVATE /usr/include)

# 可选：FetchContent 自动下载编译（离线环境/版本锁定）
# include(FetchContent)
# FetchContent_Declare(fftw
#   GIT_REPOSITORY https://github.com/FFTW/fftw3
#   GIT_TAG        fftw-3.3.10
# )
# FetchContent_MakeAvailable(fftw)

构建：

mkdir build && cd build
cmake .. && make -j
./main

3. 核心代码逐段解析

完整源码见文章代码片段，此处分关键步骤讲解。

3.1 合成信号生成（替代 CSV 读取）

constexpr double fs = 1000.0;  // 采样率 Hz
constexpr int N = 1024;        // 2 的幂，FFTW 效率最高
constexpr double f0 = 50.0;    // 基波
constexpr double f1 = 150.0;   // 3次谐波
constexpr double f2 = 500.0;   // 高频干扰

double* input = (double*)fftw_malloc(sizeof(double) * N);  // 必须用 fftw_malloc 保证对齐

std::mt19937 gen(42);
std::normal_distribution<double> noise(0.0, 0.1);

for (int i = 0; i < N; ++i) {
    double t = i / fs;
    input[i] = std::sin(2*PI*f0*t) + 0.5*std::sin(2*PI*f1*t)
             + 0.3*std::sin(2*PI*f2*t) + noise(gen);
}

关键点：必须用 fftw_malloc/fftw_free 分配输入输出数组，保证 16/32 字节对齐，否则 SSE/AVX 指令会崩溃。

3.2 先创建 Plan，再填充数据（避坑指南）

// 1. 先创建计划（FFTW_MEASURE 会覆盖输入输出数组！）
fftw_plan forward = fftw_plan_dft_r2c_1d(N, input, out, FFTW_MEASURE);
fftw_plan backward = fftw_plan_dft_c2r_1d(N, out, smoothed, FFTW_MEASURE);

// 2. 计划创建后，再填充输入数据
for (int i = 0; i < N; ++i) input[i] = ...;

// 3. 执行变换
fftw_execute(forward);

⚠️ FFTW_MEASURE/PATIENT/EXHAUSTIVE 规划阶段会读写输入输出数组，必须在 fftw_plan_* 之后再初始化数据，否则数据被覆盖导致全零谱。

3.3 实数 FFT 正变换 (`r2c`)

fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * (N/2 + 1));
fftw_plan forward = fftw_plan_dft_r2c_1d(N, input, out, FFTW_MEASURE);
fftw_execute(forward);

// 输出格式：out[0..N/2] 对应频率 0, fs/N, 2fs/N, ..., fs/2 (奈奎斯特频率)
// 共 N/2+1 个复数，利用实信号共轭对称性仅存储非冗余部分
for (int i = 0; i < N/2 + 1; ++i) {
    double re = out[i][0];
    double im = out[i][1];
    double amp = std::sqrt(re*re + im*im);  // 幅值
    double freq = i * fs / N;               // 频率轴
}

3.4 低通滤波（核心修正版）

int cutoff_bin = N / 10;  // 截止频率 bin，对应 fs/10 = 100Hz
// 保留低频 [0, cutoff_bin)，零化高频 [cutoff_bin, N/2]
for (int i = cutoff_bin; i < N/2 + 1; ++i) {
    out[i][0] = 0.0;
    out[i][1] = 0.0;
}

修正说明：原草稿中 for (int i = 0; i < cutoff; ++i) out[i]=0 是高通而非低通。正确逻辑是保留低频分量（i < cutoff_bin）、零化高频分量（i >= cutoff_bin）。

3.5 实数 FFT 反变换 (`c2r`) + 归一化

fftw_plan backward = fftw_plan_dft_c2r_1d(N, out, smoothed, FFTW_MEASURE);
fftw_execute(backward);

// c2r 结果需除以 N 恢复原始幅值
for (int i = 0; i < N; ++i) smoothed[i] /= N;

为什么要除以 N？ FFTW 正变换无归一化，反变换累加 N 次。r2c → c2r 往返等价于乘以 N，需 /N 抵消。

3.6 结果导出 CSV（供外部工具/复现）

// spectrum.csv: Frequency_Hz,Amplitude
// waveform_raw.csv: Time_s,Amplitude
// waveform_filtered.csv: Time_s,Amplitude

4. Python + Matplotlib 可视化（高质量 PNG）

核心 FFT 计算用 C++，数据导出 CSV，再用 Python 生成专业图表，避免 C++ 手写 SVG 的坐标映射、渲染兼容性问题。

# plot_from_csv.py 关键片段
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# 1. 频谱图：峰值标注 + 截止频率虚线
df = pd.read_csv('spectrum.csv')
plt.plot(df['Frequency_Hz'], df['Amplitude'], color='#2c3e50', lw=1.2)
plt.axvline(x=100, color='#e74c3c', ls='--', lw=2, label='Cutoff: 100Hz')
# 标注 50Hz/150Hz 主峰
for f, a in [(49.8, 479.8), (150.4, 191.7)]:
    plt.annotate(f'{f:.0f}Hz\n{a:.0f}', xy=(f, a), xytext=(f+15, a+50),
                 arrowprops=dict(arrowstyle='->', color='#2c3e50'))
plt.savefig('fftw/spectrum.png', dpi=150)

# 2. 时域波形对比：双线同图，透明度区分
df_raw = pd.read_csv('waveform_raw.csv')
df_filt = pd.read_csv('waveform_filtered.csv')
plt.plot(df_raw['Time_s'], df_raw['Amplitude'], color='#e74c3c', lw=0.8, alpha=0.7, label='Raw')
plt.plot(df_filt['Time_s'], df_filt['Amplitude'], color='#27ae60', lw=1.5, alpha=0.95, label='Filtered')
plt.savefig('fftw/waveform_compare.png', dpi=150)

# 3. 滤波器响应：理想矩形窗
freq = np.linspace(0, 500, 513)
gain = np.where(freq < 100, 1.0, 0.0)
plt.plot(freq, gain, color='#8e44ad', lw=2, ls='--', label='Ideal Low-pass')
plt.savefig('fftw/filter_response.png', dpi=150)

# 4. 波形局部放大 (0-0.1s)：双子图展示滤波细节
fig, axes = plt.subplots(2, 1, sharex=True, figsize=(10, 8))
axes[0].plot(df_raw['Time_s'], df_raw['Amplitude'], color='#e74c3c')
axes[1].plot(df_filt['Time_s'], df_filt['Amplitude'], color='#27ae60')
plt.savefig('fftw/waveform_zoom.png', dpi=150)

生成 4 张高质量 PNG (150 DPI)：

图表	文件	说明
频谱图	`fftw/spectrum.png`	幅值谱，峰值标注(50Hz/150Hz)，截止频率虚线
时域对比	`fftw/waveform_compare.png`	红=原始含噪/高频，绿=滤波后平滑
滤波器响应	`fftw/filter_response.png`	矩形窗理想低通响应
波形放大	`fftw/waveform_zoom.png`	0-0.1s 局部，直观展示滤波前后细节

5. 关键点深度解析（折叠面板）

实数 FFT 存储格式：为何只有 N/2+1 个复数？

实信号 FFT 满足共轭对称性：$X[k] = X^*[N-k]$。FFTW 仅存储非冗余部分：

out[0]：DC 分量 (0 Hz)，实部有效，虚部为 0
out[1] ~ out[N/2-1]：正频率分量
out[N/2]：奈奎斯特频率 ($f_s/2$)，实部有效，虚部为 0

负频率部分由共轭对称性隐含，节省一半存储与计算。

归一化因子：为什么 c2r 后要除以 N？

FFTW 定义：

正变换 (DFT)：$X[k] = \sum_{n=0}^{N-1} x[n] e^{-i 2\pi kn/N}$（无系数）
反变换 (IDFT)：$x[n] = \sum_{k=0}^{N-1} X[k] e^{+i 2\pi kn/N}$（无系数）

往返变换：$IDFT(DFT(x)) = N \cdot x[n]$，故需 /N 恢复幅值。
若用 fftw_plan_dft_r2c_1d + fftw_plan_dft_c2r_1d，归一化因子为 $N$。

FFTW_ESTIMATE vs MEASURE vs PATIENT

标志	规划耗时	运行性能	适用场景
`FFTW_ESTIMATE`	极快(启发式)	较低	一次性小规模、原型验证
`FFTW_MEASURE`	中等(实测多种实现)	最优	生产环境推荐
`FFTW_PATIENT`	慢(更彻底搜索)	极优	离线预计算、极致性能
`FFTW_EXHAUSTIVE`	极慢(穷举)	理论最优	研究/基准测试

本文用 FFTW_MEASURE：首次运行约 200ms 规划，后续执行极快（~0ms）。

Wisdom 缓存：二次运行加速规划

// 首次运行后保存
fftw_export_wisdom_to_filename("fftw_wisdom.dat");

// 后续运行加载（跳过 MEASURE 耗时）
fftw_import_wisdom_from_filename("fftw_wisdom.dat");
fftw_plan forward = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); // 瞬间完成

适合固定尺寸重复运行的场景（音频流、实时频谱仪）。

多线程并行加速

fftw_init_threads();                    // 初始化线程系统
fftw_plan_with_nthreads(4);             // 使用 4 线程
fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE);
// ...
fftw_cleanup_threads();

对大规模 FFT (N > 1<<18) 并行加速显著，小规模开销大于收益。

6. 运行结果展示

频谱图

频谱图
图 1：FFT 频谱幅值，红虚线为 100Hz 截止频率。可见 50Hz 基波(幅值~~480)、150Hz 谐波(幅值~~192)保留，500Hz 高频干扰被抑制。

时域波形对比

波形对比
图 2：红=原始含噪/高频信号，绿=低通滤波后平滑波形。滤波后保留 50Hz 基波与 150Hz 谐波，500Hz 高频分量与高频噪声被有效抑制。

滤波器频率响应

滤波器响应
波形局部放大 (0-0.1s)
图 3：矩形窗理想低通响应。实际应用中建议使用窗函数(汉宁/布莱克曼)平滑过渡带，减少吉布斯现象。

统计数据

原始信号 RMS: 0.7973
滤波后 RMS: 0.7087
RMS 降低: 11.1%

RMS 降低主要来自 500Hz 分量(幅值 0.3)与高频噪声的移除。

7. 常见坑与最佳实践

问题	现象	解决方案
输入长度非 2 的幂	速度慢、精度差	零填充至最近 2 的幂：`while(sz<N) sz*=2;`
内存未对齐	崩溃(SIGSEGV/SIGBUS)	必须用 `fftw_malloc`/`fftw_free`，不可 `new`/`malloc`/`vector`
MEASURE 覆盖数据	频谱全零	先创建 plan，再填充输入数组
忘记归一化	幅值偏大 N 倍	`c2r` 结果必须 `/N`
多线程无加速	CPU 占用不增	调用 `fftw_init_threads(); fftw_plan_with_nthreads(n);`
单精度需求	内存/带宽受限	链接 `fftwf`，API 前缀 `fftwf_`，类型 `fftwf_complex`
部署动态库缺失	运行报错 `libfftw3.so.3 not found`	打包 `.so` 或静态链接 `-Wl,-Bstatic -lfftw3 -Wl,-Bdynamic`
频谱泄漏	峰值展宽、旁瓣高	加窗(汉宁/布莱克曼)、增加 N、零填充

8. 参考资料

comment:

LiveRe
ChangYan