📚 课程目录第4章 CORDIC迭代次数与精度关系分析

4. CORDIC迭代次数与精度关系分析

任务目标

以16位定点数（Q1.14格式）实现CORDIC向量模式，计算输入角度对应的余弦与正弦值。通过改变迭代次数n（1~16），统计与IEEE 754浮点参考值之间的最大误差，绘制误差曲线，确定满足精度要求（误差<1e-3）的最小迭代次数。

实现方案

采用旋转模式（Rotation Mode）CORDIC，角度范围限定在[-π/2, π/2]内（利用三角恒等式将任意角度折叠至该区间）。每次迭代执行：

if (sigma >= 0) { x' = x - y * d; y' = y + x * d; z' = z - atan_i; }
else          { x' = x + y * d; y' = y - x * d; z' = z + atan_i; }

其中d = 2^-i，atan_i为预存的反正切值（Q1.14定点）。迭代完成后乘以增益因子K≈0.60725（定点化为9952/16384）。

关键完整代码

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

#define Q14 16384.0
#define ITER_MAX 16

// 预存atan(2^-i)的Q1.14定点值（弧度*16384）
static const int16_t atan_table[ITER_MAX] = {
    12867, 7596, 4010, 2037, 1023, 512, 256, 128,
    64, 32, 16, 8, 4, 2, 1, 1
};

// 定点CORDIC旋转模式，n为迭代次数
void cordic_rot(int32_t theta_q, uint8_t n, int32_t *cos_q, int32_t *sin_q) {
    int32_t x = 19898;  // 1/K 的Q1.14定点值（约1.64676*16384）
    int32_t y = 0;
    int32_t z = theta_q;
    for (int i = 0; i < n; i++) {
        int32_t dx = y >> i;
        int32_t dy = x >> i;
        if (z >= 0) {
            x -= dx;  y += dy;  z -= atan_table[i];
        } else {
            x += dx;  y -= dy;  z += atan_table[i];
        }
    }
    *cos_q = x;
    *sin_q = y;
}

// 角度折叠到[-pi/2, pi/2]，返回象限标志
float fold_angle(float deg, int *quadrant) {
    float rad = deg * M_PI / 180.0f;
    *quadrant = (int)(rad / (M_PI/2)) % 4;
    if (*quadrant < 0) *quadrant += 4;
    float rem = rad - (*quadrant) * (M_PI/2);
    return rem;
}

int main() {
    float test_angles[] = {0, 15, 30, 45, 60, 75, 90, 120, 150, 180, 225, 270, 315, 360};
    int num_angles = sizeof(test_angles)/sizeof(test_angles[0]);
    
    printf("迭代次数 | 最大误差(cos) | 最大误差(sin) | 是否满足<1e-3\n");
    printf("---------|-------------------|-------------------|-------------\n");
    
    for (int n = 1; n <= ITER_MAX; n++) {
        float max_err_cos = 0, max_err_sin = 0;
        for (int a = 0; a < num_angles; a++) {
            int quad;
            float folded = fold_angle(test_angles[a], &quad);
            int32_t theta_q = (int32_t)(folded * Q14);
            
            int32_t cos_q, sin_q;
            cordic_rot(theta_q, n, &cos_q, &sin_q);
            
            float cos_f = cos_q / Q14;
            float sin_f = sin_q / Q14;
            
            // 根据象限还原真实符号
            float cos_real, sin_real;
            switch (quad) {
                case 0: cos_real = cos_f;  sin_real = sin_f;  break;
                case 1: cos_real = -sin_f; sin_real = cos_f;  break;
                case 2: cos_real = -cos_f; sin_real = -sin_f; break;
                case 3: cos_real = sin_f;  sin_real = -cos_f; break;
            }
            
            float ref_cos = cos(test_angles[a] * M_PI / 180.0f);
            float ref_sin = sin(test_angles[a] * M_PI / 180.0f);
            
            float err_cos = fabs(cos_real - ref_cos);
            float err_sin = fabs(sin_real - ref_sin);
            if (err_cos > max_err_cos) max_err_cos = err_cos;
            if (err_sin > max_err_sin) max_err_sin = err_sin;
        }
        int pass = (max_err_cos < 1e-3 && max_err_sin < 1e-3);
        printf("%4d      | %.6e        | %.6e        | %s\n",
               n, max_err_cos, max_err_sin, pass ? "是" : "否");
    }
    return 0;
}

编译与运行

gcc -o cordic_analysis cordic_analysis.c -lm
./cordic_analysis

验证方法与判定标准

程序遍历14个测试角度（覆盖四个象限及边界），对每个迭代次数n计算所有角度的最大误差。观察输出表格：

误差随n增大呈指数下降趋势，n=1时误差约0.4量级，n≥10后误差降至1e-4以下；
当n=12时，cos与sin的最大误差均小于1e-3，此时输出标记为“是”；
继续增大n至16，误差不再显著下降（受限于16位定点量化精度，极限约1e-4~1e-5）。

通过条件：找到第一个满足“最大误差<1e-3”的迭代次数n_min，并确认n_min≤16。预期n_min=12（若量化误差导致边界波动，n=11或13也属正常，取决于角度集合选择）。

精度分析要点

定点CORDIC误差来源有三：

旋转角度量化误差：θ以Q1.14表示，量化步长约6.1e-5弧度，对应角度误差约0.0035°；
有限迭代截断误差：n次迭代后残余角度z_n≤atan(2^-n)，n=12时约2.4e-4弧度；
定点运算舍入误差：每次移位右移丢失低位，随迭代累积。

三者共同决定最终精度下限。当n>14时，第2项已小于第1、3项，继续增加迭代不再改善精度，此时应通过增加数据位宽（如Q1.30）而非迭代次数来提升精度。

上图示意误差随迭代次数增加而下降的趋势（非实测数据，仅用于展示误差量级与阈值关系）。实际运行程序后，可观察到cos与sin误差曲线均从约0.4量级单调下降，在n=12附近穿越1e-3阈值线。