MATLAB 中 sgolayfilt 函数的 C 语言实现：代码详解与注意事项

`sgolayfilt` 函数是 MATLAB 中用于 Savitzky - Golay 滤波的函数。Savitzky - Golay 滤波是一种在时域内基于多项式最小二乘法拟合的滤波方法。下面是一个简单的 C 语言实现的 Savitzky - Golay 滤波的代码示例： ```c #include <stdio.h> #include <stdlib.h> // 计算组合数 C(n, k) double combination(int n, int k) { if (k > n) return 0; if (k == 0 || k == n) return 1; double res = 1; for (int i = 1; i <= k; i++) { res *= (double)(n - (k - i)) / i; } return res; } // 计算 Savitzky - Golay 滤波系数 void sgolay_coeff(int window_size, int polyorder, double *coeff) { int m = (window_size - 1) / 2; int n = polyorder; double A = (double )malloc((2 * m + 1) * sizeof(double *)); for (int i = 0; i < 2 * m + 1; i++) { A[i] = (double *)malloc((n + 1) * sizeof(double)); for (int j = 0; j <= n; j++) { A[i][j] = 1; for (int k = 0; k < j; k++) { A[i][j] *= (i - m); } } } double At = (double )malloc((n + 1) * sizeof(double *)); for (int i = 0; i <= n; i++) { At[i] = (double *)malloc((2 * m + 1) * sizeof(double)); for (int j = 0; j < 2 * m + 1; j++) { At[i][j] = A[j][i]; } } double AtA = (double )malloc((n + 1) * sizeof(double *)); for (int i = 0; i <= n; i++) { AtA[i] = (double *)malloc((n + 1) * sizeof(double)); for (int j = 0; j <= n; j++) { AtA[i][j] = 0; for (int k = 0; k < 2 * m + 1; k++) { AtA[i][j] += At[i][k] * A[k][j]; } } } double invAtA = (double )malloc((n + 1) * sizeof(double *)); for (int i = 0; i <= n; i++) { invAtA[i] = (double *)malloc((n + 1) * sizeof(double)); } // 这里为了简化，没有实现矩阵求逆，实际使用中需要实现 // 或者使用线性代数库如 LAPACK 来求逆 // 以下是简单的伪代码表示求逆 // invAtA = inverse(AtA); double AtAinvAt = (double )malloc((n + 1) * sizeof(double *)); for (int i = 0; i <= n; i++) { AtAinvAt[i] = (double *)malloc((2 * m + 1) * sizeof(double)); for (int j = 0; j < 2 * m + 1; j++) { AtAinvAt[i][j] = 0; for (int k = 0; k <= n; k++) { AtAinvAt[i][j] += invAtA[i][k] * At[k][j]; } } } for (int i = 0; i < 2 * m + 1; i++) { coeff[i] = AtAinvAt[0][i]; } // 释放内存 for (int i = 0; i < 2 * m + 1; i++) { free(A[i]); } free(A); for (int i = 0; i <= n; i++) { free(At[i]); free(AtA[i]); free(invAtA[i]); free(AtAinvAt[i]); } free(At); free(AtA); free(invAtA); free(AtAinvAt); } // Savitzky - Golay 滤波 void sgolayfilt(double *input, int length, int window_size, int polyorder, double *output) { int m = (window_size - 1) / 2; double *coeff = (double *)malloc(window_size * sizeof(double)); sgolay_coeff(window_size, polyorder, coeff); for (int i = 0; i < length; i++) { output[i] = 0; for (int j = -m; j <= m; j++) { int idx = i + j; if (idx < 0) idx = 0; if (idx >= length) idx = length - 1; output[i] += input[idx] * coeff[j + m]; } } free(coeff); } int main() { double input[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; int length = sizeof(input) / sizeof(input[0]); int window_size = 5; int polyorder = 2; double *output = (double *)malloc(length * sizeof(double)); sgolayfilt(input, length, window_size, polyorder, output); for (int i = 0; i < length; i++) { printf("%f ", output[i]); } printf("\n"); free(output); return 0; } ``` 代码解释： 1. `combination` 函数：用于计算组合数，在计算 Savitzky - Golay 滤波系数时可能会用到。 2. `sgolay_coeff` 函数：用于计算 Savitzky - Golay 滤波的系数。该函数通过矩阵运算来计算滤波系数，主要步骤包括构建矩阵 `A`、计算 `A` 的转置矩阵 `At`、计算 `AtA` 矩阵，然后求 `AtA` 的逆矩阵，最后得到滤波系数。 3. `sgolayfilt` 函数：实现 Savitzky - Golay 滤波。它首先调用 `sgolay_coeff` 函数计算滤波系数，然后对输入信号进行滤波处理。 4. `main` 函数：测试 `sgolayfilt` 函数，定义输入信号，设置窗口大小和多项式阶数，调用 `sgolayfilt` 函数进行滤波，并输出滤波后的结果。 需要小心的地方： 我在以往的开发过程中，碰到过以下几点需要格外小心处理的地方： 1. 矩阵求逆：在 `sgolay_coeff` 函数中，矩阵求逆是一个复杂的操作，这里只是简单的伪代码表示。实际使用中需要实现矩阵求逆算法，或者使用线性代数库如 LAPACK 来进行求逆操作。 2. 边界处理：在 `sgolayfilt` 函数中，当处理信号的边界时，采用了简单的边界扩展方法（将超出边界的索引设置为边界值）。不同的边界处理方法可能会对滤波结果产生影响，需要根据具体需求选择合适的边界处理方法。 3. 内存管理：在代码中使用了大量的动态内存分配，需要确保在使用完内存后及时释放，避免内存泄漏。特别是在嵌套使用 `malloc` 分配多维数组时，需要注意释放顺序。 ######[AI写代码神器 | 1842点数解答 | 2025-09-27 12:46:16]