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