GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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_FFTW3_H) || defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H)
29
30#if defined(HAVE_FFTW3_H)
31#include <fftw3.h>
32#define c_re(c) ((c)[0])
33#define c_im(c) ((c)[1])
34#elif defined(HAVE_FFTW_H)
35#include <fftw.h>
36#elif defined(HAVE_DFFTW_H)
37#include <dfftw.h>
38#endif
39
40#include <stdlib.h>
41#include <stdio.h>
42#include <math.h>
43#include <grass/gmath.h>
44#include <grass/gis.h>
45
46/**
47 * \fn int fft2(int i_sign, double (*data)[2], int NN, int dimc, int dimr)
48 *
49 * \brief Fast Fourier Transform for two-dimensional array.
50 *
51 * Fast Fourier Transform for two-dimensional array.<br>
52 * <bNote:</b> If passing real data to fft() forward transform
53 * (especially when using fft() in a loop), explicitly (re-)initialize
54 * the imaginary part to zero (DATA[1][i] = 0.0). Returns 0.
55 *
56 * \param[in] i_sign Direction of transform -1 is normal, +1 is inverse
57 * \param[in,out] data Pointer to complex linear array in row major order
58 * containing data and result
59 * \param[in] NN Value of DATA dimension (dimc * dimr)
60 * \param[in] dimc Value of image column dimension (max power of 2)
61 * \param[in] dimr Value of image row dimension (max power of 2)
62 * \return int always returns 0
63 */
64
65int fft2(int i_sign, double (*data)[2], int NN, int dimc, int dimr)
66{
67#ifdef HAVE_FFTW3_H
68 fftw_plan plan;
69#else
70 fftwnd_plan plan;
71#endif
72 double norm;
73 int i;
74
75 norm = 1.0 / sqrt(NN);
76
77#ifdef HAVE_FFTW3_H
78 plan = fftw_plan_dft_2d(dimr, dimc, data, data,
79 (i_sign < 0) ? FFTW_FORWARD : FFTW_BACKWARD,
80 FFTW_ESTIMATE);
81
82 fftw_execute(plan);
83
84 fftw_destroy_plan(plan);
85#else
86 plan = fftw2d_create_plan(dimc, dimr,
87 (i_sign < 0) ? FFTW_FORWARD : FFTW_BACKWARD,
88 FFTW_ESTIMATE | FFTW_IN_PLACE);
89
90 fftwnd_one(plan, data, data);
91
92 fftwnd_destroy_plan(plan);
93#endif
94
95 for (i = 0; i < NN; i++) {
96 data[i][0] *= norm;
97 data[i][1] *= norm;
98 }
99
100 return 0;
101}
102
103/**
104 * \fn int fft(int i_sign, double *DATA[2], int NN, int dimc, int dimr)
105 *
106 * \brief Fast Fourier Transform for two-dimensional array.
107 *
108 * Fast Fourier Transform for two-dimensional array.<br>
109 * <bNote:</b> If passing real data to fft() forward transform
110 * (especially when using fft() in a loop), explicitly (re-)initialize
111 * the imaginary part to zero (DATA[1][i] = 0.0). Returns 0.
112 *
113 * \param[in] i_sign Direction of transform -1 is normal, +1 is inverse
114 * \param[in,out] DATA Pointer to complex linear array in row major order
115 * containing data and result
116 * \param[in] NN Value of DATA dimension (dimc * dimr)
117 * \param[in] dimc Value of image column dimension (max power of 2)
118 * \param[in] dimr Value of image row dimension (max power of 2)
119 * \return int always returns 0
120 */
121
122int fft(int i_sign, double *DATA[2], int NN, int dimc, int dimr)
123{
124 fftw_complex *data;
125 int i;
126
127 data = (fftw_complex *)G_malloc(NN * sizeof(fftw_complex));
128
129 for (i = 0; i < NN; i++) {
130 c_re(data[i]) = DATA[0][i];
131 c_im(data[i]) = DATA[1][i];
132 }
133
134 fft2(i_sign, data, NN, dimc, dimr);
135
136 for (i = 0; i < NN; i++) {
137 DATA[0][i] = c_re(data[i]);
138 DATA[1][i] = c_im(data[i]);
139 }
140
141 G_free(data);
142
143 return 0;
144}
145
146#endif /* HAVE_FFT */
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150