请用C语言写一份完整的Savitzky-Golay filter算法

# Savitzky-Golay Filter in C Below is a complete implementation of the Savitzky-Golay filter in C. This implementation includes functions for calculating the filter coefficients and applying the filter to a data set. ```c #include <stdio.h> #include <stdlib.h> #include <math.h> // Function to calculate the Savitzky-Golay coefficients void savitzky_golay_coeff(int window_size, int poly_order, double* coeff) { int half_window = window_size / 2; double A[(2*poly_order+1) * window_size]; double AtA[(2*poly_order+1) * (2*poly_order+1)]; double invAtA[(2*poly_order+1) * (2*poly_order+1)]; // Construct the design matrix A for (int i = -half_window; i <= half_window; i++) { for (int j = 0; j <= poly_order; j++) { A[(i + half_window) * (poly_order + 1) + j] = pow(i, j); } } // Compute A^T * A for (int i = 0; i <= poly_order; i++) { for (int j = 0; j <= poly_order; j++) { AtA[i * (poly_order + 1) + j] = 0.0; for (int k = 0; k < window_size; k++) { AtA[i * (poly_order + 1) + j] += A[k * (poly_order + 1) + i] * A[k * (poly_order + 1) + j]; } } } // Compute inverse of A^T * A (using simple inversion for small matrices) double det = 0.0; if (poly_order == 2) { // For quadratic polynomial det = AtA[0]*(AtA[4]*AtA[8]-AtA[5]*AtA[7]) - AtA[1]*(AtA[3]*AtA[8]-AtA[5]*AtA[6]) + AtA[2]*(AtA[3]*AtA[7]-AtA[4]*AtA[6]); invAtA[0] = (AtA[4]*AtA[8] - AtA[5]*AtA[7]) / det; invAtA[1] = (AtA[2]*AtA[7] - AtA[1]*AtA[8]) / det; invAtA[2] = (AtA[1]*AtA[5] - AtA[2]*AtA[4]) / det; invAtA[3] = (AtA[5]*AtA[6] - AtA[3]*AtA[8]) / det; invAtA[4] = (AtA[0]*AtA[8] - AtA[2]*AtA[6]) / det; invAtA[5] = (AtA[2]*AtA[3] - AtA[0]*AtA[5]) / det; invAtA[6] = (AtA[3]*AtA[7] - AtA[4]*AtA[6]) / det; invAtA[7] = (AtA[1]*AtA[6] - AtA[0]*AtA[7]) / det; invAtA[8] = (AtA[0]*AtA[4] - AtA[1]*AtA[3]) / det; } else { // For linear polynomial (simpler case) det = AtA[0]*AtA[3] - AtA[1]*AtA[1]; invAtA[0] = AtA[3] / det; invAtA[1] = -AtA[1] / det; invAtA[2] = -AtA[1] / det; invAtA[3] = AtA[0] / det; } // Compute the coefficients for the middle point (smoothing) for (int i = 0; i < window_size; i++) { coeff[i] = 0.0; for (int j = 0; j <= poly_order; j++) { for (int k = 0; k <= poly_order; k++) { coeff[i] += A[i * (poly_order + 1) + k] * invAtA[k * (poly_order + 1) + j] * A[half_window * (poly_order + 1) + j]; } } } } // Function to apply Savitzky-Golay filter to data void apply_savitzky_golay(double* input, double* output, int data_size, int window_size, int poly_order) { double* coeff = (double*)malloc(window_size * sizeof(double)); savitzky_golay_coeff(window_size, poly_order, coeff); int half_window = window_size / 2; // Handle boundaries by simple copying (could be improved) for (int i = 0; i < half_window; i++) { output[i] = input[i]; output[data_size - 1 - i] = input[data_size - 1 - i]; } // Apply filter to middle points for (int i = half_window; i < data_size - half_window; i++) { output[i] = 0.0; for (int j = -half_window; j <= half_window; j++) { output[i] += coeff[j + half_window] * input[i + j]; } } free(coeff); } int main() { // Example usage const int data_size = 20; double input[data_size] = {1.2, 1.5, 1.7, 1.4, 1.6, 1.8, 2.0, 2.1, 2.3, 2.2, 2.0, 1.9, 1.7, 1.6, 1.5, 1.4, 1.3, 1.2, 1.1, 1.0}; double output[data_size]; int window_size = 5; // Must be odd int poly_order = 2; // Typically 2 (quadratic) or 3 (cubic) apply_savitzky_golay(input, output, data_size, window_size, poly_order); printf("Original Data\tSmoothed Data\n"); for (int i = 0; i < data_size; i++) { printf("%.2f\t\t%.2f\n", input[i], output[i]); } return 0; } ``` ## 代码说明 1. `savitzky_golay_coeff`函数计算Savitzky-Golay滤波器的系数矩阵 2. `apply_savitzky_golay`函数将滤波器应用于输入数据 3. 主函数演示了如何使用这些函数处理示例数据 ## 使用注意事项 1. 窗口大小(window_size)必须是奇数 2. 多项式阶数(poly_order)通常选择2(二次)或3(三次) 3. 边界处理采用了简单的复制方法，可以根据需要改进 4. 矩阵求逆部分目前只实现了线性和二次多项式的情况 ## 额外需要注意的逻辑点 1. 窗口大小和多项式阶数的选择会影响滤波效果，窗口越大平滑效果越强但可能丢失细节 2. 边界处理可以改进为使用较小的窗口或镜像数据等方法 3. 对于更高阶多项式，需要实现更通用的矩阵求逆算法 [2025-07-22 17:58:35 | AI写代码神器 | 1730点数解答]