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_dft.c
Go to the documentation of this file.
1 /***************************************************************************
2  File : nsl_dft.c
3  Project : LabPlot
4  Description : NSL discrete Fourier transform functions
5  --------------------------------------------------------------------
6  Copyright : (C) 2016 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_dft.h"
30 #include "nsl_common.h"
31 #include <gsl/gsl_fft_real.h>
32 #include <gsl/gsl_fft_halfcomplex.h>
33 #ifdef HAVE_FFTW3
34 #include <fftw3.h>
35 #endif
36 
37 const char* nsl_dft_result_type_name[] = {i18n("Magnitude"), i18n("Amplitude"), i18n("real part"), i18n("imaginary part"), i18n("Power"), i18n("Phase"),
38  i18n("Amplitude in dB"), i18n("normalized amplitude in dB"), i18n("Magnitude squared"), i18n("Amplitude squared"), i18n("raw")};
39 const char* nsl_dft_xscale_name[] = {i18n("Frequency"), i18n("Index"), i18n("Period")};
40 
41 int nsl_dft_transform_window(double data[], size_t stride, size_t n, int two_sided, nsl_dft_result_type type, nsl_sf_window_type window_type) {
42  /* apply window function */
43  if (window_type != nsl_sf_window_uniform)
44  nsl_sf_apply_window(data, n, window_type);
45 
46  /* transform */
47  int status = nsl_dft_transform(data, stride, n, two_sided, type);
48 
49  return status;
50 }
51 
52 int nsl_dft_transform(double data[], size_t stride, size_t n, int two_sided, nsl_dft_result_type type) {
53  if (n < 2) // we need at least 2 points
54  return 1;
55  size_t i;
56  double* result = (double*)malloc(2*n*sizeof(double));
57  size_t N = n/2; /* number of resulting data points */
58  if(two_sided)
59  N = n;
60 #ifdef HAVE_FFTW3
61  /* stride ignored */
62  (void)stride;
63 
64  fftw_plan plan = fftw_plan_dft_r2c_1d(n, data, (fftw_complex *) result, FFTW_ESTIMATE);
65  fftw_execute(plan);
66  fftw_destroy_plan(plan);
67 
68  /* 2. unpack data */
69  if(two_sided) {
70  for(i = 1; i < n-i; i++) {
71  result[2*(n - i)] = result[2*i];
72  result[2*(n - i)+1] = -result[2*i+1];
73  }
74  if (i == n - i) {
75  result[2*i] = result[2*(n-1)];
76  result[2*i+1] = 0;
77  }
78  }
79 #else
80  /* 1. transform */
81  gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(n);
82  gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(n);
83 
84  gsl_fft_real_transform(data, stride, n, real, work);
85  gsl_fft_real_wavetable_free(real);
86  gsl_fft_real_workspace_free(work);
87 
88  /* 2. unpack data */
89  gsl_fft_halfcomplex_unpack(data, result, stride, n);
90 #endif
91 
92  /* 3. write result */
93  switch(type) {
95  for (i = 0; i < N; i++)
96  data[i] = sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]));
97  break;
99  for (i = 0; i < N; i++) {
100  data[i] = sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/n;
101  if(i>0)
102  data[i] *= 2.;
103  }
104  break;
105  case nsl_dft_result_real:
106  for (i = 0; i < N; i++)
107  data[i] = result[2*i];
108  break;
109  case nsl_dft_result_imag:
110  for (i = 0; i < N; i++)
111  data[i] = result[2*i+1];
112  break;
114  for (i = 0; i < N; i++) {
115  data[i] = (gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/n;
116  if (i > 0)
117  data[i] *= 2.;
118  }
119  break;
121  for (i = 0; i < N; i++)
122  data[i] = -atan2(result[2*i+1],result[2*i]);
123  break;
124  case nsl_dft_result_dB:
125  for (i = 0; i < N; i++) {
126  if (i == 0)
127  data[i] = 20.*log10(sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/(double)n);
128  else
129  data[i] = 20.*log10(2.*sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/(double)n);
130  }
131  break;
132  case nsl_dft_result_normdB: {
133  double maxdB=0;
134  for (i = 0; i < N; i++) {
135  if (i == 0)
136  data[i] = 20.*log10(sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/(double)n);
137  else
138  data[i] = 20.*log10(2.*sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/(double)n);
139 
140  if (i == 0 || maxdB < data[i])
141  maxdB = data[i];
142  }
143  for (i = 0; i < N; i++)
144  data[i] -= maxdB;
145  break;
146  }
148  for (i = 0; i < N; i++)
149  data[i] = gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]);
150  break;
152  for (i = 0; i < N; i++) {
153  data[i] = (gsl_pow_2(result[2*i]) + gsl_pow_2(result[2*i+1]))/gsl_pow_2((double)n);
154  if (i > 0)
155  data[i] *= 4.;
156  }
157  break;
158  case nsl_dft_result_raw:
159  break;
160  }
161 
162  free(result);
163 
164  return 0;
165 }
166 
#define i18n(m)
Definition: nsl_common.h:38
const char * nsl_dft_xscale_name[]
Definition: nsl_dft.c:39
int nsl_dft_transform_window(double data[], size_t stride, size_t n, int two_sided, nsl_dft_result_type type, nsl_sf_window_type window_type)
Definition: nsl_dft.c:41
const char * nsl_dft_result_type_name[]
Definition: nsl_dft.c:37
int nsl_dft_transform(double data[], size_t stride, size_t n, int two_sided, nsl_dft_result_type type)
Definition: nsl_dft.c:52
nsl_dft_result_type
Definition: nsl_dft.h:50
@ nsl_dft_result_squareamplitude
Definition: nsl_dft.h:51
@ nsl_dft_result_raw
Definition: nsl_dft.h:52
@ nsl_dft_result_real
Definition: nsl_dft.h:50
@ nsl_dft_result_phase
Definition: nsl_dft.h:51
@ nsl_dft_result_normdB
Definition: nsl_dft.h:51
@ nsl_dft_result_imag
Definition: nsl_dft.h:50
@ nsl_dft_result_dB
Definition: nsl_dft.h:51
@ nsl_dft_result_magnitude
Definition: nsl_dft.h:50
@ nsl_dft_result_power
Definition: nsl_dft.h:50
@ nsl_dft_result_amplitude
Definition: nsl_dft.h:50
@ nsl_dft_result_squaremagnitude
Definition: nsl_dft.h:51
int nsl_sf_apply_window(double data[], size_t N, nsl_sf_window_type type)
Definition: nsl_sf_window.c:38
nsl_sf_window_type
Definition: nsl_sf_window.h:35
@ nsl_sf_window_uniform
Definition: nsl_sf_window.h:35