File indexing completed on 2024-05-19 15:01:44
0001 /*************************************************************************** 0002 File : nsl_conv.c 0003 Project : LabPlot 0004 Description : NSL discrete convolution functions 0005 -------------------------------------------------------------------- 0006 Copyright : (C) 2018 by Stefan Gerlach (stefan.gerlach@uni.kn) 0007 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_conv.h" 0030 #include "nsl_common.h" 0031 #include <gsl/gsl_cblas.h> 0032 #include <gsl/gsl_fft_halfcomplex.h> 0033 #ifdef HAVE_FFTW3 0034 #include <fftw3.h> 0035 #endif 0036 #include "backend/nsl/nsl_stats.h" 0037 0038 const char* nsl_conv_direction_name[] = {i18n("forward (convolution)"), i18n("backward (deconvolution)")}; 0039 const char* nsl_conv_type_name[] = {i18n("linear (zero-padded)"), i18n("circular")}; 0040 const char* nsl_conv_method_name[] = {i18n("auto"), i18n("direct"), i18n("FFT")}; 0041 const char* nsl_conv_norm_name[] = {i18n("none"), i18n("sum"), i18n("Euclidean")}; 0042 const char* nsl_conv_wrap_name[] = {i18n("none"), i18n("maximum"), i18n("center (acausal)")}; 0043 const char* nsl_conv_kernel_name[] = {i18n("sliding average"), i18n("triangular smooth"), i18n("pseudo-Gaussian smooth"), i18n("first derivative"), i18n("smooth first derivative"), 0044 i18n("second derivative"), i18n("third derivative"), i18n("fourth derivative"), i18n("Gaussian"), i18n("Lorentzian") }; 0045 0046 int nsl_conv_standard_kernel(double k[], size_t n, nsl_conv_kernel_type type) { 0047 size_t i; 0048 0049 int validn = 1; 0050 switch (type) { 0051 case nsl_conv_kernel_avg: 0052 for (i = 0; i < n; i++) 0053 k[i] = 1.; 0054 break; 0055 case nsl_conv_kernel_smooth_triangle: 0056 for (i = 0; i < n/2; i++) 0057 k[i] = i + 1.; 0058 for (i = n/2; i < n; i++) 0059 k[i] = (double)(n - i); 0060 break; 0061 case nsl_conv_kernel_smooth_gaussian: 0062 switch (n) { 0063 case 5: 0064 k[0]=1;k[1]=4;k[2]=6;k[3]=4;k[4]=1; 0065 break; 0066 case 7: 0067 k[0]=1;k[1]=4;k[2]=8;k[3]=10;k[4]=8;k[5]=4;k[6]=1; 0068 break; 0069 case 9: 0070 k[0]=1;k[1]=4;k[2]=9;k[3]=14;k[4]=17;k[5]=14;k[6]=9;k[7]=4;k[8]=1; 0071 break; 0072 default: 0073 validn = 0; 0074 } 0075 break; 0076 case nsl_conv_kernel_first_derivative: 0077 if (n == 2) { 0078 k[0]=-1;k[1]=1; 0079 } else 0080 validn = 0; 0081 break; 0082 case nsl_conv_kernel_smooth_first_derivative: 0083 if (n%2) { 0084 for (i = 0; i < n; i++) 0085 k[i] = (int)(i - n/2); 0086 } else 0087 validn = 0; 0088 break; 0089 case nsl_conv_kernel_second_derivative: 0090 if (n == 3) { 0091 k[0]=1;k[1]=-2;k[2]=1; 0092 } else 0093 validn = 0; 0094 break; 0095 case nsl_conv_kernel_third_derivative: 0096 if (n == 4) { 0097 k[0]=1;k[1]=-3;k[2]=3;k[3]=-1; 0098 } else 0099 validn = 0; 0100 break; 0101 case nsl_conv_kernel_fourth_derivative: 0102 if (n == 5) { 0103 k[0]=1;k[1]=-4;k[2]=6;k[3]=-4;k[4]=1; 0104 } else 0105 validn = 0; 0106 break; 0107 case nsl_conv_kernel_gaussian: { 0108 double s = n/5.; /* relative width */ 0109 for (i = 0;i < n; i++) { 0110 double x = i - (n-1.)/2.; 0111 k[i] = M_SQRT1_2/M_SQRTPI/s * exp(-x*x/2./s/s); 0112 } 0113 break; 0114 } 0115 case nsl_conv_kernel_lorentzian: { 0116 double s = n/5.; /* relative width */ 0117 for (i = 0;i < n; i++) { 0118 double x = i - (n-1.)/2.; 0119 k[i] = s/M_PI/(s*s + x*x); 0120 } 0121 break; 0122 } 0123 } 0124 0125 if (!validn) { 0126 printf("ERROR: kernel size %lu not supported for kernel %s\n", (unsigned long)n, nsl_conv_kernel_name[type]); 0127 return -1; 0128 } 0129 0130 /* debug */ 0131 printf("["); 0132 for (i = 0; i < n; i++) 0133 printf("%g ", k[i]); 0134 puts("]"); 0135 0136 return 0; 0137 } 0138 0139 int nsl_conv_convolution_direction(double s[], size_t n, double r[], size_t m, nsl_conv_direction_type dir, nsl_conv_type_type type, nsl_conv_method_type method, 0140 nsl_conv_norm_type normalize, nsl_conv_wrap_type wrap, double out[]) { 0141 if (dir == nsl_conv_direction_forward) 0142 return nsl_conv_convolution(s, n, r, m, type, method, normalize, wrap, out); 0143 else 0144 return nsl_conv_deconvolution(s, n, r, m, type, normalize, wrap, out); 0145 } 0146 0147 int nsl_conv_convolution(double s[], size_t n, double r[], size_t m, nsl_conv_type_type type, nsl_conv_method_type method, nsl_conv_norm_type normalize, nsl_conv_wrap_type wrap, double out[]) { 0148 if (method == nsl_conv_method_direct || (method == nsl_conv_method_auto && GSL_MAX_INT(n, m) <= NSL_CONV_METHOD_BORDER)) { 0149 if (type == nsl_conv_type_linear) 0150 return nsl_conv_linear_direct(s, n, r, m, normalize, wrap, out); 0151 else if (type == nsl_conv_type_circular) 0152 return nsl_conv_circular_direct(s, n, r, m, normalize, wrap, out); 0153 } else { 0154 return nsl_conv_fft_type(s, n, r, m, nsl_conv_direction_forward, type, normalize, wrap, out); 0155 } 0156 0157 return 0; 0158 } 0159 0160 int nsl_conv_deconvolution(double s[], size_t n, double r[], size_t m, nsl_conv_type_type type, nsl_conv_norm_type normalize, nsl_conv_wrap_type wrap, double out[]) { 0161 /* only supported by FFT method */ 0162 return nsl_conv_fft_type(s, n, r, m, nsl_conv_direction_backward, type, normalize, wrap, out); 0163 } 0164 0165 int nsl_conv_linear_direct(double s[], size_t n, double r[], size_t m, nsl_conv_norm_type normalize, nsl_conv_wrap_type wrap, double out[]) { 0166 size_t i, j, size = n + m - 1, wi = 0; 0167 double norm = 1; 0168 if (normalize == nsl_conv_norm_euclidean) { 0169 if ((norm = cblas_dnrm2((int)m, r, 1)) == 0) 0170 norm = 1.; 0171 } else if (normalize == nsl_conv_norm_sum) { 0172 if ((norm = cblas_dasum((int)m, r, 1)) == 0) 0173 norm = 1.; 0174 } 0175 0176 if (wrap == nsl_conv_wrap_max) 0177 nsl_stats_maximum(r, m, &wi); 0178 else if (wrap == nsl_conv_wrap_center) 0179 wi = m/2; 0180 0181 for (j = 0; j < size; j++) { 0182 int index; // can be negative 0183 double res = 0; 0184 for (i = 0; i < n; i++) { 0185 index = (int)(j - i); 0186 if (index >= 0 && index < (int)m) 0187 res += s[i] * r[index]/norm; 0188 } 0189 index = (int)(j - wi); 0190 if (index < 0) 0191 index += (int)size; 0192 out[index] = res; 0193 } 0194 0195 return 0; 0196 } 0197 0198 int nsl_conv_circular_direct(double s[], size_t n, double r[], size_t m, nsl_conv_norm_type normalize, nsl_conv_wrap_type wrap, double out[]) { 0199 size_t i, j, size = GSL_MAX(n,m), wi = 0; 0200 double norm = 1; 0201 if (normalize == nsl_conv_norm_euclidean) { 0202 if ((norm = cblas_dnrm2((int)m, r, 1)) == 0) 0203 norm = 1.; 0204 } else if (normalize == nsl_conv_norm_sum) { 0205 if ((norm = cblas_dasum((int)m, r, 1)) == 0) 0206 norm = 1.; 0207 } 0208 0209 if (wrap == nsl_conv_wrap_max) 0210 nsl_stats_maximum(r, m, &wi); 0211 else if (wrap == nsl_conv_wrap_center) 0212 wi = m/2; 0213 0214 for (j = 0; j < size; j++) { 0215 int index; // can be negative 0216 double res = 0; 0217 for (i = 0; i < n; i++) { 0218 index = (int)(j - i); 0219 if (index < 0) 0220 index += (int)size; 0221 if (index < (int)m) 0222 res += s[i]*r[index]/norm; 0223 } 0224 index = (int)(j - wi); 0225 if (index < 0) 0226 index += (int)size; 0227 out[index] = res; 0228 } 0229 0230 return 0; 0231 } 0232 0233 int nsl_conv_fft_type(double s[], size_t n, double r[], size_t m, nsl_conv_direction_type dir, nsl_conv_type_type type, nsl_conv_norm_type normalize, nsl_conv_wrap_type wrap, double out[]) { 0234 size_t i, size, wi = 0; 0235 if (type == nsl_conv_type_linear) 0236 size = n + m - 1; 0237 else // circular 0238 size = GSL_MAX(n, m); 0239 0240 double norm = 1.; 0241 if (normalize == nsl_conv_norm_euclidean) { 0242 if ((norm = cblas_dnrm2((int)m, r, 1)) == 0) 0243 norm = 1.; 0244 } else if (normalize == nsl_conv_norm_sum) { 0245 if ((norm = cblas_dasum((int)m, r, 1)) == 0) 0246 norm = 1.; 0247 } 0248 0249 if (wrap == nsl_conv_wrap_max) 0250 nsl_stats_maximum(r, m, &wi); 0251 else if (wrap == nsl_conv_wrap_center) 0252 wi = m/2; 0253 0254 #ifdef HAVE_FFTW3 0255 // already zero-pad here for FFT method and FFTW r2c 0256 size_t oldsize = size; 0257 size = 2*(size/2+1); 0258 #endif 0259 0260 // zero-padded arrays 0261 double *stmp = (double*)malloc(size*sizeof(double)); 0262 if (stmp == NULL) { 0263 printf("nsl_conv_fft_type(): ERROR allocating memory for 'stmp'!\n"); 0264 return -1; 0265 } 0266 0267 double *rtmp = (double*)malloc(size*sizeof(double)); 0268 if (rtmp == NULL) { 0269 free(stmp); 0270 printf("nsl_corr_fft_type(): ERROR allocating memory for 'rtmp'!\n"); 0271 return -1; 0272 } 0273 0274 for (i = 0; i < n; i++) 0275 stmp[i] = s[i]; 0276 for (i = n; i < size; i++) 0277 stmp[i] = 0; 0278 for (i = 0; i < m; i++) { 0279 rtmp[i] = r[i]/norm; /* norm response */ 0280 } 0281 for (i = m; i < size; i++) 0282 rtmp[i] = 0; 0283 0284 int status; 0285 #ifdef HAVE_FFTW3 // already wraps output 0286 status = nsl_conv_fft_FFTW(stmp, rtmp, oldsize, dir, wi, out); 0287 #else // GSL 0288 status = nsl_conv_fft_GSL(stmp, rtmp, size, dir, out); 0289 for (i = 0; i < size; i++) { // wrap output 0290 size_t index = (i + wi) % size; 0291 stmp[i] = out[index]; 0292 } 0293 memcpy(out, stmp, size * sizeof(double)); 0294 #endif 0295 0296 free(stmp); 0297 free(rtmp); 0298 0299 return status; 0300 } 0301 0302 #ifdef HAVE_FFTW3 0303 int nsl_conv_fft_FFTW(double s[], double r[], size_t n, nsl_conv_direction_type dir, size_t wi, double out[]) { 0304 size_t i; 0305 const size_t size = 2*(n/2+1); 0306 double* in = (double*)malloc(size*sizeof(double)); 0307 fftw_plan rpf = fftw_plan_dft_r2c_1d((int)n, in, (fftw_complex*)in, FFTW_ESTIMATE); 0308 0309 fftw_execute_dft_r2c(rpf, s, (fftw_complex*)s); 0310 fftw_execute_dft_r2c(rpf, r, (fftw_complex*)r); 0311 fftw_destroy_plan(rpf); 0312 free(in); 0313 0314 // multiply/divide 0315 if (dir == nsl_conv_direction_forward) { 0316 for (i = 0; i < size; i += 2) { 0317 double re = s[i]*r[i] - s[i+1]*r[i+1]; 0318 double im = s[i]*r[i+1] + s[i+1]*r[i]; 0319 0320 s[i] = re; 0321 s[i+1] = im; 0322 } 0323 } else { 0324 for (i = 0; i < size; i += 2) { 0325 double norm = r[i]*r[i] + r[i+1]*r[i+1]; 0326 if (norm < DBL_MIN) 0327 norm = 1.; 0328 double re = (s[i]*r[i] + s[i+1]*r[i+1])/norm; 0329 double im = (s[i+1]*r[i] - s[i]*r[i+1])/norm; 0330 0331 s[i] = re; 0332 s[i+1] = im; 0333 } 0334 } 0335 0336 // back transform 0337 double* o = (double*)malloc(size*sizeof(double)); 0338 fftw_plan rpb = fftw_plan_dft_c2r_1d((int)n, (fftw_complex*)o, o, FFTW_ESTIMATE); 0339 0340 fftw_execute_dft_c2r(rpb, (fftw_complex*)s, s); 0341 fftw_destroy_plan(rpb); 0342 0343 for (i = 0; i < n; i++) { 0344 size_t index = (i + wi) % n; 0345 out[i] = s[index]/n; 0346 } 0347 free(o); 0348 0349 return 0; 0350 } 0351 #endif 0352 0353 int nsl_conv_fft_GSL(double s[], double r[], size_t n, nsl_conv_direction_type dir, double out[]) { 0354 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(n); 0355 gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(n); 0356 0357 /* FFT s and r */ 0358 gsl_fft_real_transform(s, 1, n, real, work); 0359 gsl_fft_real_transform(r, 1, n, real, work); 0360 gsl_fft_real_wavetable_free(real); 0361 0362 size_t i; 0363 /* calculate halfcomplex product/quotient depending on direction */ 0364 if (dir == nsl_conv_direction_forward) { 0365 out[0] = s[0]*r[0]; 0366 for (i = 1; i < n; i++) { 0367 if (i%2) { /* Re */ 0368 out[i] = s[i]*r[i]; 0369 if (i < n-1) /* when n is even last value is real */ 0370 out[i] -= s[i+1]*r[i+1]; 0371 } else /* Im */ 0372 out[i] = s[i-1]*r[i] + s[i]*r[i-1]; 0373 } 0374 } else { 0375 out[0] = s[0]/r[0]; 0376 for (i = 1; i < n; i++) { 0377 if (i%2) { /* Re */ 0378 if (i == n-1) /* when n is even last value is real */ 0379 out[i] = s[i]/r[i]; 0380 else { 0381 double norm = r[i]*r[i] + r[i+1]*r[i+1]; 0382 if (norm < DBL_MIN) 0383 norm = 1.; 0384 out[i] = (s[i]*r[i] + s[i+1]*r[i+1])/norm; 0385 } 0386 } else { /* Im */ 0387 double norm = r[i-1]*r[i-1] + r[i]*r[i]; 0388 if (norm < DBL_MIN) 0389 norm = 1.; 0390 out[i] = (s[i]*r[i-1] - s[i-1]*r[i])/norm; 0391 } 0392 } 0393 } 0394 0395 /* back transform */ 0396 gsl_fft_halfcomplex_wavetable *hc = gsl_fft_halfcomplex_wavetable_alloc(n); 0397 gsl_fft_halfcomplex_inverse(out, 1, n, hc, work); 0398 gsl_fft_halfcomplex_wavetable_free(hc); 0399 gsl_fft_real_workspace_free(work); 0400 0401 return 0; 0402 } 0403