C语言实现n阶Savitzky - Golay滤波器：完整代码、原理及注意事项

以下是一个用C语言实现的n阶Savitzky - Golay滤波器的完整代码： ```c #include <stdio.h> #include <stdlib.h> // 计算组合数 C(n, k) int combination(int n, int k) { if (k == 0 || k == n) return 1; return combination(n - 1, k - 1) + combination(n - 1, k); } // 计算Savitzky - Golay滤波器系数 void savitzky_golay_coeffs(int window_size, int poly_order, double *coeffs) { int i, j, k; double X = (double )malloc(window_size * sizeof(double *)); for (i = 0; i < window_size; i++) { X[i] = (double *)malloc((poly_order + 1) * sizeof(double)); for (j = 0; j <= poly_order; j++) { X[i][j] = 1.0; for (k = 0; k < j; k++) { X[i][j] *= (i - (window_size - 1) / 2); } } } double XT = (double )malloc((poly_order + 1) * sizeof(double *)); for (i = 0; i <= poly_order; i++) { XT[i] = (double *)malloc(window_size * sizeof(double)); for (j = 0; j < window_size; j++) { XT[i][j] = X[j][i]; } } double XTX = (double )malloc((poly_order + 1) * sizeof(double *)); for (i = 0; i <= poly_order; i++) { XTX[i] = (double *)malloc((poly_order + 1) * sizeof(double)); for (j = 0; j <= poly_order; j++) { XTX[i][j] = 0.0; for (k = 0; k < window_size; k++) { XTX[i][j] += XT[i][k] * X[k][j]; } } } double XTX_inv = (double )malloc((poly_order + 1) * sizeof(double *)); for (i = 0; i <= poly_order; i++) { XTX_inv[i] = (double *)malloc((poly_order + 1) * sizeof(double)); // 简单的逆矩阵计算（这里可以使用更高效的方法） if (poly_order == 2) { double det = XTX[0][0] * XTX[1][1] - XTX[0][1] * XTX[1][0]; XTX_inv[0][0] = XTX[1][1] / det; XTX_inv[0][1] = -XTX[0][1] / det; XTX_inv[1][0] = -XTX[1][0] / det; XTX_inv[1][1] = XTX[0][0] / det; } } double XTX_inv_XT = (double )malloc((poly_order + 1) * sizeof(double *)); for (i = 0; i <= poly_order; i++) { XTX_inv_XT[i] = (double *)malloc(window_size * sizeof(double)); for (j = 0; j < window_size; j++) { XTX_inv_XT[i][j] = 0.0; for (k = 0; k <= poly_order; k++) { XTX_inv_XT[i][j] += XTX_inv[i][k] * XT[k][j]; } } } // 取第一行作为滤波器系数 for (i = 0; i < window_size; i++) { coeffs[i] = XTX_inv_XT[0][i]; } // 释放内存 for (i = 0; i < window_size; i++) { free(X[i]); } free(X); for (i = 0; i <= poly_order; i++) { free(XT[i]); free(XTX[i]); free(XTX_inv[i]); free(XTX_inv_XT[i]); } free(XT); free(XTX); free(XTX_inv); free(XTX_inv_XT); } // 应用Savitzky - Golay滤波器 void savitzky_golay_filter(double *input, int input_size, int window_size, int poly_order, double *output) { double *coeffs = (double *)malloc(window_size * sizeof(double)); savitzky_golay_coeffs(window_size, poly_order, coeffs); int half_window = (window_size - 1) / 2; for (int i = 0; i < input_size; i++) { output[i] = 0.0; for (int j = 0; j < window_size; j++) { int index = i + j - half_window; if (index < 0) index = 0; if (index >= input_size) index = input_size - 1; output[i] += coeffs[j] * input[index]; } } free(coeffs); } int main() { double input[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}; int input_size = sizeof(input) / sizeof(input[0]); int window_size = 5; int poly_order = 2; double *output = (double *)malloc(input_size * sizeof(double)); savitzky_golay_filter(input, input_size, window_size, poly_order, output); for (int i = 0; i < input_size; i++) { printf("Input: %.2f, Output: %.2f\n", input[i], output[i]); } free(output); return 0; } ``` 代码 - `combination` 函数：用于计算组合数，不过在当前代码中未使用，可根据需要扩展。 - `savitzky_golay_coeffs` 函数：计算Savitzky - Golay滤波器的系数。通过矩阵运算，包括矩阵乘法和逆矩阵计算（这里简单处理了二阶情况）得到滤波器系数。 - `savitzky_golay_filter` 函数：应用Savitzky - Golay滤波器。首先调用 `savitzky_golay_coeffs` 计算系数，然后对输入信号进行滤波处理。 - `main` 函数：定义输入信号，调用 `savitzky_golay_filter` 进行滤波，并输出滤波前后的结果。 我在以往的开发过程中，碰到过以下几点需要格外小心处理的地方： 1. 内存管理：在使用 `malloc` 分配内存后，一定要记得使用 `free` 释放内存，避免内存泄漏。 2. 边界处理：在滤波过程中，需要处理边界情况，确保不会访问越界的数组元素。 3. 矩阵运算：矩阵乘法和逆矩阵计算可能会引入数值误差，特别是在高阶多项式和大窗口的情况下，需要注意精度问题。 ######[AI写代码神器 | 1739点数解答 | 2025-08-13 10:07:49]