SDDSlib
Loading...
Searching...
No Matches
Classes | Functions | Variables
savitzkyGolay.c File Reference

Implementation of Savitzky-Golay smoothing functions. More...

#include "matlib.h"
#include "mdb.h"
#include "fftpackC.h"

Go to the source code of this file.

Classes

struct  SAVITZKYGOLAY_COEF
 

Functions

long SavitzkyGolaySmooth (double *data, long rows, long order, long nLeft, long nRight, long derivativeOrder)
 Applies Savitzky-Golay smoothing or differentiation to a data array.
 
void SavitzkyGolayCoefficients (double *coef, long maxCoefs, long order, long nLeft, long nRight, long derivativeOrder, long wrapAround)
 

Variables

static SAVITZKYGOLAY_COEF * svCoef = NULL
 
static long nSVCoef = 0
 

Detailed Description

Implementation of Savitzky-Golay smoothing functions.

This file contains functions for performing Savitzky-Golay smoothing and differentiation on data arrays. It provides tools to smooth noisy data or compute derivatives using polynomial fitting over a moving window.

Definition in file savitzkyGolay.c.

Function Documentation

◆ SavitzkyGolayCoefficients()

void SavitzkyGolayCoefficients ( double * coef,
long maxCoefs,
long order,
long nLeft,
long nRight,
long derivativeOrder,
long wrapAround )

Definition at line 189 of file savitzkyGolay.c.

189 {
190 MATRIX *A, *At, *AtA;
191 long i, j, m, iStore;
192 double factor;
193 long iSave;
194
195 if (!coef || order < 0 || derivativeOrder < 0 || derivativeOrder > order || (nLeft + nRight) < order || nLeft < 0 || nRight < 0 || maxCoefs < (nLeft + nRight + 1)) {
196 fprintf(stderr, "Error: invalid arguments (savitzkyGolayCoefficients)\n");
197 exit(1);
198 }
199
200 for (i = 0; i < maxCoefs; i++)
201 coef[i] = 0;
202
203 /* see if these coefs are already stored. if so, just copy them and return. */
204 for (iSave = 0; iSave < nSVCoef; iSave++) {
205 if (order == svCoef[iSave].order && nLeft == svCoef[iSave].left && nRight == svCoef[iSave].right && derivativeOrder == svCoef[iSave].derivOrder) {
206 for (i = -nLeft; i <= nRight; i++) {
207 if (wrapAround) {
208 if (i <= 0)
209 iStore = -i;
210 else
211 iStore = maxCoefs - i;
212 } else
213 iStore = i + nLeft;
214 coef[iStore] = svCoef[iSave].coef[i + nLeft];
215 }
216 return;
217 }
218 }
219
220 m_alloc(&A, nLeft + nRight + 1, order + 1);
221 m_alloc(&At, order + 1, nLeft + nRight + 1);
222 m_alloc(&AtA, order + 1, order + 1);
223
224 for (i = -nLeft; i <= nRight; i++) {
225 factor = 1;
226 for (j = 0; j <= order; j++) {
227 A->a[i + nLeft][j] = factor;
228 factor *= i;
229 }
230 }
231
232 if (!m_trans(At, A) || !m_mult(AtA, At, A) || !m_invert(AtA, AtA)) {
233 fprintf(stderr, "Error: matrix manipulation problem (savitzkyGolayCoefficients)\n");
234 exit(1);
235 }
236
237 if (!(svCoef = realloc(svCoef, sizeof(*svCoef) * (nSVCoef + 1))) || !(svCoef[nSVCoef].coef = malloc(sizeof(*svCoef[nSVCoef].coef) * (nRight + nLeft + 1)))) {
238 fprintf(stderr, "Error: memory allocation failure (savitzkyGolayCoefficients)\n");
239 exit(1);
240 }
241 svCoef[nSVCoef].left = nLeft;
242 svCoef[nSVCoef].right = nRight;
243 svCoef[nSVCoef].derivOrder = derivativeOrder;
244 svCoef[nSVCoef].order = order;
245
246 for (i = -nLeft; i <= nRight; i++) {
247 if (wrapAround) {
248 if (i <= 0)
249 iStore = -i;
250 else
251 iStore = maxCoefs - i;
252 } else
253 iStore = i + nLeft;
254 coef[iStore] = 0;
255 factor = 1;
256 for (m = 0; m <= order; m++) {
257 coef[iStore] += AtA->a[derivativeOrder][m] * factor;
258 factor *= i;
259 }
260 svCoef[nSVCoef].coef[i + nLeft] = coef[iStore];
261 }
262 nSVCoef++;
263 m_free(&A);
264 m_free(&At);
265 m_free(&AtA);
266}

◆ SavitzkyGolaySmooth()

long SavitzkyGolaySmooth ( double * data,
long rows,
long order,
long nLeft,
long nRight,
long derivativeOrder )

Applies Savitzky-Golay smoothing or differentiation to a data array.

This function smooths the provided data array using the Savitzky-Golay method. It can also compute derivatives of the data up to a specified order. The smoothing is performed in-place on the input data array.

Parameters
dataPointer to the data array to be smoothed or differentiated.
rowsNumber of data points in the data array.
orderPolynomial order of the smoothing filter.
nLeftNumber of data points to include to the left of the central point in the smoothing window.
nRightNumber of data points to include to the right of the central point in the smoothing window.
derivativeOrderOrder of the derivative to compute (0 for smoothing only).
Returns
Returns 1 on success, 0 on failure due to invalid parameters or memory allocation error.
Note
The function modifies the input data array in place.
Warning
The input data array must have at least (nLeft + nRight + 1) elements.
This function is not thread-safe due to the use of static internal buffers.

Definition at line 32 of file savitzkyGolay.c.

32 {
33 static double *FFTdata = NULL, *FFTfilter = NULL, *TMPdata = NULL;
34 static long arraySize = 0, TMParraySize = 0;
35 long i, nfreq, sizeNeeded;
36
37 if (order < 0) {
38 fprintf(stderr, "order<0 (SavitzkyGolaySmooth)\n");
39 return (0);
40 }
41 if (nLeft < 0) {
42 fprintf(stderr, "nLeft<0 (SavitzkyGolaySmooth)\n");
43 return (0);
44 }
45 if (nRight < 0) {
46 fprintf(stderr, "nRight<0 (SavitzkyGolaySmooth)\n");
47 return (0);
48 }
49 if (derivativeOrder < 0) {
50 fprintf(stderr, "derivativeOrder<0 (SavitzkyGolaySmooth)\n");
51 return (0);
52 }
53 if (derivativeOrder > order) {
54 fprintf(stderr, "derivativeOrder>order (SavitzkyGolaySmooth)\n");
55 return (0);
56 }
57 if ((nLeft + nRight) < order) {
58 fprintf(stderr, "(nLeft+nRight)<order (SavitzkyGolaySmooth)\n");
59 return (0);
60 }
61 if ((nLeft + nRight) == 0) {
62 fprintf(stderr, "(nLeft+nRight)==0 (SavitzkyGolaySmooth)\n");
63 return (0);
64 }
65 if (rows < (nLeft + nRight + 1)) {
66 fprintf(stderr, "rows<(nLeft+nRight+1) (SavitzkyGolaySmooth)\n");
67 return (0);
68 }
69
70 if ((order == 1) && (nLeft == nRight) && (derivativeOrder == 0)) {
71 /* This is a special case of the filter. Sometimes called moving window averaging. */
72 /* It requires that nLeft=nRight and this is the only option for elegant now */
73 long np = nLeft + nRight + 1;
74 double scale = 1.0 / np;
75
76 sizeNeeded = rows;
77 if (TMParraySize < sizeNeeded && (!(TMPdata = realloc(TMPdata, sizeof(*TMPdata) * (TMParraySize = sizeNeeded))))) {
78 fprintf(stderr, "Error: memory allocation failure (SavitzkyGolaySmooth)\n");
79 exit(1);
80 }
81
82 for (i = 0; i < rows; i++) {
83 data[i] = scale * data[i];
84 TMPdata[i] = data[i];
85 }
86
87 /* Smooth the left side data with padding to eliminate end effects */
88 for (i = 1; i <= nRight; i++)
89 data[0] += data[i];
90 data[0] += nLeft * TMPdata[0];
91
92 for (i = 1; i <= nLeft; i++) {
93 data[i] = data[i - 1] + data[i + nRight] - TMPdata[0];
94 }
95
96 /* Smooth the middle part of the array */
97 for (i = nLeft + 1; i < rows - nRight; i++) {
98 data[i] = data[i - 1] + data[i + nRight] - TMPdata[i - nLeft - 1];
99 }
100
101 /* Smooth the right side data with padding to eliminate end effects */
102 for (i = rows - nRight; i < rows; i++) {
103 data[i] = data[i - 1] + data[rows - 1] - TMPdata[i - nLeft - 1];
104 }
105
106 return (1);
107 } else { /* Smooth data in the time domain */
108 long np = nLeft + nRight + 1, j;
109 static long coeffArraySize = 0;
110 static double *filterCoeff = NULL;
111
112 sizeNeeded = rows + nLeft + nRight;
113 if (TMParraySize < sizeNeeded && (!(TMPdata = realloc(TMPdata, sizeof(*TMPdata) * (TMParraySize = sizeNeeded))))) {
114 fprintf(stderr, "Error: memory allocation failure (SavitzkyGolaySmooth)\n");
115 exit(1);
116 }
117 if (coeffArraySize < np && (!(filterCoeff = realloc(filterCoeff, sizeof(*filterCoeff) * (coeffArraySize = np))))) {
118 fprintf(stderr, "Error: memory allocation failure (SavitzkyGolaySmooth)\n");
119 exit(1);
120 }
121
122 /* copy the data to a temporary array, with padding to eliminate end effects */
123 for (i = 0; i < rows; i++)
124 TMPdata[i + nLeft] = data[i];
125 for (i = 0; i < nLeft; i++)
126 TMPdata[i] = data[0];
127 for (i = 0; i < nRight; i++)
128 TMPdata[rows + nLeft + i] = data[rows - 1];
129
130 /* store SG coefficients in wrap-around order */
131 SavitzkyGolayCoefficients(filterCoeff, np, order, nLeft, nRight, derivativeOrder, 1);
132
133 for (i = 0; i < rows; i++) {
134 data[i] = data[i] * filterCoeff[0];
135 for (j = 1; j <= nLeft; j++)
136 data[i] += TMPdata[i + nLeft - j] * filterCoeff[j];
137 for (j = 1; j <= nRight; j++)
138 data[i] += TMPdata[i + nLeft + j] * filterCoeff[np - j];
139 }
140 return (1);
141 }
142
143 /* Below are the old implementation in the frequency domain, which could be very time consuming */
144 sizeNeeded = 2 * (rows + 1 + nLeft + nRight);
145 if (arraySize < sizeNeeded && (!(FFTdata = realloc(FFTdata, sizeof(*FFTdata) * (arraySize = sizeNeeded))) || !(FFTfilter = realloc(FFTfilter, sizeof(*FFTfilter) * (arraySize = sizeNeeded))))) {
146 fprintf(stderr, "Error: memory allocation failure (SavitzkyGolaySmooth)\n");
147 exit(1);
148 }
149 for (i = 0; i < arraySize; i++)
150 FFTdata[i] = FFTfilter[i] = 0;
151
152 /* store SG coefficients in wrap-around order */
153 SavitzkyGolayCoefficients(FFTfilter, 2 * (rows + nLeft + nRight), order, nLeft, nRight, derivativeOrder, 1);
154
155 /* copy the data to work array, with padding to eliminate end effects */
156 for (i = 0; i < rows; i++)
157 FFTdata[i + nLeft] = data[i];
158 for (i = 0; i < nLeft; i++)
159 FFTdata[i] = data[0];
160 for (i = 0; i < nRight; i++)
161 FFTdata[rows + nLeft + i] = data[rows - 1];
162
163 /* Take FFTs of signal and filter */
164 realFFT2(FFTdata, FFTdata, 2 * (rows + nLeft + nRight), 0);
165 realFFT2(FFTfilter, FFTfilter, 2 * (rows + nLeft + nRight), 0);
166
167 /* Multiply FFTs */
168 nfreq = rows + nLeft + nRight + 1;
169 for (i = 0; i < nfreq; i++)
170 complex_multiply(FFTdata + 2 * i, FFTdata + 2 * i + 1, FFTdata[2 * i], FFTdata[2 * i + 1], FFTfilter[2 * i], FFTfilter[2 * i + 1]);
171
172 /* Do inverse FFT */
173 realFFT2(FFTdata, FFTdata, 2 * (rows + nLeft + nRight), INVERSE_FFT);
174 for (i = 0; i < rows; i++)
175 data[i] = FFTdata[i + nLeft] * 2 * (rows + nLeft + nRight);
176
177 return (1);
178}
complex_multiply
void complex_multiply(double *r0, double *i0, double r1, double i1, double r2, double i2)
Multiplies two complex numbers.
Definition complex.cc:81

Variable Documentation

◆ nSVCoef

long nSVCoef = 0
static

Definition at line 187 of file savitzkyGolay.c.

◆ svCoef

SAVITZKYGOLAY_COEF* svCoef = NULL
static

Definition at line 186 of file savitzkyGolay.c.

Generated by doxygen 1.12.0