嵌入式整数信号变换库：纯定点FFT/DCT实现

1. 项目概述“Transform”是一个专为嵌入式系统设计的轻量级整数信号变换库其核心目标是在无浮点硬件支持、内存受限、实时性敏感的MCU平台上以纯整数运算实现高精度、低开销的常用信号变换算法。该库不依赖标准数学库如math.h、不使用float或double类型所有计算均基于int16_t、int32_t及定点Q格式如Q15、Q31完成适用于STM32F0/F1/F3、NXP Kinetis L系列、ESP32non-RTOS裸机模式、RISC-V GD32VF103等资源紧张的主流MCU平台。与通用DSP库如ARM CMSIS-DSP不同“Transform”并非对浮点FFT的简单整数化封装而是从底层算法结构出发针对整数域重新设计了蝶形运算单元、旋转因子生成机制、位宽饱和控制策略及内存访问模式。其设计哲学可概括为三点确定性所有函数执行周期严格可预测无分支预测失败、无动态内存分配、无递归调用可验证性所有中间结果均可通过查表法或手工推演复现便于在安全关键场景如电机控制环路、医疗传感器前端中进行形式化验证可裁剪性通过预编译宏精细控制功能子集最小配置下仅需约1.2KB Flash与128B RAM即可运行128点基2-FFT。该库并非学术研究型工具而是面向量产嵌入式产品的工程实践产物——它解决的是真实产线中反复出现的问题ADC采样数据需在2ms内完成频谱分析以触发过流保护8位MCU需在不增加BOM成本的前提下实现振动频谱识别低功耗SoC在深度睡眠唤醒后300μs内必须完成一次DCT压缩以降低无线传输带宽。2. 核心变换算法实现原理2.1 整数FFT基2-DIT与旋转因子定点化“Transform”采用基2按时间抽取DIT-FFT算法但其旋转因子Twiddle Factor生成与传统实现存在本质差异。标准浮点FFT中$W_N^k e^{-j2\pi k/N}$ 通过cos()/sin()查表或CORDIC计算而本库采用预计算分段线性插值Q31定点缩放三重保障机制所有$N \leq 1024$的旋转因子均在编译期通过Python脚本tools/gen_twiddle.py离线生成以const int32_t twiddle_q31[512]形式存储于Flash对于$N 1024$启用运行时分段线性插值将单位圆划分为64段每段存储首尾两点Q31值插值系数由k % 64直接索引误差0.0015°关键蝶形运算统一采用Q31输入、Q31输出乘加过程使用__SMLALDARM Cortex-M3/M4或mulhRISC-V指令保障32×32→64位中间结果不溢出最终通过__SSAT饱和截断强制回写为Q31。典型蝶形运算代码如下transform_fft.c// Q31蝶形a a W*b; b a - W*b a,b为Q31输入W为Q31旋转因子 static inline void fft_butterfly_q31(int32_t *a, int32_t *b, int32_t w) { int64_t prod (int64_t)(*a) * w; // a*W: Q31*Q31 Q62 int32_t t (int32_t)((prod (1LL 30)) 31); // Q62 - Q31 (四舍五入) int32_t a_new __SSAT(*a t, 32); // 饱和加法 int32_t b_new __SSAT(*a - t, 32); *a a_new; *b b_new; }此实现避免了传统整数FFT中常见的“缩放因子漂移”问题——每次蝶形后不进行除法缩放而是在整个FFT结束时统一右移log2(N)位既消除累积误差又节省12个周期/蝶形。2.2 DCT-II快速递归结构与整数化优化库中DCT-II离散余弦变换采用Leng–Jain快速算法其核心是将$N$点DCT分解为两个$N/2$点DCT与一个$N/2$点FFT。为适配整数域“Transform”对此结构进行三项改造DCT旋转因子整数化$c_k \cos(\pi k / (2N))$ 转换为Q15格式利用对称性仅存储前$N/4$个值FFT子模块复用直接调用库内已验证的整数FFT模块避免重复实现奇偶分离预处理输入序列先执行位反转重排与符号翻转x[n] x[n] * (-1)^n使后续DCT计算完全由实数运算构成消除虚部管理开销。DCT初始化函数签名明确体现其工程约束/** * brief 初始化DCT-II变换器 * param handle DCT句柄指针 * param len 变换长度必须为2的幂且4 ≤ len ≤ 2048 * param scale_mode 缩放模式0无缩放1正交归一化输出乘sqrt(2/N) * param work_buf 工作缓冲区大小2*len*sizeof(int32_t)用于FFT中间存储 * return TRANSFORM_OK 成功TRANSFORM_ERR_INVALID_LEN 长度非法 */ transform_status_t transform_dct_init(transform_dct_handle_t *handle, uint16_t len, uint8_t scale_mode, int32_t *work_buf);scale_mode1时库不计算sqrt()而是查表const uint32_t sqrt2_by_n_q31[11] { /* 2/sqrt(4), 2/sqrt(8), ..., 2/sqrt(2048) */ }确保正交DCT输出能量守恒。2.3 其他变换DST-I与Hilbert变换的整数实现DST-I离散正弦变换通过输入序列奇延拓$(2N)$点FFT实现库内提供专用transform_dst_i_init()函数自动配置延拓缓冲区避免用户手动复制数据Hilbert变换非频域方法采用9阶FIR滤波器逼近理想90°相移系数经Parks-McClellan算法优化并量化为Q15整数脉冲响应存于ROM执行时调用transform_fir_q15()完成卷积。所有变换均遵循同一内存模型输入/输出缓冲区地址由用户传入库内不申请任何动态内存。典型调用流程为#define FFT_LEN 256 int32_t adc_data[FFT_LEN]; // 原始ADC采样Q15 int32_t fft_out[FFT_LEN]; // FFT输出Q31实部/虚部交错 int32_t twiddle_buf[FFT_LEN/2]; // 旋转因子缓冲区仅基2需要 transform_fft_handle_t fft_h; transform_fft_init(fft_h, FFT_LEN, TRANSFORM_FFT_RADIX2, twiddle_buf); // 数据准备ADC值左移16位转Q31 for(uint16_t i 0; i FFT_LEN; i) { fft_out[i] (int32_t)adc_data[i] 16; } // 执行原位FFTfft_out被覆盖为频域结果 transform_fft_execute(fft_h, fft_out); // 幅值计算|X[k]| sqrt(Re² Im²)使用Q31定点sqrt for(uint16_t k 0; k FFT_LEN/2; k) { int32_t re fft_out[2*k]; int32_t im fft_out[2*k1]; mag_spectrum[k] transform_sqrt_q31(re*re im*im); // Q31输入Q31输出 }3. API接口详解3.1 全局状态与错误码库定义统一返回类型transform_status_t强制用户检查关键操作结果枚举值含义典型触发场景TRANSFORM_OK操作成功所有函数正常完成TRANSFORM_ERR_INVALID_LEN长度非法FFT长度非2的幂或超出[4,2048]范围TRANSFORM_ERR_NULL_PTR空指针传入NULL工作缓冲区或句柄TRANSFORM_ERR_NOT_INIT未初始化调用execute前未调用initTRANSFORM_ERR_BUF_SIZE缓冲区不足work_buf长度小于算法要求3.2 FFT模块API函数参数说明关键约束transform_fft_init()handle: 句柄指针len: 变换长度radix:TRANSFORM_FFT_RADIX2或TRANSFORM_FFT_RADIX4twiddle_buf: 旋转因子缓冲区Radix-2需len/2个int32_tRadix-4需len/4个Radix-4仅支持len≥16且len%40twiddle_buf必须4字节对齐transform_fft_execute()handle: 已初始化句柄data: 输入/输出缓冲区原位运算长度len个int32_tdata必须为int32_t数组实部/虚部交错存储data[0]Re0,data[1]Im0, ...transform_fft_scale()data: FFT输出缓冲区len: 长度shift: 右移位数通常log2(len)手动缩放用于避免execute中统一缩放导致的精度损失3.3 DCT/DST模块API函数功能要点transform_dct_init()支持DCT-I/II/III/IV通过type参数选择work_buf需2*len*sizeof(int32_t)容纳FFT中间数据transform_dct_execute()原位运算输入为int16_tQ15输出为int32_tQ31若scale_mode1输出已正交归一化transform_dst_i_init()内部自动分配2*len长度的延拓缓冲区用户只需提供len长度原始数据指针3.4 辅助工具函数函数用途定点格式transform_sqrt_q31()Q31输入的32位整数平方根输入Q31输出Q31即$\sqrt{x} \times 2^{31}$transform_atan2_q31()Q31坐标系下的atan2(y,x)输出Q31弧度$-\pi$ to $\pi$transform_bitrev_table_gen()生成位反转查找表表项为uint16_t用于加速FFT重排4. 工程实践指南4.1 内存与性能优化配置在transform_config.h中用户可通过宏精细控制资源占用// 启用/禁用特定变换注释掉则不编译节省Flash #define TRANSFORM_ENABLE_FFT 1 #define TRANSFORM_ENABLE_DCT 1 #define TRANSFORM_ENABLE_DST 0 // 禁用DST节省~800B Flash // 选择FFT基数影响速度与代码大小 #define TRANSFORM_FFT_DEFAULT_RADIX 2 // Radix-2: 代码小速度中等 // #define TRANSFORM_FFT_DEFAULT_RADIX 4 // Radix-4: 速度快15%代码300B // 启用调试断言仅开发阶段 // #define TRANSFORM_DEBUG_ASSERT 1性能实测数据STM32F407VG 168MHz变换类型长度执行周期Flash占用RAM占用静态FFT25614,2002.1 KB12 B句柄len/2twiddleDCT-II12818,5003.4 KB12 B 2*len*4work_bufHilbert643,8001.7 KB12 B 9*4FIR系数注work_buf为栈/全局变量不计入静态RAMtwiddle_buf可置于Flash只读或RAM可写支持动态长度。4.2 与HAL/FreeRTOS集成示例在FreeRTOS任务中安全使用FFT需注意缓冲区隔离与临界区保护// FreeRTOS任务每10ms执行一次频谱分析 void vSpectrumTask(void *pvParameters) { static int32_t adc_raw[256]; static int32_t fft_result[256]; static int32_t twiddle_rom[128]; // 存于Flash只读 transform_fft_handle_t fft_h; transform_fft_init(fft_h, 256, TRANSFORM_FFT_RADIX2, twiddle_rom); for(;;) { // 1. 从DMA缓冲区拷贝最新256点ADC数据Q15 memcpy(adc_raw, hdma_adc.Instance-NDTR, 256 * sizeof(int16_t)); // 2. 转Q31并执行FFT临界区保护共享缓冲区 taskENTER_CRITICAL(); for(uint16_t i 0; i 256; i) { fft_result[i] (int32_t)adc_raw[i] 16; } transform_fft_execute(fft_h, fft_result); taskEXIT_CRITICAL(); // 3. 计算幅值谱可于非临界区执行 for(uint16_t k 0; k 128; k) { int32_t re fft_result[2*k]; int32_t im fft_result[2*k1]; mag_spectrum[k] transform_sqrt_q31(re*re im*im); } vTaskDelay(10); // 10ms周期 } }4.3 常见问题与解决方案问题FFT输出幅值随长度增大而衰减原因未调用transform_fft_scale()或shift参数错误。方案执行transform_fft_scale(fft_result, 256, 8)256点需右移8位。问题DCT输出含明显直流偏移原因输入数据未去均值。方案在transform_dct_execute()前计算均值并减去int32_t mean 0; for(uint16_t i 0; i 128; i) mean (int32_t)dct_input[i]; mean 7; // /128 for(uint16_t i 0; i 128; i) dct_input[i] - (int16_t)mean;问题编译报错“undefined reference to__aeabi_idiv”原因链接器未包含整数除法库。方案在GCC链接选项中添加-u _printf_float禁用浮点printf并确保-lc已链接或定义#define TRANSFORM_NO_DIVISION 1启用查表除法。5. 实际应用案例电机振动频谱监测某BLDC电机控制器需在STM32G071上实现振动异常检测。方案如下硬件ADXL345加速度计I2C±2g12-bit采样率1kHz软件每秒截取4段256点窗口重叠率50%对每段执行FFT提取0-200Hz频带对应0-25.5 bin能量关键代码// 预计算200Hz对应bin (200 * 256) / 1000 51.2 → 取51 const uint16_t FREQ_BIN_MAX 51; // 计算0-200Hz能量Q31累加 int64_t energy_0_200 0; for(uint16_t k 0; k FREQ_BIN_MAX; k) { int32_t mag transform_sqrt_q31( fft_result[2*k] * fft_result[2*k] fft_result[2*k1] * fft_result[2*k1] ); energy_0_200 mag; } energy_0_200 31; // Q31 → 纯整数 // 触发阈值energy_0_200 15000标定值 if(energy_0_200 15000) { HAL_GPIO_WritePin(ALERT_GPIO_Port, ALERT_Pin, GPIO_PIN_SET); }此方案在G07164MHz下单次FFT能量计算耗时800μs满足实时性要求且因全程整数运算避免了浮点单元不可用G0系列无FPU导致的性能崩溃。6. 源码结构与构建说明库目录结构严格遵循嵌入式项目惯例transform/ ├── Inc/ # 公共头文件 │ ├── transform.h # 主接口声明 │ ├── transform_fft.h # FFT专用接口 │ └── transform_config.h # 用户可配置宏 ├── Src/ # 核心源码 │ ├── transform.c # 全局初始化与错误处理 │ ├── transform_fft.c # FFT实现含蝶形、位反转 │ ├── transform_dct.c # DCT/DST实现 │ └── transform_math.c # sqrt/atan2等定点数学函数 ├── Tools/ # 辅助工具 │ └── gen_twiddle.py # 旋转因子生成脚本Python 3.6 └── Examples/ # STM32CubeIDE工程模板 └── STM32G071RB_Transform_FFT构建要求编译器ARM GCC 10.2 或 IAR EWARM 8.50必须启用-O2或-O3优化否则定点乘法性能下降40%transform_config.h需置于用户工程Inc/路径优先于库内默认配置。最后强调一个硬性工程纪律所有变换函数的输入数据必须经过ADC校准与零点补偿否则整数运算的量化误差会被指数级放大。我们曾在一个风电变桨控制器项目中发现未补偿的ADC零点偏移±3LSB导致FFT基波幅值测量误差达12%远超整数FFT本身的理论误差0.5%。这提醒我们再精妙的算法也必须扎根于扎实的硬件前端设计。