单片机数字信号处理：FIR滤波器与Goertzel算法实战

1. 单片机数字信号处理概述

数字信号处理（DSP）在现代嵌入式系统中扮演着至关重要的角色。从智能手机的语音识别到工业传感器的噪声过滤，DSP算法无处不在。传统上，这类任务由专用DSP芯片或FPGA完成，但现代单片机通过集成MAC（乘积累加）引擎和高性能CPU内核，已经能够胜任许多实时信号处理任务。

C8051F12x和C8051F36x系列单片机就是这类"DSP-enabled MCU"的典型代表。它们具有以下关键特性：

• 100 MIPS的8051兼容CPU内核
• 2周期完成16x16位乘法的硬件MAC引擎
• 片上高精度ADC和DAC
• 仅5x5mm的极小封装尺寸

这些特性使得在8位单片机架构上实现实时数字信号处理成为可能。与专用DSP芯片相比，这种方案具有明显的成本和空间优势，同时保留了足够的CPU带宽来处理其他系统任务。

2. 数字FIR滤波器实现

2.1 FIR滤波器基础原理

FIR（有限脉冲响应）滤波器是数字信号处理中最常用的滤波器类型之一。其输出仅依赖于当前和过去的输入值，数学表达式为：

y(n) = C₀x(n) + C₁x(n-1) + C₂x(n-2) + ... + Cₙx(n-N)

其中：

• y(n)是当前输出
• x(n-k)是当前及过去的输入样本
• Cₖ是滤波器系数，决定滤波器的频率响应特性

FIR滤波器的主要优势包括：

1. 绝对稳定性：由于没有反馈回路，FIR滤波器永远不会振荡
2. 线性相位特性：保持信号中各频率成分的时序关系
3. 对系数量化误差不敏感：适合16位定点实现

2.2 单片机上的优化实现

在C8051F系列单片机上实现FIR滤波器时，我们采用了多项优化技术：

2.2.1 循环缓冲区技术

传统实现需要每次采样后移动所有数据，消耗大量CPU周期。我们使用环形缓冲区管理输入样本：

c复制#define TAP_SIZE 32
int16_t sampleBuffer[TAP_SIZE];
uint8_t bufferIndex = 0;

void FIR_ProcessSample(int16_t newSample) {
    sampleBuffer[bufferIndex] = newSample;
    bufferIndex = (bufferIndex + 1) % TAP_SIZE;
}

这种方法完全消除了数据搬移开销，仅需更新索引指针。

2.2.2 系数对称性利用

大多数FIR滤波器系数具有中心对称特性（C₀=Cₙ, C₁=Cₙ₋₁等）。我们可以利用这一特性减少50%的乘法运算：

assembly复制MOV DPTR, #COEFF_START  ; 系数指针
MOV R0, #BUFFER_START   ; 缓冲区起始
MOV R1, #BUFFER_END     ; 缓冲区末尾

MAC_Loop:
MOVX A, @DPTR          ; 加载系数Ck
INC DPTR
MOVX B, @R0            ; 加载样本x[k]
INC R0
MOVX C, @R1            ; 加载样本x[N-k]
DEC R1
MAC A, B               ; Ck*x[k]
MAC A, C               ; Ck*x[N-k]
DJNZ R7, MAC_Loop      ; 循环计数

2.2.3 半带滤波器优化

对于半带滤波器（每隔一个系数为零的特殊FIR设计），可以跳过零系数乘法，进一步提升效率。

2.3 性能实测数据

我们在C8051F360上测试了不同阶数FIR滤波器的性能：

即使处理20阶FIR滤波器，CPU仍有90%以上的带宽可用于其他任务，充分展示了这种方案的实用性。

3. Goertzel算法与DTMF解码

3.1 Goertzel算法原理

Goertzel算法是DFT的一种高效实现形式，专门用于检测信号中特定频率成分的能量。其核心方程为：

Q₀ = coefₖ × Q₁ - Q₂ + x[n]
Q₁ = Q₀[n-1]
Q₂ = Q₁[n-1]

其中：

• x[n]是输入样本
• coefₖ = 2cos(2πk/N)，k为目标频率对应的DFT bin号
• 最终能量计算：Power = Q₁² + Q₂² - coefₖ × Q₁ × Q₂

与传统FFT相比，Goertzel算法在检测少量频率时效率更高，非常适合DTMF（双音多频）解码等应用。

3.2 DTMF解码实现

DTMF信号由两组频率组合而成：

• 低频组：697Hz, 770Hz, 852Hz, 941Hz
• 高频组：1209Hz, 1336Hz, 1477Hz, 1633Hz

我们的实现方案包含以下关键技术点：

3.2.1 两级检测机制

为防止语音误触发，系统同时检测基频和二次谐波：

1. 首先检测8个基频（4低频+4高频）
2. 仅在检测到有效基频组合后，才启动二次谐波检测
3. 只有基频能量显著大于谐波时，才判定为有效DTMF信号

3.2.2 动态增益控制

ADC输入信号强度可能变化很大，我们在ISR中实现自动增益控制：

c复制#define TARGET_LEVEL 2000
int16_t AGC_Gain = 256;

void ADC_ISR() {
    int32_t sample = ADC_RESULT;
    sample = (sample * AGC_Gain) >> 8;
    
    // 更新增益
    if(abs(sample) > TARGET_LEVEL) {
        AGC_Gain = (AGC_Gain * TARGET_LEVEL) / abs(sample);
    }
}

3.2.3 实时处理优化

通过合理安排计算顺序，所有Goertzel滤波计算均在ADC中断服务例程中完成：

1. 每125μs（8kHz采样率）触发一次ADC中断
2. 在中断中：
   • 读取ADC样本并应用AGC
   • 更新8个Goertzel滤波器状态
   • 当收集足够样本后（通常N=205），计算各频率能量
3. 主循环仅需检查检测结果

3.3 性能对比

MAC引擎使计算速度提升8倍以上，实现了真正的实时DTMF检测。

4. FFT算法实现

4.1 FFT基本原理

FFT通过蝶形运算高效计算DFT。N点FFT的计算复杂度从DFT的O(N²)降低到O(NlogN)。基本蝶形运算：

A' = A + Wₙᵏ × B
B' = A - Wₙᵏ × B

其中Wₙᵏ = e^(-j2πk/N)是旋转因子。

4.2 单片机优化实现

4.2.1 定点数运算

为节省资源，我们使用Q15格式定点数表示旋转因子：

c复制typedef struct {
    int16_t real;
    int16_t imag;
} Complex;

const Complex W32[16] = {
    {32767, 0}, {32138, -6393}, {30273, -12540}, 
    {27245, -18205}, {23170, -23170}, {18205, -27245},
    // ... 其他旋转因子
};

4.2.2 位反转寻址

FFT要求输入数据按位反转顺序排列。我们预先计算位反转索引表：

c复制const uint8_t bitRevTable[64] = {
    0x00, 0x20, 0x10, 0x30, 0x08, 0x28, 0x18, 0x38,
    // ... 其他位反转索引
};

void BitRevReorder(Complex* data) {
    for(uint8_t i=0; i<64; i++) {
        uint8_t j = bitRevTable[i];
        if(i < j) {
            Complex temp = data[i];
            data[i] = data[j];
            data[j] = temp;
        }
    }
}

4.2.3 窗函数应用

为减少频谱泄漏，我们使用汉宁窗预处理数据：

c复制const int16_t hanningWindow[64] = {
    0, 83, 331, 741, 1309, 2027, 2886, 3876,
    // ... 其他窗系数
};

void ApplyWindow(Complex* data) {
    for(uint8_t i=0; i<64; i++) {
        data[i].real = (data[i].real * hanningWindow[i]) >> 15;
        data[i].imag = (data[i].imag * hanningWindow[i]) >> 15;
    }
}

4.3 性能考量

在C8051F360上实现64点FFT的关键指标：

• 需要约12KB代码空间（包含旋转因子表）
• 执行时间约15ms（6.6% CPU占用率@100MHz）
• 动态范围约72dB（16位定点实现）

虽然性能不及专用DSP，但已能满足许多嵌入式应用的频域分析需求。

5. 实际应用建议

5.1 算法选择指南

5.2 资源优化技巧

1. 存储器管理：

   • 将系数表放在CODE区（Flash）
   • 使用xdata关键字将大型缓冲区放在外部RAM
   • 复用缓冲区空间（如FFT的输入/输出可共用同一数组）

2. 中断优先级设置：

   c复制void InitInterrupts() {
       IP = 0x10;    // 将ADC中断设为高优先级
       EADC = 1;     // 使能ADC中断
       EA = 1;       // 全局中断使能
   }
   

3. 低功耗设计：

   • 在采样间隔期间进入IDLE模式
   • 动态调整CPU时钟频率
   • 关闭未使用的外设时钟

5.3 调试与验证方法

1. 时域验证：

   • 通过DAC输出中间信号
   • 使用逻辑分析仪捕获处理流程

2. 频域验证：

   matlab复制% 在MATLAB中对比理论响应与实际输出
   [h,f] = freqz(b,1,512,8000);
   plot(f,20*log10(abs(h))); 
   hold on;
   plot(measuredFreq, measuredGain, 'r--');
   

3. 性能分析：

   • 使用定时器测量关键代码段执行时间
   • 通过UART输出资源使用统计信息

通过合理选择算法并充分利用单片机硬件特性，开发者可以在资源受限的嵌入式系统中实现高效的实时信号处理。C8051F系列的MAC引擎和高速CPU内核为这类应用提供了理想的平台，兼顾了性能、成本和功耗的平衡。