File indexing completed on 2024-05-12 03:47:54
0001 /* 0002 File : nsl_smooth.c 0003 Project : LabPlot 0004 Description : NSL smooth functions 0005 -------------------------------------------------------------------- 0006 SPDX-FileCopyrightText: 2010 Knut Franke <knut.franke@gmx.de> 0007 SPDX-FileCopyrightText: 2016 Stefan Gerlach <stefan.gerlach@uni.kn> 0008 0009 SPDX-License-Identifier: GPL-2.0-or-later 0010 */ 0011 0012 #include "nsl_smooth.h" 0013 #include "nsl_common.h" 0014 #include "nsl_sf_kernel.h" 0015 #include "nsl_stats.h" 0016 #include <gsl/gsl_blas.h> 0017 #include <gsl/gsl_linalg.h> 0018 #include <gsl/gsl_math.h> 0019 #include <gsl/gsl_sf_gamma.h> /* gsl_sf_choose */ 0020 0021 const char* nsl_smooth_type_name[] = {i18n("Moving Average (Central)"), i18n("Moving Average (Lagged)"), i18n("Percentile"), i18n("Savitzky-Golay")}; 0022 const char* nsl_smooth_pad_mode_name[] = {i18n("None"), i18n("Interpolating"), i18n("Mirror"), i18n("Nearest"), i18n("Constant"), i18n("Periodic")}; 0023 const char* nsl_smooth_weight_type_name[] = {i18n("Uniform (Rectangular)"), 0024 i18n("Triangular"), 0025 i18n("Binomial"), 0026 i18n("Parabolic (Epanechnikov)"), 0027 i18n("Quartic (Biweight)"), 0028 i18n("Triweight"), 0029 i18n("Tricube"), 0030 i18n("Cosine")}; 0031 double nsl_smooth_pad_constant_lvalue = 0.0, nsl_smooth_pad_constant_rvalue = 0.0; 0032 0033 int nsl_smooth_moving_average(double* data, size_t n, size_t points, nsl_smooth_weight_type weight, nsl_smooth_pad_mode mode) { 0034 if (n == 0 || points == 0) 0035 return -1; 0036 0037 size_t i, j; 0038 double* result = (double*)malloc(n * sizeof(double)); 0039 for (i = 0; i < n; i++) 0040 result[i] = 0; 0041 0042 for (i = 0; i < n; i++) { 0043 size_t np = points; 0044 size_t half = (points - 1) / 2; 0045 if (mode == nsl_smooth_pad_none) { /* reduce points */ 0046 half = GSL_MIN(GSL_MIN((points - 1) / 2, i), n - i - 1); 0047 np = 2 * half + 1; 0048 } 0049 0050 /* weight */ 0051 double sum = 0.0; 0052 double* w = (double*)malloc(np * sizeof(double)); 0053 switch (weight) { 0054 case nsl_smooth_weight_uniform: 0055 for (j = 0; j < np; j++) 0056 w[j] = 1. / np; 0057 break; 0058 case nsl_smooth_weight_triangular: 0059 sum = gsl_pow_2((double)(np + 1) / 2); 0060 for (j = 0; j < np; j++) 0061 w[j] = GSL_MIN(j + 1, np - j) / sum; 0062 break; 0063 case nsl_smooth_weight_binomial: 0064 sum = (np - 1) / 2.; 0065 for (j = 0; j < np; j++) 0066 w[j] = gsl_sf_choose((unsigned int)(2 * sum), (unsigned int)((sum + fabs(j - sum)) / pow(4., sum))); 0067 break; 0068 case nsl_smooth_weight_parabolic: 0069 for (j = 0; j < np; j++) { 0070 w[j] = nsl_sf_kernel_parabolic(2. * (j - (np - 1) / 2.) / (np + 1)); 0071 sum += w[j]; 0072 } 0073 for (j = 0; j < np; j++) 0074 w[j] /= sum; 0075 break; 0076 case nsl_smooth_weight_quartic: 0077 for (j = 0; j < np; j++) { 0078 w[j] = nsl_sf_kernel_quartic(2. * (j - (np - 1) / 2.) / (np + 1)); 0079 sum += w[j]; 0080 } 0081 for (j = 0; j < np; j++) 0082 w[j] /= sum; 0083 break; 0084 case nsl_smooth_weight_triweight: 0085 for (j = 0; j < np; j++) { 0086 w[j] = nsl_sf_kernel_triweight(2. * (j - (np - 1) / 2.) / (np + 1)); 0087 sum += w[j]; 0088 } 0089 for (j = 0; j < np; j++) 0090 w[j] /= sum; 0091 break; 0092 case nsl_smooth_weight_tricube: 0093 for (j = 0; j < np; j++) { 0094 w[j] = nsl_sf_kernel_tricube(2. * (j - (np - 1) / 2.) / (np + 1)); 0095 sum += w[j]; 0096 } 0097 for (j = 0; j < np; j++) 0098 w[j] /= sum; 0099 break; 0100 case nsl_smooth_weight_cosine: 0101 for (j = 0; j < np; j++) { 0102 w[j] = nsl_sf_kernel_cosine((j - (np - 1) / 2.) / ((np + 1) / 2.)); 0103 sum += w[j]; 0104 } 0105 for (j = 0; j < np; j++) 0106 w[j] /= sum; 0107 break; 0108 } 0109 0110 /*printf("(%d) w:",i); 0111 for(j=0;j<np;j++) 0112 printf(" %g",w[j]); 0113 printf(" (half=%d) index = ",half);*/ 0114 0115 /* calculate weighted average */ 0116 for (j = 0; j < np; j++) { 0117 int index = (int)(i - half + j); 0118 switch (mode) { 0119 case nsl_smooth_pad_none: 0120 result[i] += w[j] * data[index]; 0121 break; 0122 case nsl_smooth_pad_interp: 0123 printf("not implemented yet\n"); 0124 break; 0125 case nsl_smooth_pad_mirror: 0126 index = abs((int)(i - half + j)); 0127 /*printf(" %d",GSL_MIN(index,2*(n-1)-index));*/ 0128 result[i] += w[j] * data[GSL_MIN(index, 2 * ((int)n - 1) - index)]; 0129 break; 0130 case nsl_smooth_pad_nearest: 0131 /*printf(" %d",GSL_MIN(n-1,GSL_MAX(0,index)));*/ 0132 result[i] += w[j] * data[GSL_MIN((int)n - 1, GSL_MAX(0, index))]; 0133 break; 0134 case nsl_smooth_pad_constant: 0135 if (index < 0) 0136 result[i] += w[j] * nsl_smooth_pad_constant_lvalue; 0137 else if (index > (int)n - 1) 0138 result[i] += w[j] * nsl_smooth_pad_constant_rvalue; 0139 else 0140 result[i] += w[j] * data[index]; 0141 break; 0142 case nsl_smooth_pad_periodic: 0143 if (index < 0) 0144 index = index + (int)n; 0145 else if (index > (int)n - 1) 0146 index = index - (int)n; 0147 result[i] += w[j] * data[index]; 0148 break; 0149 } 0150 } 0151 /*puts("");*/ 0152 free(w); 0153 } 0154 0155 for (i = 0; i < n; i++) 0156 data[i] = result[i]; 0157 free(result); 0158 0159 return 0; 0160 } 0161 0162 int nsl_smooth_moving_average_lagged(double* data, size_t n, size_t points, nsl_smooth_weight_type weight, nsl_smooth_pad_mode mode) { 0163 if (n == 0 || points == 0) 0164 return -1; 0165 0166 size_t i, j; 0167 double* result = (double*)malloc(n * sizeof(double)); 0168 for (i = 0; i < n; i++) 0169 result[i] = 0; 0170 0171 for (i = 0; i < n; i++) { 0172 size_t np = points; 0173 size_t half = (points - 1) / 2; 0174 if (mode == nsl_smooth_pad_none) { /* reduce points */ 0175 np = GSL_MIN(points, i + 1); 0176 half = np - 1; 0177 } 0178 0179 /* weight */ 0180 double sum = 0.0, *w = (double*)malloc(np * sizeof(double)); 0181 switch (weight) { 0182 case nsl_smooth_weight_uniform: 0183 for (j = 0; j < np; j++) 0184 w[j] = 1. / np; 0185 break; 0186 case nsl_smooth_weight_triangular: 0187 sum = np * (double)(np + 1) / 2; 0188 for (j = 0; j < np; j++) 0189 w[j] = (j + 1) / sum; 0190 break; 0191 case nsl_smooth_weight_binomial: 0192 for (j = 0; j < np; j++) { 0193 w[j] = gsl_sf_choose((unsigned int)(2 * (np - 1)), (unsigned int)j); 0194 sum += w[j]; 0195 } 0196 for (j = 0; j < np; j++) 0197 w[j] /= sum; 0198 break; 0199 case nsl_smooth_weight_parabolic: 0200 for (j = 0; j < np; j++) { 0201 w[j] = nsl_sf_kernel_parabolic(1. - (1 + j) / (double)np); 0202 sum += w[j]; 0203 } 0204 for (j = 0; j < np; j++) 0205 w[j] /= sum; 0206 break; 0207 case nsl_smooth_weight_quartic: 0208 for (j = 0; j < np; j++) { 0209 w[j] = nsl_sf_kernel_quartic(1. - (1 + j) / (double)np); 0210 sum += w[j]; 0211 } 0212 for (j = 0; j < np; j++) 0213 w[j] /= sum; 0214 break; 0215 case nsl_smooth_weight_triweight: 0216 for (j = 0; j < np; j++) { 0217 w[j] = nsl_sf_kernel_triweight(1. - (1 + j) / (double)np); 0218 sum += w[j]; 0219 } 0220 for (j = 0; j < np; j++) 0221 w[j] /= sum; 0222 break; 0223 case nsl_smooth_weight_tricube: 0224 for (j = 0; j < np; j++) { 0225 w[j] = nsl_sf_kernel_tricube(1. - (1 + j) / (double)np); 0226 sum += w[j]; 0227 } 0228 for (j = 0; j < np; j++) 0229 w[j] /= sum; 0230 break; 0231 case nsl_smooth_weight_cosine: 0232 for (j = 0; j < np; j++) { 0233 w[j] = nsl_sf_kernel_cosine((np - 1 - j) / (double)np); 0234 sum += w[j]; 0235 } 0236 for (j = 0; j < np; j++) 0237 w[j] /= sum; 0238 break; 0239 } 0240 0241 /*printf("(%d) w:",i); 0242 for(j=0;j<np;j++) 0243 printf(" %g",w[j]); 0244 printf(" (half=%d) index = ",half);*/ 0245 0246 /* calculate weighted average */ 0247 for (j = 0; j < np; j++) { 0248 int index = (int)(i - np + 1 + j); 0249 switch (mode) { 0250 case nsl_smooth_pad_none: 0251 result[i] += w[j] * data[i - half + j]; 0252 /*printf(" %d",index);*/ 0253 break; 0254 case nsl_smooth_pad_interp: 0255 printf("not implemented yet\n"); 0256 break; 0257 case nsl_smooth_pad_mirror: 0258 index = abs(index); 0259 /*printf(" %d",index);*/ 0260 result[i] += w[j] * data[index]; 0261 break; 0262 case nsl_smooth_pad_nearest: 0263 /*printf(" %d", index);*/ 0264 result[i] += w[j] * data[GSL_MAX(0, index)]; 0265 break; 0266 case nsl_smooth_pad_constant: 0267 if (index < 0) 0268 result[i] += w[j] * nsl_smooth_pad_constant_lvalue; 0269 else 0270 result[i] += w[j] * data[index]; 0271 0272 break; 0273 case nsl_smooth_pad_periodic: 0274 if (index < 0) 0275 index += (int)n; 0276 /*printf(" %d",index);*/ 0277 result[i] += w[j] * data[index]; 0278 break; 0279 } 0280 } 0281 /*puts("");*/ 0282 free(w); 0283 } 0284 0285 for (i = 0; i < n; i++) 0286 data[i] = result[i]; 0287 free(result); 0288 0289 return 0; 0290 } 0291 0292 int nsl_smooth_percentile(double* data, size_t n, size_t points, double percentile, nsl_smooth_pad_mode mode) { 0293 if (n == 0 || points == 0) 0294 return -1; 0295 0296 size_t i, j; 0297 double* result = (double*)malloc(n * sizeof(double)); 0298 0299 for (i = 0; i < n; i++) { 0300 size_t np = points; 0301 size_t half = (points - 1) / 2; 0302 if (mode == nsl_smooth_pad_none) { /* reduce points */ 0303 half = GSL_MIN(GSL_MIN((points - 1) / 2, i), n - i - 1); 0304 np = 2 * half + 1; 0305 } 0306 0307 double* values = (double*)malloc(np * sizeof(double)); 0308 for (j = 0; j < np; j++) { 0309 int index = (int)(i - half + j); 0310 switch (mode) { 0311 case nsl_smooth_pad_none: 0312 /*printf(" %d",index);*/ 0313 values[j] = data[index]; 0314 break; 0315 case nsl_smooth_pad_interp: 0316 printf("not implemented yet\n"); 0317 break; 0318 case nsl_smooth_pad_mirror: 0319 index = abs(index); 0320 /*printf(" %d",GSL_MIN(index,2*(n-1)-index));*/ 0321 values[j] = data[GSL_MIN(index, 2 * ((int)n - 1) - index)]; 0322 break; 0323 case nsl_smooth_pad_nearest: 0324 /*printf(" %d",GSL_MIN(n-1,GSL_MAX(0,index)));*/ 0325 values[j] = data[GSL_MIN((int)n - 1, GSL_MAX(0, index))]; 0326 break; 0327 case nsl_smooth_pad_constant: 0328 if (index < 0) 0329 values[j] = nsl_smooth_pad_constant_lvalue; 0330 else if (index > (int)n - 1) 0331 values[j] = nsl_smooth_pad_constant_rvalue; 0332 else 0333 values[j] = data[index]; 0334 break; 0335 case nsl_smooth_pad_periodic: 0336 if (index < 0) 0337 index = index + (int)n; 0338 else if (index > (int)n - 1) 0339 index = index - (int)n; 0340 /*printf(" %d",index);*/ 0341 values[j] = data[index]; 0342 break; 0343 } 0344 } 0345 /*puts("");*/ 0346 0347 /*using type 7 as default */ 0348 result[i] = nsl_stats_quantile(values, 1, np, percentile, nsl_stats_quantile_type7); 0349 free(values); 0350 } 0351 0352 for (i = 0; i < n; i++) 0353 data[i] = result[i]; 0354 free(result); 0355 0356 return 0; 0357 } 0358 0359 /* taken from SciDAVis */ 0360 int nsl_smooth_savgol_coeff(size_t points, int order, gsl_matrix* h) { 0361 size_t i; 0362 int j, error = 0; 0363 0364 /* compute Vandermonde matrix */ 0365 gsl_matrix* vandermonde = gsl_matrix_alloc(points, order + 1); 0366 for (i = 0; i < points; ++i) { 0367 gsl_matrix_set(vandermonde, i, 0, 1.0); 0368 for (j = 1; j <= order; ++j) 0369 gsl_matrix_set(vandermonde, i, j, gsl_matrix_get(vandermonde, i, j - 1) * i); 0370 } 0371 0372 /* compute V^TV */ 0373 gsl_matrix* vtv = gsl_matrix_alloc(order + 1, order + 1); 0374 error = gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, vandermonde, vandermonde, 0.0, vtv); 0375 0376 if (!error) { 0377 /* compute (V^TV)^(-1) using LU decomposition */ 0378 gsl_permutation* p = gsl_permutation_alloc(order + 1); 0379 int signum; 0380 error = gsl_linalg_LU_decomp(vtv, p, &signum); 0381 0382 if (!error) { 0383 gsl_matrix* vtv_inv = gsl_matrix_alloc(order + 1, order + 1); 0384 error = gsl_linalg_LU_invert(vtv, p, vtv_inv); 0385 if (!error) { 0386 /* compute (V^TV)^(-1)V^T */ 0387 gsl_matrix* vtv_inv_vt = gsl_matrix_alloc(order + 1, points); 0388 error = gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, vtv_inv, vandermonde, 0.0, vtv_inv_vt); 0389 0390 if (!error) { 0391 /* finally, compute H = V(V^TV)^(-1)V^T */ 0392 error = gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, vandermonde, vtv_inv_vt, 0.0, h); 0393 } 0394 gsl_matrix_free(vtv_inv_vt); 0395 } 0396 gsl_matrix_free(vtv_inv); 0397 } 0398 gsl_permutation_free(p); 0399 } 0400 gsl_matrix_free(vtv); 0401 gsl_matrix_free(vandermonde); 0402 0403 return error; 0404 } 0405 0406 void nsl_smooth_pad_constant_set(double lvalue, double rvalue) { 0407 nsl_smooth_pad_constant_lvalue = lvalue; 0408 nsl_smooth_pad_constant_rvalue = rvalue; 0409 } 0410 0411 int nsl_smooth_savgol(double* data, size_t n, size_t points, int order, nsl_smooth_pad_mode mode) { 0412 size_t i, k; 0413 int error = 0; 0414 size_t half = (points - 1) / 2; /* n//2 */ 0415 0416 if (points > n) { 0417 printf("Tried to smooth over more points (points=%d) than given as input (%d).", (int)points, (int)n); 0418 return -1; 0419 } 0420 if (order < 1 || (size_t)order > points - 1) { 0421 printf("The polynomial order must be between 1 and %d (%d given).", (int)(points - 1), order); 0422 return -2; 0423 } 0424 0425 /* Savitzky-Golay coefficient matrix, y' = H y */ 0426 gsl_matrix* h = gsl_matrix_alloc(points, points); 0427 error = nsl_smooth_savgol_coeff(points, order, h); 0428 if (error) { 0429 printf("Internal error in Savitzky-Golay algorithm:\n%s", gsl_strerror(error)); 0430 gsl_matrix_free(h); 0431 return error; 0432 } 0433 0434 double* result = (double*)malloc(n * sizeof(double)); 0435 for (i = 0; i < n; i++) 0436 result[i] = 0; 0437 0438 /* left edge */ 0439 if (mode == nsl_smooth_pad_none) { 0440 for (i = 0; i < half; i++) { 0441 /*reduce points and order*/ 0442 size_t rpoints = 2 * i + 1; 0443 int rorder = GSL_MIN(order, (int)(rpoints - GSL_MIN(rpoints, 2))); 0444 0445 gsl_matrix* rh = gsl_matrix_alloc(rpoints, rpoints); 0446 error = nsl_smooth_savgol_coeff(rpoints, rorder, rh); 0447 if (error) { 0448 printf("Internal error in Savitzky-Golay algorithm:\n%s", gsl_strerror(error)); 0449 gsl_matrix_free(rh); 0450 free(result); 0451 return error; 0452 } 0453 0454 for (k = 0; k < rpoints; k++) 0455 result[i] += gsl_matrix_get(rh, i, k) * data[k]; 0456 } 0457 } else { 0458 for (i = 0; i < half; i++) { 0459 for (k = 0; k < points; k++) 0460 switch (mode) { 0461 case nsl_smooth_pad_interp: 0462 result[i] += gsl_matrix_get(h, i, k) * data[k]; 0463 break; 0464 case nsl_smooth_pad_mirror: 0465 result[i] += gsl_matrix_get(h, half, k) * data[abs((int)(k + i - half))]; 0466 break; 0467 case nsl_smooth_pad_nearest: 0468 result[i] += gsl_matrix_get(h, half, k) * data[i + k - GSL_MIN(half, i + k)]; 0469 break; 0470 case nsl_smooth_pad_constant: 0471 if (k < half - i) 0472 result[i] += gsl_matrix_get(h, half, k) * nsl_smooth_pad_constant_lvalue; 0473 else 0474 result[i] += gsl_matrix_get(h, half, k) * data[i - half + k]; 0475 break; 0476 case nsl_smooth_pad_periodic: 0477 result[i] += gsl_matrix_get(h, half, k) * data[k < half - i ? n + i + k - half : i - half + k]; 0478 case nsl_smooth_pad_none: 0479 break; 0480 } 0481 } 0482 } 0483 0484 /* central part: convolve with fixed row of h */ 0485 for (i = half; i < n - half; i++) 0486 for (k = 0; k < points; k++) 0487 result[i] += gsl_matrix_get(h, half, k) * data[i - half + k]; 0488 0489 /* right edge */ 0490 if (mode == nsl_smooth_pad_none) { 0491 for (i = n - half; i < n; i++) { 0492 /*reduce points and order*/ 0493 size_t rpoints = 2 * (n - 1 - i) + 1; 0494 int rorder = (int)GSL_MIN((size_t)order, rpoints - GSL_MIN(2, rpoints)); 0495 0496 gsl_matrix* rh = gsl_matrix_alloc(rpoints, rpoints); 0497 error = nsl_smooth_savgol_coeff(rpoints, rorder, rh); 0498 if (error) { 0499 printf("Internal error in Savitzky-Golay algorithm:\n%s", gsl_strerror(error)); 0500 gsl_matrix_free(rh); 0501 free(result); 0502 return error; 0503 } 0504 0505 for (k = 0; k < rpoints; k++) 0506 result[i] += gsl_matrix_get(rh, n - 1 - i, k) * data[n - rpoints + k]; 0507 } 0508 } else { 0509 for (i = n - half; i < n; i++) { 0510 for (k = 0; k < points; k++) 0511 switch (mode) { 0512 case nsl_smooth_pad_interp: 0513 result[i] += gsl_matrix_get(h, points - n + i, k) * data[n - points + k]; 0514 break; 0515 case nsl_smooth_pad_mirror: 0516 result[i] += gsl_matrix_get(h, half, k) * data[n - 1 - abs((int)(k + 1 + i - n - half))]; 0517 break; 0518 case nsl_smooth_pad_nearest: 0519 result[i] += gsl_matrix_get(h, half, k) * data[GSL_MIN(i - half + k, n - 1)]; 0520 break; 0521 case nsl_smooth_pad_constant: 0522 if (k < n - i + half) 0523 result[i] += gsl_matrix_get(h, half, k) * data[i - half + k]; 0524 else 0525 result[i] += gsl_matrix_get(h, half, k) * nsl_smooth_pad_constant_rvalue; 0526 break; 0527 case nsl_smooth_pad_periodic: 0528 result[i] += gsl_matrix_get(h, half, k) * data[(i - half + k) % n]; 0529 case nsl_smooth_pad_none: 0530 break; 0531 } 0532 } 0533 } 0534 0535 gsl_matrix_free(h); 0536 0537 for (i = 0; i < n; i++) 0538 data[i] = result[i]; 0539 free(result); 0540 0541 return 0; 0542 } 0543 0544 int nsl_smooth_savgol_default(double* data, size_t n, size_t points, int order) { 0545 return nsl_smooth_savgol(data, n, points, order, nsl_smooth_pad_constant); 0546 }