File indexing completed on 2024-05-12 03:47:50
0001 /* 0002 File : nsl_corr.c 0003 Project : LabPlot 0004 Description : NSL discrete correlation functions 0005 -------------------------------------------------------------------- 0006 SPDX-FileCopyrightText: 2018 Stefan Gerlach <stefan.gerlach@uni.kn> 0007 SPDX-License-Identifier: GPL-2.0-or-later 0008 */ 0009 0010 #include "nsl_corr.h" 0011 #include "nsl_common.h" 0012 #include <gsl/gsl_cblas.h> 0013 #include <gsl/gsl_fft_halfcomplex.h> 0014 #ifdef HAVE_FFTW3 0015 #include <fftw3.h> 0016 #endif 0017 0018 const char* nsl_corr_type_name[] = {i18n("Linear (Zero-padded)"), i18n("Circular")}; 0019 const char* nsl_corr_norm_name[] = {i18n("None"), i18n("Biased"), i18n("Unbiased"), i18n("Coeff")}; 0020 0021 int nsl_corr_correlation(double s[], size_t n, double r[], size_t m, nsl_corr_type_type type, nsl_corr_norm_type normalize, double out[]) { 0022 return nsl_corr_fft_type(s, n, r, m, type, normalize, out); 0023 } 0024 0025 int nsl_corr_fft_type(double s[], size_t n, double r[], size_t m, nsl_corr_type_type type, nsl_corr_norm_type normalize, double out[]) { 0026 size_t i, size, N = GSL_MAX(n, m), maxlag = N - 1; 0027 if (type == nsl_corr_type_linear) 0028 size = maxlag + N; 0029 else // circular 0030 size = N; 0031 0032 size_t oldsize = size; 0033 #ifdef HAVE_FFTW3 0034 // already zero-pad here for FFT method and FFTW r2c 0035 size = 2 * (size / 2 + 1); 0036 #endif 0037 0038 // zero-padded arrays 0039 double* stmp = (double*)malloc(size * sizeof(double)); 0040 if (stmp == NULL) { 0041 printf("nsl_corr_fft_type(): ERROR allocating memory for 'stmp'!\n"); 0042 return -1; 0043 } 0044 0045 double* rtmp = (double*)malloc(size * sizeof(double)); 0046 if (rtmp == NULL) { 0047 free(stmp); 0048 printf("nsl_corr_fft_type(): ERROR allocating memory for 'rtmp'!\n"); 0049 return -1; 0050 } 0051 0052 if (type == nsl_corr_type_linear) { 0053 for (i = 0; i < maxlag; i++) 0054 stmp[i] = 0.; 0055 for (i = 0; i < n; i++) 0056 stmp[i + maxlag] = s[i]; 0057 for (i = n + maxlag; i < size; i++) 0058 stmp[i] = 0; 0059 for (i = 0; i < m; i++) 0060 rtmp[i] = r[i]; 0061 for (i = m; i < size; i++) 0062 rtmp[i] = 0; 0063 } else { // circular 0064 for (i = 0; i < n; i++) 0065 stmp[i] = s[i]; 0066 for (i = n; i < N; i++) 0067 stmp[i] = 0.; 0068 for (i = 0; i < m; i++) 0069 rtmp[i] = r[i]; 0070 for (i = m; i < N; i++) 0071 rtmp[i] = 0.; 0072 } 0073 0074 int status; 0075 #ifdef HAVE_FFTW3 // already wraps output 0076 status = nsl_corr_fft_FFTW(stmp, rtmp, oldsize, out); 0077 #else // GSL 0078 status = nsl_corr_fft_GSL(stmp, rtmp, size, out); 0079 #endif 0080 0081 free(stmp); 0082 free(rtmp); 0083 0084 /* normalization */ 0085 switch (normalize) { 0086 case nsl_corr_norm_none: 0087 break; 0088 case nsl_corr_norm_biased: 0089 for (i = 0; i < oldsize; i++) 0090 out[i] = out[i] / N; 0091 break; 0092 case nsl_corr_norm_unbiased: 0093 for (i = 0; i < oldsize; i++) { 0094 size_t norm = i < oldsize / 2 ? i + 1 : oldsize - i; 0095 out[i] = out[i] / norm; 0096 } 0097 break; 0098 case nsl_corr_norm_coeff: { 0099 double snorm = cblas_dnrm2((int)n, s, 1); 0100 double rnorm = cblas_dnrm2((int)m, r, 1); 0101 for (i = 0; i < oldsize; i++) 0102 out[i] = out[i] / snorm / rnorm; 0103 break; 0104 } 0105 } 0106 0107 // reverse array for circular type 0108 if (type == nsl_corr_type_circular) { 0109 for (i = 0; i < N / 2; i++) { 0110 double tmp = out[i]; 0111 out[i] = out[N - i - 1]; 0112 out[N - i - 1] = tmp; 0113 } 0114 } 0115 0116 return status; 0117 } 0118 0119 #ifdef HAVE_FFTW3 0120 int nsl_corr_fft_FFTW(double s[], double r[], size_t n, double out[]) { 0121 if (n == 0) 0122 return -1; 0123 0124 const size_t size = 2 * (n / 2 + 1); 0125 double* in = (double*)malloc(size * sizeof(double)); 0126 fftw_plan rpf = fftw_plan_dft_r2c_1d((int)n, in, (fftw_complex*)in, FFTW_ESTIMATE); 0127 0128 fftw_execute_dft_r2c(rpf, s, (fftw_complex*)s); 0129 fftw_execute_dft_r2c(rpf, r, (fftw_complex*)r); 0130 fftw_destroy_plan(rpf); 0131 0132 size_t i; 0133 0134 // multiply 0135 for (i = 0; i < size; i += 2) { 0136 double re = s[i] * r[i] + s[i + 1] * r[i + 1]; 0137 double im = s[i + 1] * r[i] - s[i] * r[i + 1]; 0138 0139 s[i] = re; 0140 s[i + 1] = im; 0141 } 0142 0143 // back transform 0144 double* o = (double*)malloc(size * sizeof(double)); 0145 fftw_plan rpb = fftw_plan_dft_c2r_1d((int)n, (fftw_complex*)o, o, FFTW_ESTIMATE); 0146 fftw_execute_dft_c2r(rpb, (fftw_complex*)s, s); 0147 fftw_destroy_plan(rpb); 0148 0149 for (i = 0; i < n; i++) 0150 out[i] = s[i] / n; 0151 0152 free(in); 0153 free(o); 0154 return 0; 0155 } 0156 #endif 0157 0158 int nsl_corr_fft_GSL(double s[], double r[], size_t n, double out[]) { 0159 gsl_fft_real_workspace* work = gsl_fft_real_workspace_alloc(n); 0160 gsl_fft_real_wavetable* real = gsl_fft_real_wavetable_alloc(n); 0161 0162 /* FFT s and r */ 0163 gsl_fft_real_transform(s, 1, n, real, work); 0164 gsl_fft_real_transform(r, 1, n, real, work); 0165 gsl_fft_real_wavetable_free(real); 0166 0167 size_t i; 0168 /* calculate halfcomplex product */ 0169 out[0] = s[0] * r[0]; 0170 for (i = 1; i < n; i++) { 0171 if (i % 2) { /* Re */ 0172 out[i] = s[i] * r[i]; 0173 if (i < n - 1) /* when n is even last value is real */ 0174 out[i] += s[i + 1] * r[i + 1]; 0175 } else /* Im */ 0176 out[i] = s[i] * r[i - 1] - s[i - 1] * r[i]; 0177 } 0178 0179 /* back transform */ 0180 gsl_fft_halfcomplex_wavetable* hc = gsl_fft_halfcomplex_wavetable_alloc(n); 0181 gsl_fft_halfcomplex_inverse(out, 1, n, hc, work); 0182 gsl_fft_halfcomplex_wavetable_free(hc); 0183 gsl_fft_real_workspace_free(work); 0184 0185 return 0; 0186 }