基于Kiss FFT的嵌入式设备FFT算法实现

引言

快速傅里叶变换（Fast Fourier Transform, FFT）是数字信号处理领域最核心的算法之一，它将时域信号转换为频域表示，广泛应用于音频处理、频谱分析、通信系统等领域。

在资源受限的嵌入式设备上实现FFT面临着诸多挑战：

• 内存限制：嵌入式MCU通常只有几KB到几十KB的RAM
• 计算能力有限：没有浮点单元（FPU）或主频较低
• 实时性要求：需要在有限的时间内完成计算

本文将介绍如何使用轻量级的开源FFT库Kiss FFT在嵌入式设备上高效实现FFT算法，并提供完整的代码示例和优化建议。

为什么选择Kiss FFT

在众多FFT实现中，Kiss FFT特别适合嵌入式应用，主要原因包括：

核心优势

1. 极简设计：整个库只有几个文件，代码量少，易于理解和移植
2. 内存友好：支持固定点数运算，无需动态内存分配
3. 可配置精度：支持float、double、定点数等多种数据类型
4. MIT许可证：商业友好，无版权顾虑
5. 跨平台：纯C语言实现，可在任何支持C编译器的平台上运行

与其他FFT库对比

特性	Kiss FFT	FFTW	CMSIS-DSP
代码量	~500行	庞大	中等
内存占用	极低	高	中等
定点支持	✅	❌	✅
易用性	简单	复杂	中等
性能优化	一般	极致	针对ARM优化

对于资源受限的嵌入式项目，Kiss FFT是一个理想的平衡选择。

Kiss FFT的核心原理

从DFT到FFT

离散傅里叶变换（DFT）

DFT是将时域离散信号转换为频域表示的数学工具。对于N点序列$x[n]$，其DFT定义为：

$$
X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j2\pi kn/N}, \quad k = 0, 1, …, N-1
$$

其中：

• $x[n]$：时域采样值
• $X[k]$：频域系数（复数）
• $e^{-j2\pi kn/N}$：旋转因子（也称为twiddle factor），记为$W_N^{kn}$

直接计算DFT需要N²次复数乘法和N(N-1)次复数加法，时间复杂度为O(N²)。

FFT的核心思想：分治法

FFT通过分治策略将大问题分解为小问题。Cooley-Tukey算法的关键洞察是：

如果N是2的幂，可以将N点DFT分解为两个N/2点DFT。

Radix-2 DIT算法详解

第一步：按奇偶分解

将输入序列$x[n]$按索引的奇偶性分成两组：

• 偶数索引：$x[0], x[2], x[4], …, x[N-2]$
• 奇数索引：$x[1], x[3], x[5], …, x[N-1]$

定义两个新序列：

$$
\begin{gathered} x_e[m] = x[2m] \quad \text{(偶数部分)} \cr x_o[m] = x[2m+1] \quad \text{(奇数部分)} \end{gathered}
$$

其中$m = 0, 1, …, N/2-1$。

第二步：重写DFT公式

将DFT公式按奇偶分解：

$$
\begin{gathered} X[k] = \sum_{n=0}^{N-1} x[n] W_N^{kn} \cr = \sum_{m=0}^{N/2-1} x[2m] W_N^{k(2m)} + \sum_{m=0}^{N/2-1} x[2m+1] W_N^{k(2m+1)} \cr = \sum_{m=0}^{N/2-1} x_e[m] W_{N/2}^{km} + W_N^k \sum_{m=0}^{N/2-1} x_o[m] W_{N/2}^{km} \end{gathered}
$$

注意到$W_N^{2km} = e^{-j2\pi (2km)/N} = e^{-j2\pi km/(N/2)} = W_{N/2}^{km}$。

第三步：定义子DFT

令：

$$
\begin{gathered} E[k] = \sum_{m=0}^{N/2-1} x_e[m] W_{N/2}^{km} \quad \text{(偶数部分的N/2点DFT)} \cr O[k] = \sum_{m=0}^{N/2-1} x_o[m] W_{N/2}^{km} \quad \text{(奇数部分的N/2点DFT)} \end{gathered}
$$

则原DFT可以表示为：

$$
X[k] = E[k] + W_N^k \cdot O[k]
$$

第四步：利用周期性

由于$E[k]$和$O[k]$都是N/2点DFT，它们具有周期性：

$$
E[k + N/2] = E[k], \quad O[k + N/2] = O[k]
$$

同时，旋转因子满足：

$$
W_N^{k+N/2} = W_N^k \cdot W_N^{N/2} = W_N^k \cdot e^{-j\pi} = -W_N^k
$$

因此，对于$k = 0, 1, …, N/2-1$：

$$
\boxed{\begin{gathered} X[k] = E[k] + W_N^k \cdot O[k] \cr X[k + N/2] = E[k] - W_N^k \cdot O[k] \end{gathered}}
$$

这就是著名的蝶形运算公式！

蝶形运算（Butterfly Operation）

蝶形运算的结构

蝶形运算是FFT的基本计算单元，它结合了两个输入产生两个输出：

输入: E[k] 和 O[k]
      ↓
   ┌───────┐
   │  ×W   │ ← 旋转因子乘法
   └───────┘
      ↓
   ╱     ╲
  ╱       ╲
 +         -
  ╲       ╱
   ╲     ╱
    ↓   ↓
输出: X[k]  X[k+N/2]

数学表达

对于每一对$(E[k], O[k])$，蝶形运算执行：

$$
\begin{gathered} X[k] = E[k] + W_N^k \cdot O[k] \cr X[k + N/2] = E[k] - W_N^k \cdot O[k] \end{gathered}
$$

其中$W_N^k = e^{-j2\pi k/N} = \cos(2\pi k/N) - j\sin(2\pi k/N)$是旋转因子。

复数运算细节

设$W_N^k = a - jb$（其中$a = \cos(2\pi k/N)$，$b = \sin(2\pi k/N)$），$O[k] = c + jd$，则：

$$
\begin{gathered} W_N^k \cdot O[k] = (a - jb)(c + jd) \cr = (ac + bd) + j(ad - bc) \end{gathered}
$$

因此：

$$
\begin{gathered} X[k] = E[k] + (ac + bd) + j(ad - bc) \cr X[k + N/2] = E[k] - (ac + bd) - j(ad - bc) \end{gathered}
$$

每次蝶形运算需要：

• 1次复数乘法（4次实数乘法 + 2次实数加法）
• 2次复数加法（4次实数加法）
• 总计：4次实数乘法 + 6次实数加法

递归分解过程

8点FFT示例

让我们通过一个8点FFT的例子来理解整个分解过程：

第0层（原始输入）:        x[0] x[1] x[2] x[3] x[4] x[5] x[6] x[7]
                            ↓   ↓   ↓   ↓   ↓   ↓   ↓   ↓
第1层（位反转排序）:      x[0] x[4] x[2] x[6] x[1] x[5] x[3] x[7]
                          ┌───┴───┐ ┌───┴───┐ ┌───┴───┐ ┌───┴───┐
第2层（4个2点DFT）:       │ 蝶形  │ │ 蝶形  │ │ 蝶形  │ │ 蝶形  │
                          └───┬───┘ └───┬───┘ └───┬───┘ └───┬───┘
                          ┌─────┴─────┐         ┌─────┴─────┐
第3层（2个4点DFT）:       │  蝶形组   │         │  蝶形组   │
                          └─────┬─────┘         └─────┬─────┘
                          ┌───────────┴───────────┐
第4层（1个8点DFT）:       │      蝶形组           │
                          └───────────┬───────────┘
                            ↓   ↓   ↓   ↓   ↓   ↓   ↓   ↓
输出:                      X[0] X[1] X[2] X[3] X[4] X[5] X[6] X[7]

分解步骤

第1步：8点DFT → 2个4点DFT

• 偶数部分：$x[0], x[2], x[4], x[6]$ → 4点DFT得到$E_0, E_1, E_2, E_3$
• 奇数部分：$x[1], x[3], x[5], x[7]$ → 4点DFT得到$O_0, O_1, O_2, O_3$
• 蝶形组合：
  • $X[0] = E_0 + W_8^0 O_0$, $X[4] = E_0 - W_8^0 O_0$
  • $X[1] = E_1 + W_8^1 O_1$, $X[5] = E_1 - W_8^1 O_1$
  • $X[2] = E_2 + W_8^2 O_2$, $X[6] = E_2 - W_8^2 O_2$
  • $X[3] = E_3 + W_8^3 O_3$, $X[7] = E_3 - W_8^3 O_3$

第2步：4点DFT → 2个2点DFT

每个4点DFT继续分解为两个2点DFT，依此类推。

第3步：2点DFT（基情况）

2点DFT是最简单的情况：

$$
\begin{gathered} X[0] = x[0] + x[1] \cr X[1] = x[0] - x[1] \end{gathered}
$$

计算复杂度分析

对于N点FFT（N = 2^M）：

• 层数：$\log_2 N$层
• 每层蝶形运算数：$N/2$个
• 总蝶形运算数：$(N/2) \cdot \log_2 N$
• 总复数乘法：$(N/2) \cdot \log_2 N$次
• 总复数加法：$N \cdot \log_2 N$次

与直接DFT对比（N=1024）：

• DFT：1024² = 1,048,576次复数乘法
• FFT：(1024/2) × 10 = 5,120次复数乘法
• 加速比：约200倍！

位反转排序（Bit-Reversal Permutation）

为什么需要位反转？

在递归分解过程中，输入数据的顺序会被打乱。观察8点FFT的输入顺序：

原始索引:  0   1   2   3   4   5   6   7
二进制:    000 001 010 011 100 101 110 111
位反转:    000 100 010 110 001 101 011 111
十进制:    0   4   2   6   1   5   3   7

FFT算法要求输入数据按照位反转顺序排列，这样蝶形运算才能正确进行。

位反转算法

// 位反转函数
uint16_t bit_reverse(uint16_t x, int bits) {
    uint16_t result = 0;
    for (int i = 0; i < bits; i++) {
        result = (result << 1) | (x & 1);
        x >>= 1;
    }
    return result;
}

// 对数组进行位反转重排
void bit_reverse_permutation(kiss_fft_cpx *data, int n) {
    int bits = log2(n);
    for (int i = 0; i < n; i++) {
        int j = bit_reverse(i, bits);
        if (i < j) {
            // 交换data[i]和data[j]
            kiss_fft_cpx temp = data[i];
            data[i] = data[j];
            data[j] = temp;
        }
    }
}

旋转因子表（Twiddle Factors）

预计算旋转因子

为了避免重复计算三角函数，Kiss FFT预先计算并存储所有需要的旋转因子：

// 预计算旋转因子表
for (int i = 0; i < nfft/2; i++) {
    double angle = -2.0 * M_PI * i / nfft;
    twiddles[i].r = cos(angle);
    twiddles[i].i = sin(angle);
}

旋转因子的对称性

旋转因子具有以下对称性质，可以减少存储：

$$
\begin{gathered} W_N^{k+N/2} = -W_N^k \cr W_N^{N-k} = W_N^{-k} = (W_N^k)^* \end{gathered}
$$

其中$*$表示复共轭。Kiss FFT利用这些性质优化存储和计算。

关键数据结构

// Kiss FFT的核心配置结构
typedef struct {
    int nfft;              // FFT点数（必须是2的幂）
    int inverse;           // 是否为逆变换
    kf_cpx_t *twiddles;    // 旋转因子表
} kiss_fft_cfg;

// 复数数据类型
typedef struct {
    kf_scalar r;  // 实部
    kf_scalar i;  // 虚部
} kf_cpx_t;

在嵌入式项目中集成Kiss FFT

步骤1：获取源代码

从GitHub获取Kiss FFT源码：

git clone https://github.com/mborgerding/kissfft.git

或者直接将以下文件复制到您的项目中：

• kiss_fft.h
• kiss_fft.c
• _kiss_fft_guts.h

步骤2：配置数据类型

根据您的需求在 kiss_fft.h中配置数据类型：

// 选项1：使用浮点数（如果有FPU）
#define KISS_FFT_USE_FLOAT 1

// 选项2：使用定点数（推荐用于无FPU的MCU）
#define KISS_FFT_FIXED_POINT 1
#define KISS_FFT_SCALAR int16_t  // 16位定点数

// 选项3：使用双精度（高精度需求）
#define KISS_FFT_USE_DOUBLE 1

步骤3：基本使用示例

正向FFT（时域→频域）

#include "kiss_fft.h"

#define FFT_SIZE 256  // FFT点数，必须是2的幂

void perform_fft(void)
{
    // 1. 准备输入数据（时域采样）
    kiss_fft_cpx input[FFT_SIZE];
    kiss_fft_cpx output[FFT_SIZE];

    // 初始化输入数据（示例：正弦波）
    for (int i = 0; i < FFT_SIZE; i++) {
        float sample = sinf(2.0f * M_PI * 50.0f * i / 1000.0f);  // 50Hz正弦波
        input[i].r = sample;
        input[i].i = 0;
    }

    // 2. 创建FFT配置
    kiss_fft_cfg cfg = kiss_fft_alloc(FFT_SIZE, 0, NULL, NULL);

    // 3. 执行FFT
    kiss_fft(cfg, input, output);

    // 4. 处理结果（计算幅度谱）
    float magnitude[FFT_SIZE / 2];
    for (int i = 0; i < FFT_SIZE / 2; i++) {
        magnitude[i] = sqrtf(output[i].r * output[i].r +
                            output[i].i * output[i].i);
    }

    // 5. 释放配置
    free(cfg);
}

逆向FFT（频域→时域）

void perform_ifft(void)
{
    kiss_fft_cpx freq_data[FFT_SIZE];
    kiss_fft_cpx time_data[FFT_SIZE];

    // ... 填充频域数据 ...

    // 创建逆FFT配置（第二个参数为1表示逆变换）
    kiss_fft_cfg cfg = kiss_fft_alloc(FFT_SIZE, 1, NULL, NULL);

    // 执行逆FFT
    kiss_fft(cfg, freq_data, time_data);

    // 归一化（IFFT需要除以N）
    for (int i = 0; i < FFT_SIZE; i++) {
        time_data[i].r /= FFT_SIZE;
        time_data[i].i /= FFT_SIZE;
    }

    free(cfg);
}

步骤4：嵌入式优化技巧

优化1：预分配内存（避免动态分配）

// 在栈上或全局区预分配内存
static uint8_t fft_work_buffer[4096];  // 工作缓冲区
static kiss_fft_cpx input_buf[FFT_SIZE];
static kiss_fft_cpx output_buf[FFT_SIZE];

void optimized_fft(void)
{
    // 使用预分配的缓冲区
    kiss_fft_cfg cfg = kiss_fft_alloc(FFT_SIZE, 0,
                                      fft_work_buffer,
                                      &sizeof(fft_work_buffer));

    kiss_fft(cfg, input_buf, output_buf);

    // 不需要free，因为使用的是静态缓冲区
}

优化2：使用定点数运算

// 在kiss_fft.h中启用定点数
#define KISS_FFT_FIXED_POINT 1
#define KISS_FFT_SCALAR int16_t

// 使用时注意定标
#define Q_FORMAT 15  // Q15格式

void fixed_point_fft(void)
{
    kiss_fft_cpx input[FFT_SIZE];
    kiss_fft_cpx output[FFT_SIZE];

    // 将浮点采样转换为Q15定点数
    for (int i = 0; i < FFT_SIZE; i++) {
        float sample = read_adc();  // 读取ADC采样
        input[i].r = (int16_t)(sample * (1 << Q_FORMAT));
        input[i].i = 0;
    }

    kiss_fft_cfg cfg = kiss_fft_alloc(FFT_SIZE, 0, NULL, NULL);
    kiss_fft(cfg, input, output);

    // 处理定点结果
    for (int i = 0; i < FFT_SIZE / 2; i++) {
        int32_t magnitude_sq = (int32_t)output[i].r * output[i].r +
                               (int32_t)output[i].i * output[i].i;
        // 使用整数平方根或查表法
    }

    free(cfg);
}

优化3：重用FFT配置

// 全局配置，只创建一次
static kiss_fft_cfg fft_cfg = NULL;

void init_fft_system(void)
{
    if (fft_cfg == NULL) {
        fft_cfg = kiss_fft_alloc(FFT_SIZE, 0, NULL, NULL);
    }
}

void process_audio_frame(const float *input, float *magnitude)
{
    kiss_fft_cpx in[FFT_SIZE];
    kiss_fft_cpx out[FFT_SIZE];

    // 转换输入
    for (int i = 0; i < FFT_SIZE; i++) {
        in[i].r = input[i];
        in[i].i = 0;
    }

    // 使用重用的配置
    kiss_fft(fft_cfg, in, out);

    // 计算幅度
    for (int i = 0; i < FFT_SIZE / 2; i++) {
        magnitude[i] = sqrtf(out[i].r * out[i].r + out[i].i * out[i].i);
    }
}

void cleanup_fft_system(void)
{
    if (fft_cfg) {
        free(fft_cfg);
        fft_cfg = NULL;
    }
}

实际应用案例：音频频谱分析

下面是一个完整的音频频谱分析仪实现：

#include "kiss_fft.h"
#include <math.h>

#define SAMPLE_RATE 8000      // 采样率8kHz
#define FFT_SIZE 256          // FFT点数
#define BIN_COUNT (FFT_SIZE / 2)  // 频域bin数量
#define FREQ_RESOLUTION (SAMPLE_RATE / (float)FFT_SIZE)  // 频率分辨率

typedef struct {
    kiss_fft_cfg cfg;
    kiss_fft_cpx input[FFT_SIZE];
    kiss_fft_cpx output[FFT_SIZE];
    float magnitude[BIN_COUNT];
    float frequency[BIN_COUNT];
} SpectrumAnalyzer;

// 初始化频谱分析仪
SpectrumAnalyzer* spectrum_analyzer_create(void)
{
    SpectrumAnalyzer *analyzer = malloc(sizeof(SpectrumAnalyzer));
    if (!analyzer) return NULL;

    analyzer->cfg = kiss_fft_alloc(FFT_SIZE, 0, NULL, NULL);
    if (!analyzer->cfg) {
        free(analyzer);
        return NULL;
    }

    // 预计算每个bin对应的频率
    for (int i = 0; i < BIN_COUNT; i++) {
        analyzer->frequency[i] = i * FREQ_RESOLUTION;
    }

    return analyzer;
}

// 处理一帧音频数据
void spectrum_analyzer_process(SpectrumAnalyzer *analyzer,
                               const float *audio_frame)
{
    // 1. 应用窗函数（Hanning窗）
    for (int i = 0; i < FFT_SIZE; i++) {
        float window = 0.5f * (1.0f - cosf(2.0f * M_PI * i / (FFT_SIZE - 1)));
        analyzer->input[i].r = audio_frame[i] * window;
        analyzer->input[i].i = 0;
    }

    // 2. 执行FFT
    kiss_fft(analyzer->cfg, analyzer->input, analyzer->output);

    // 3. 计算幅度谱（dB）
    for (int i = 0; i < BIN_COUNT; i++) {
        float mag = sqrtf(analyzer->output[i].r * analyzer->output[i].r +
                         analyzer->output[i].i * analyzer->output[i].i);
        // 转换为dB，避免log(0)
        analyzer->magnitude[i] = 20.0f * log10f(fmaxf(mag, 1e-6f));
    }
}

// 获取指定频率范围的平均能量
float spectrum_analyzer_get_band_energy(SpectrumAnalyzer *analyzer,
                                        float low_freq, float high_freq)
{
    int low_bin = (int)(low_freq / FREQ_RESOLUTION);
    int high_bin = (int)(high_freq / FREQ_RESOLUTION);

    low_bin = fmax(0, fmin(low_bin, BIN_COUNT - 1));
    high_bin = fmax(low_bin, fmin(high_bin, BIN_COUNT - 1));

    float energy = 0;
    for (int i = low_bin; i <= high_bin; i++) {
        energy += analyzer->magnitude[i];
    }

    return energy / (high_bin - low_bin + 1);
}

// 销毁频谱分析仪
void spectrum_analyzer_destroy(SpectrumAnalyzer *analyzer)
{
    if (analyzer) {
        if (analyzer->cfg) {
            free(analyzer->cfg);
        }
        free(analyzer);
    }
}

// 使用示例
void audio_processing_example(void)
{
    SpectrumAnalyzer *analyzer = spectrum_analyzer_create();
    if (!analyzer) return;

    float audio_buffer[FFT_SIZE];

    while (1) {
        // 从ADC或I2S获取音频数据
        capture_audio(audio_buffer, FFT_SIZE);

        // 处理频谱
        spectrum_analyzer_process(analyzer, audio_buffer);

        // 分析不同频段
        float bass = spectrum_analyzer_get_band_energy(analyzer, 20, 250);
        float mid = spectrum_analyzer_get_band_energy(analyzer, 250, 4000);
        float treble = spectrum_analyzer_get_band_energy(analyzer, 4000, 4000);

        // 根据频谱进行LED显示或其他处理
        update_led_display(bass, mid, treble);
    }

    spectrum_analyzer_destroy(analyzer);
}

常见问题与解决方案

Q1: FFT结果不正确或出现噪声

原因：

• 未应用窗函数导致频谱泄漏
• 输入数据未正确归一化

解决：

// 始终应用窗函数
for (int i = 0; i < FFT_SIZE; i++) {
    float window = 0.5f * (1.0f - cosf(2.0f * M_PI * i / (FFT_SIZE - 1)));
    input[i].r *= window;
}

Q2: 内存不足

解决：

• 减小FFT_SIZE
• 使用定点数而非浮点数
• 重用配置和缓冲区

Q3: 计算速度太慢

解决：

• 启用编译器优化（-O2或-O3）
• 使用硬件FPU
• 考虑使用CMSIS-DSP库（针对ARM Cortex-M优化）
• 降低采样率或FFT大小

总结

Kiss FFT为嵌入式设备提供了一个轻量、高效的FFT实现方案。通过合理配置和优化，可以在资源受限的MCU上实现实时的频谱分析、音频处理等应用。

关键要点：

1. ✅ Kiss FFT代码简洁，易于移植和理解
2. ✅ 支持定点数运算，适合无FPU的设备
3. ✅ 预分配内存避免动态分配的开销
4. ✅ 应用窗函数减少频谱泄漏
5. ✅ 重用FFT配置提高性能

在实际项目中，可以根据具体需求选择合适的FFT库：

• 超资源受限：Kiss FFT
• ARM Cortex-M系列：CMSIS-DSP
• 高性能需求：FFTW（桌面/服务器）

希望本文能帮助您在嵌入式项目中成功实现FFT算法！

参考文献

1. Kiss FFT官方仓库: https://github.com/mborgerding/kissfft
2. Cooley, J. W., & Tukey, J. W. (1965). An algorithm for the machine calculation of complex Fourier series
3. ARM CMSIS-DSP文档: https://arm-software.github.io/CMSIS_5/DSP/html/index.html