GRASS GIS 8 Programmer's Manual  8.3.2(2024)-exported
fft.c
Go to the documentation of this file.
1 /**
2  * \file fft.c
3  *
4  * \brief Fast Fourier Transformation of Two Dimensional Satellite Data
5  * functions.
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or (at
10  * your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful, but
13  * WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15  * General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20  *
21  * \author GRASS GIS Development Team
22  *
23  * \date 2001-2006
24  */
25 
26 #include <grass/config.h>
27 
28 #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
29 
30 #ifdef HAVE_FFTW_H
31 #include <fftw.h>
32 #endif
33 
34 #ifdef HAVE_DFFTW_H
35 #include <dfftw.h>
36 #endif
37 
38 #ifdef HAVE_FFTW3_H
39 #include <fftw3.h>
40 #define c_re(c) ((c)[0])
41 #define c_im(c) ((c)[1])
42 #endif
43 
44 #include <stdlib.h>
45 #include <stdio.h>
46 #include <math.h>
47 #include <grass/gmath.h>
48 #include <grass/gis.h>
49 
50 /**
51  * \fn int fft2(int i_sign, double (*data)[2], int NN, int dimc, int dimr)
52  *
53  * \brief Fast Fourier Transform for two-dimensional array.
54  *
55  * Fast Fourier Transform for two-dimensional array.<br>
56  * <bNote:</b> If passing real data to fft() forward transform
57  * (especially when using fft() in a loop), explicitly (re-)initialize
58  * the imaginary part to zero (DATA[1][i] = 0.0). Returns 0.
59  *
60  * \param[in] i_sign Direction of transform -1 is normal, +1 is inverse
61  * \param[in,out] data Pointer to complex linear array in row major order
62  * containing data and result
63  * \param[in] NN Value of DATA dimension (dimc * dimr)
64  * \param[in] dimc Value of image column dimension (max power of 2)
65  * \param[in] dimr Value of image row dimension (max power of 2)
66  * \return int always returns 0
67  */
68 
69 int fft2(int i_sign, double (*data)[2], int NN, int dimc, int dimr)
70 {
71 #ifdef HAVE_FFTW3_H
72  fftw_plan plan;
73 #else
74  fftwnd_plan plan;
75 #endif
76  double norm;
77  int i;
78 
79  norm = 1.0 / sqrt(NN);
80 
81 #ifdef HAVE_FFTW3_H
82  plan = fftw_plan_dft_2d(dimr, dimc, data, data,
83  (i_sign < 0) ? FFTW_FORWARD : FFTW_BACKWARD,
84  FFTW_ESTIMATE);
85 
86  fftw_execute(plan);
87 
88  fftw_destroy_plan(plan);
89 #else
90  plan = fftw2d_create_plan(dimc, dimr,
91  (i_sign < 0) ? FFTW_FORWARD : FFTW_BACKWARD,
92  FFTW_ESTIMATE | FFTW_IN_PLACE);
93 
94  fftwnd_one(plan, data, data);
95 
96  fftwnd_destroy_plan(plan);
97 #endif
98 
99  for (i = 0; i < NN; i++) {
100  data[i][0] *= norm;
101  data[i][1] *= norm;
102  }
103 
104  return 0;
105 }
106 
107 /**
108  * \fn int fft(int i_sign, double *DATA[2], int NN, int dimc, int dimr)
109  *
110  * \brief Fast Fourier Transform for two-dimensional array.
111  *
112  * Fast Fourier Transform for two-dimensional array.<br>
113  * <bNote:</b> If passing real data to fft() forward transform
114  * (especially when using fft() in a loop), explicitly (re-)initialize
115  * the imaginary part to zero (DATA[1][i] = 0.0). Returns 0.
116  *
117  * \param[in] i_sign Direction of transform -1 is normal, +1 is inverse
118  * \param[in,out] DATA Pointer to complex linear array in row major order
119  * containing data and result
120  * \param[in] NN Value of DATA dimension (dimc * dimr)
121  * \param[in] dimc Value of image column dimension (max power of 2)
122  * \param[in] dimr Value of image row dimension (max power of 2)
123  * \return int always returns 0
124  */
125 
126 int fft(int i_sign, double *DATA[2], int NN, int dimc, int dimr)
127 {
128  fftw_complex *data;
129  int i;
130 
131  data = (fftw_complex *)G_malloc(NN * sizeof(fftw_complex));
132 
133  for (i = 0; i < NN; i++) {
134  c_re(data[i]) = DATA[0][i];
135  c_im(data[i]) = DATA[1][i];
136  }
137 
138  fft2(i_sign, data, NN, dimc, dimr);
139 
140  for (i = 0; i < NN; i++) {
141  DATA[0][i] = c_re(data[i]);
142  DATA[1][i] = c_im(data[i]);
143  }
144 
145  G_free(data);
146 
147  return 0;
148 }
149 
150 #endif /* HAVE_FFT */
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:150