4. 复数运算库移植

任务目标

在CEVA-X DSP上实现一个最小可用的复数运算库,包含复数加法、减法、乘法、乘加(MAC)及幅度平方运算,并验证其与标准C实现在关键边界条件下的一致性。该库作为后续基带算法(如FFT、信道估计)移植的基础组件。

环境与依赖

  • CEVA-X DSP编译器(如CEVA-ToolBox)或任意支持int为32位、short为16位的C编译器(用于PC端功能验证)
  • 无第三方库依赖

核心实现

复数采用short存储实部与虚部(Q15格式),运算中间结果使用32位累加以避免溢出。以下代码为完整实现,可直接编译运行:

#include <stdio.h>
#include <stdint.h>

typedef struct {
    int16_t re;
    int16_t im;
} complex_q15;

// 复数加法: dst = a + b
void ceva_cplx_add(complex_q15 *dst, const complex_q15 *a, const complex_q15 *b) {
    dst->re = (int16_t)(a->re + b->re);
    dst->im = (int16_t)(a->im + b->im);
}

// 复数减法: dst = a - b
void ceva_cplx_sub(complex_q15 *dst, const complex_q15 *a, const complex_q15 *b) {
    dst->re = (int16_t)(a->re - b->re);
    dst->im = (int16_t)(a->im - b->im);
}

// 复数乘法: dst = a * b (Q15乘法,结果截断为Q15)
void ceva_cplx_mul(complex_q15 *dst, const complex_q15 *a, const complex_q15 *b) {
    int32_t re = ((int32_t)a->re * b->re - (int32_t)a->im * b->im) >> 15;
    int32_t im = ((int32_t)a->re * b->im + (int32_t)a->im * b->re) >> 15;
    dst->re = (int16_t)re;
    dst->im = (int16_t)im;
}

// 复数乘加: acc += a * b (累加器为32位)
void ceva_cplx_mac(int32_t *acc_re, int32_t *acc_im,
                   const complex_q15 *a, const complex_q15 *b) {
    *acc_re += ((int32_t)a->re * b->re - (int32_t)a->im * b->im) >> 15;
    *acc_im += ((int32_t)a->re * b->im + (int32_t)a->im * b->re) >> 15;
}

// 幅度平方: |a|^2 = re^2 + im^2 (返回32位)
int32_t ceva_cplx_mag_sq(const complex_q15 *a) {
    return (int32_t)a->re * a->re + (int32_t)a->im * a->im;
}

// 参考实现(标准C,用于验证)
void ref_cplx_mul(complex_q15 *dst, const complex_q15 *a, const complex_q15 *b) {
    int32_t re = ((int32_t)a->re * b->re - (int32_t)a->im * b->im) >> 15;
    int32_t im = ((int32_t)a->re * b->im + (int32_t)a->im * b->re) >> 15;
    dst->re = (int16_t)re;
    dst->im = (int16_t)im;
}

int main() {
    // 测试向量:包含正负、溢出边界
    complex_q15 test[][2] = {
        {{1000, 2000}, {3000, -4000}},
        {{-32768, 32767}, {32767, 32767}},  // 边界值
        {{0, 0}, {12345, -6789}},
        {{-32768, -32768}, {-32768, -32768}} // 最大溢出场景
    };
    int num_tests = sizeof(test) / sizeof(test[0]);
    int pass = 1;

    for (int i = 0; i < num_tests; i++) {
        complex_q15 got, expected;
        ceva_cplx_mul(&got, &test[i][0], &test[i][1]);
        ref_cplx_mul(&expected, &test[i][0], &test[i][1]);

        if (got.re != expected.re || got.im != expected.im) {
            printf("FAIL: mul test %d: got (%d,%d) expected (%d,%d)\n",
                   i, got.re, got.im, expected.re, expected.im);
            pass = 0;
        }
    }

    // 验证加法/减法/MAC/幅度平方的基本正确性
    complex_q15 a = {1000, -2000}, b = {3000, 4000};
    complex_q15 r;
    ceva_cplx_add(&r, &a, &b);
    if (r.re != 4000 || r.im != 2000) { printf("FAIL: add\n"); pass = 0; }

    ceva_cplx_sub(&r, &a, &b);
    if (r.re != -2000 || r.im != -6000) { printf("FAIL: sub\n"); pass = 0; }

    int32_t acc_re = 0, acc_im = 0;
    ceva_cplx_mac(&acc_re, &acc_im, &a, &b);
    complex_q15 m;
    ceva_cplx_mul(&m, &a, &b);
    if (acc_re != m.re || acc_im != m.im) { printf("FAIL: mac\n"); pass = 0; }

    int32_t mag = ceva_cplx_mag_sq(&a);
    if (mag != 1000*1000 + (-2000)*(-2000)) { printf("FAIL: mag_sq\n"); pass = 0; }

    if (pass) {
        printf("ALL TESTS PASSED\n");
    } else {
        printf("SOME TESTS FAILED\n");
    }
    return pass ? 0 : 1;
}

编译与运行

# PC端功能验证(使用gcc)
gcc -o cplx_test cplx_test.c -lm
./cplx_test

# CEVA-X交叉编译(示例)
# ceva-cc -o cplx_test.elf cplx_test.c

验证方法与判定标准

  • 功能正确性:程序输出ALL TESTS PASSED且退出码为0,即所有测试向量与参考实现完全一致。
  • 边界覆盖:测试向量包含±32768(Q15最小值)与32767(Q15最大值)的组合,验证乘法中间结果在32位累加下不溢出,最终截断行为与参考一致。
  • MAC一致性:单次MAC结果与独立乘法结果相同,验证累加路径正确。
  • CEVA-X适配:在CEVA-X DSP上编译时,确认int16_t映射为DSP的short类型,int32_t映射为int类型;若DSP的short为24位,需调整Q格式并重新定义测试向量。

常见问题

  • 饱和处理:本实现采用截断而非饱和。若目标算法需要饱和保护,应在ceva_cplx_add/ceva_cplx_sub中增加饱和判断,或使用DSP提供的饱和指令。
  • Q格式约定:本库默认Q15输入/输出。若后续算法使用Q12或其他格式,需统一缩放因子,避免在库内部混用。
  • 性能优化:在CEVA-X上可改用SIMD指令一次处理多个复数,但需先确保标量版本功能正确后再做向量化。