labplot  2.8.2
About: LabPlot is an application for plotting and analysis of 2D and 3D functions and data. It is a complete rewrite of LabPlot1 and lacks in the first release a lot of features available in the predecessor. On the other hand, the GUI and the usability is more superior.
  Fossies Dox: labplot-2.8.2.tar.gz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

nsl_corr.c
Go to the documentation of this file.
1 /***************************************************************************
2  File : nsl_corr.c
3  Project : LabPlot
4  Description : NSL discrete correlation functions
5  --------------------------------------------------------------------
6  Copyright : (C) 2018 by Stefan Gerlach (stefan.gerlach@uni.kn)
7 
8  ***************************************************************************/
9 
10 /***************************************************************************
11  * *
12  * This program is free software; you can redistribute it and/or modify *
13  * it under the terms of the GNU General Public License as published by *
14  * the Free Software Foundation; either version 2 of the License, or *
15  * (at your option) any later version. *
16  * *
17  * This program is distributed in the hope that it will be useful, *
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
20  * GNU General Public License for more details. *
21  * *
22  * You should have received a copy of the GNU General Public License *
23  * along with this program; if not, write to the Free Software *
24  * Foundation, Inc., 51 Franklin Street, Fifth Floor, *
25  * Boston, MA 02110-1301 USA *
26  * *
27  ***************************************************************************/
28 
29 #include "nsl_corr.h"
30 #include "nsl_common.h"
31 #include <gsl/gsl_cblas.h>
32 #include <gsl/gsl_fft_halfcomplex.h>
33 #ifdef HAVE_FFTW3
34 #include <fftw3.h>
35 #endif
36 
37 const char* nsl_corr_type_name[] = {i18n("linear (zero-padded)"), i18n("circular")};
38 const char* nsl_corr_norm_name[] = {i18n("none"), i18n("biased"), i18n("unbiased"), i18n("coeff")};
39 
40 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[]) {
41  return nsl_corr_fft_type(s, n, r, m, type, normalize, out);
42 }
43 
44 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[]) {
45  size_t i, size, N = GSL_MAX(n, m), maxlag = N - 1;
47  size = maxlag + N;
48  else // circular
49  size = N;
50 
51  size_t oldsize = size;
52 #ifdef HAVE_FFTW3
53  // already zero-pad here for FFT method and FFTW r2c
54  size = 2*(size/2+1);
55 #endif
56 
57  // zero-padded arrays
58  double *stmp = (double*)malloc(size*sizeof(double));
59  if (stmp == NULL) {
60  printf("nsl_corr_fft_type(): ERROR allocating memory for 'stmp'!\n");
61  return -1;
62  }
63 
64  double *rtmp = (double*)malloc(size*sizeof(double));
65  if (rtmp == NULL) {
66  free(stmp);
67  printf("nsl_corr_fft_type(): ERROR allocating memory for 'rtmp'!\n");
68  return -1;
69  }
70 
71  if (type == nsl_corr_type_linear) {
72  for (i = 0; i < maxlag; i++)
73  stmp[i] = 0.;
74  for (i = 0; i < n; i++)
75  stmp[i+maxlag] = s[i];
76  for (i = n+maxlag; i < size; i++)
77  stmp[i] = 0;
78  for (i = 0; i < m; i++)
79  rtmp[i] = r[i];
80  for (i = m; i < size; i++)
81  rtmp[i] = 0;
82  } else { // circular
83  for (i = 0; i < n; i++)
84  stmp[i] = s[i];
85  for (i = n; i < N; i++)
86  stmp[i] = 0.;
87  for (i = 0; i < m; i++)
88  rtmp[i] = r[i];
89  for (i = m; i < N; i++)
90  rtmp[i] = 0.;
91  }
92 
93  int status;
94 #ifdef HAVE_FFTW3 // already wraps output
95  status = nsl_corr_fft_FFTW(stmp, rtmp, oldsize, out);
96 #else // GSL
97  status = nsl_corr_fft_GSL(stmp, rtmp, size, out);
98 #endif
99 
100  free(stmp);
101  free(rtmp);
102 
103  /* normalization */
104  switch (normalize) {
105  case nsl_corr_norm_none:
106  break;
108  for (i = 0; i < oldsize; i++)
109  out[i] = out[i]/N;
110  break;
112  for (i = 0; i < oldsize; i++) {
113  size_t norm = i < oldsize/2 ? i+1 : oldsize - i;
114  out[i] = out[i]/norm;
115  }
116  break;
117  case nsl_corr_norm_coeff: {
118  double snorm = cblas_dnrm2((int)n, s, 1);
119  double rnorm = cblas_dnrm2((int)m, r, 1);
120  for (i = 0; i < oldsize; i++)
121  out[i] = out[i]/snorm/rnorm;
122  break;
123  }
124  }
125 
126  // reverse array for circular type
127  if (type == nsl_corr_type_circular) {
128  for (i = 0; i < N/2; i++) {
129  double tmp = out[i];
130  out[i] = out[N - i - 1];
131  out[N -i - 1] = tmp;
132  }
133  }
134 
135  return status;
136 }
137 
138 #ifdef HAVE_FFTW3
139 int nsl_corr_fft_FFTW(double s[], double r[], size_t n, double out[]) {
140  if (n <= 0)
141  return -1;
142 
143  const size_t size = 2*(n/2+1);
144  double* in = (double*)malloc(size*sizeof(double));
145  fftw_plan rpf = fftw_plan_dft_r2c_1d((int)n, in, (fftw_complex*)in, FFTW_ESTIMATE);
146 
147  fftw_execute_dft_r2c(rpf, s, (fftw_complex*)s);
148  fftw_execute_dft_r2c(rpf, r, (fftw_complex*)r);
149  fftw_destroy_plan(rpf);
150 
151  size_t i;
152 
153  // multiply
154  for (i = 0; i < size; i += 2) {
155  double re = s[i]*r[i] + s[i+1]*r[i+1];
156  double im = s[i+1]*r[i] - s[i]*r[i+1];
157 
158  s[i] = re;
159  s[i+1] = im;
160  }
161 
162  // back transform
163  double* o = (double*)malloc(size*sizeof(double));
164  fftw_plan rpb = fftw_plan_dft_c2r_1d((int)n, (fftw_complex*)o, o, FFTW_ESTIMATE);
165  fftw_execute_dft_c2r(rpb, (fftw_complex*)s, s);
166  fftw_destroy_plan(rpb);
167 
168  for (i = 0; i < n; i++)
169  out[i] = s[i]/n;
170 
171  free(in);
172  free(o);
173  return 0;
174 }
175 #endif
176 
177 int nsl_corr_fft_GSL(double s[], double r[], size_t n, double out[]) {
178  gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(n);
179  gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(n);
180 
181  /* FFT s and r */
182  gsl_fft_real_transform(s, 1, n, real, work);
183  gsl_fft_real_transform(r, 1, n, real, work);
184  gsl_fft_real_wavetable_free(real);
185 
186  size_t i;
187  /* calculate halfcomplex product */
188  out[0] = s[0]*r[0];
189  for (i = 1; i < n; i++) {
190  if (i%2) { /* Re */
191  out[i] = s[i]*r[i];
192  if (i < n-1) /* when n is even last value is real */
193  out[i] += s[i+1]*r[i+1];
194  } else /* Im */
195  out[i] = s[i]*r[i-1] - s[i-1]*r[i];
196  }
197 
198  /* back transform */
199  gsl_fft_halfcomplex_wavetable *hc = gsl_fft_halfcomplex_wavetable_alloc(n);
200  gsl_fft_halfcomplex_inverse(out, 1, n, hc, work);
201  gsl_fft_halfcomplex_wavetable_free(hc);
202  gsl_fft_real_workspace_free(work);
203 
204  return 0;
205 }
#define i18n(m)
Definition: nsl_common.h:38
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[])
Definition: nsl_corr.c:40
int nsl_corr_fft_GSL(double s[], double r[], size_t n, double out[])
Definition: nsl_corr.c:177
const char * nsl_corr_type_name[]
Definition: nsl_corr.c:37
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[])
Definition: nsl_corr.c:44
const char * nsl_corr_norm_name[]
Definition: nsl_corr.c:38
nsl_corr_type_type
Definition: nsl_corr.h:36
@ nsl_corr_type_circular
Definition: nsl_corr.h:36
@ nsl_corr_type_linear
Definition: nsl_corr.h:36
nsl_corr_norm_type
Definition: nsl_corr.h:45
@ nsl_corr_norm_coeff
Definition: nsl_corr.h:45
@ nsl_corr_norm_biased
Definition: nsl_corr.h:45
@ nsl_corr_norm_none
Definition: nsl_corr.h:45
@ nsl_corr_norm_unbiased
Definition: nsl_corr.h:45