C语言卷积积计算异常：点数增多结果劣化问题排查

正弦波生成代码 #

我用以下C代码生成正弦波信号：

typedef struct {
  double amplitude;
  double offset;
  double frequency;        // Hz
  double phase;            // radians
} sin_wave_config_t;

typedef struct {
  size_t n_points;            // 点数
  double fs;                  // 采样频率 (Hz)
  double *data;
} signal_t;

int generate_sin_wave(signal_t *signal, sin_wave_config_t *wave_config,
    double x_start, double x_end, size_t n_points) {
  assert((signal != NULL) && (wave_config != NULL));
  assert(x_end > x_start);
  assert(n_points > 1);

  if (signal_set_size(signal, n_points) != EXIT_SUCCESS) {
    return EXIT_FAILURE;
  }

  double step = (x_end - x_start) / (n_points - 1);
  signal->fs = 1 / step;

  double x = x_start;
  for (size_t i = 0; i < n_points; i++) {
    signal->data[i] = (wave_config->amplitude / 2)
        * sin(2 * M_PI * wave_config->frequency * x + wave_config->phase)
        + wave_config->offset;
    x += step;
  }
  return EXIT_SUCCESS;
}

初始卷积实现 #

用于计算input_signal1 * input_signal2并将结果存入output_signal_copy的卷积算法：

size_t p = input_signal_1->n_points;
size_t q = input_signal_2->n_points;
output_signal_copy->n_points = p+q-1;
size_t r, s;
float t;
for(size_t n = 0; n < output_signal_copy->n_points; n++){
    t = 0;
    if(n+1 > q){
        r = n+1-q;
    }
    else{
        r = 0;
    }

    s = __min__(n,p-1);
    for(size_t i = r; i <= s; i++){
        t += input_signal_1->data[i] * input_signal_2->data[n-i];
    }
    output_signal_copy->data[n] = t;
}

问题现象 #

输入信号为1kHz正弦波（振幅A）与25kHz正弦波（振幅A/2）的混合信号，与29点低通滤波器冲激响应卷积：

• 500点输入时结果正常，输出信号符合低通滤波预期
• 增加到5000点时，输出结果劣化，输入信号本身也出现异常
• 输出与输入信号的差值呈现周期性模式，振幅约为输入信号的7%

简化卷积实现后问题依旧 #

我简化了卷积算法，但问题仍然存在：

signal_set_to_zero(output_signal_copy); // 将所有数据置0

size_t p = input_signal_1->n_points;
size_t q = input_signal_2->n_points;
for(size_t i = 0; i < p; i++){

    for(size_t j = 0; j < q; j++){
        output_signal_copy->data[i+j] += input_signal_1->data[i]*input_signal_2->data[j];
    }

}

可能的问题原因 #

• 浮点数精度不匹配：输入信号的data是double类型，但卷积累加变量t使用了精度更低的float。高点数下累加次数剧增，精度损失被不断放大，最终导致输出误差。
• 正弦波生成的累积误差：生成正弦波时，x通过x += step逐步累加，当点数极大时，step的浮点误差会累积，导致x的实际值与理论值偏差越来越大，正弦波相位偏移，混合信号频谱畸变，进而影响卷积结果。
• 内存分配隐患：5000点输入时，输出信号需要分配5000+29-1=5028个double元素的内存，需确认signal_set_size函数是否能正确处理大内存分配，是否存在分配失败但返回成功的情况。
• 数据类型赋值损失：输出信号的data是double类型，但卷积中用float类型的t赋值给它，进一步加剧精度损失。

内容的提问来源于stack exchange，提问作者user12312425