1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
|
/*
fft.c
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: A program that uses FFTs
*
* Author: E.BERTIN (IAP)
*
* Contents: Routines dealing with double precision FFT.
*
* Last modify: 26/03/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>
#include "define.h"
#include "globals.h"
#include "fft.h"
#include "prefs.h"
#ifdef USE_THREADS
#include "threads.h"
#endif
int firsttimeflag;
#ifdef USE_THREADS
pthread_mutex_t fftmutex;
#endif
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
/****** fft_init ************************************************************
PROTO void fft_init(void)
PURPOSE Initialize the FFT routines
INPUT -.
OUTPUT -.
NOTES Global preferences are used for multhreading.
AUTHOR E. Bertin (IAP)
VERSION 29/11/2006
***/
void fft_init(void)
{
if (!firsttimeflag)
{
#ifdef USE_THREADS
if (!fftw_init_threads())
error(EXIT_FAILURE, "*Error*: thread initialization failed in ", "FFTW");
fftw_plan_with_nthreads(prefs.nthreads);
QPTHREAD_MUTEX_INIT(&fftmutex, NULL);
#endif
firsttimeflag = 1;
}
return;
}
/****** fft_end ************************************************************
PROTO void fft_init(void)
PURPOSE Clear up stuff set by FFT routines
INPUT -.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 29/11/2006
***/
void fft_end(void)
{
if (firsttimeflag)
{
firsttimeflag = 0;
#ifdef USE_THREADS
fftw_cleanup_threads();
QPTHREAD_MUTEX_DESTROY(&fftmutex);
#endif
fftw_cleanup();
}
return;
}
/****** fft_conv ************************************************************
PROTO void fft_conv(double *data1, double *fdata2, int *size)
PURPOSE Optimized 2-dimensional FFT convolution using the FFTW library.
INPUT ptr to the first image,
ptr to the Fourier transform of the second image,
image size vector.
OUTPUT -.
NOTES For data1 and fdata2, memory must be allocated for
size[0]* ... * 2*(size[naxis-1]/2+1) doubles (padding required).
AUTHOR E. Bertin (IAP)
VERSION 26/03/2007
***/
void fft_conv(double *data1, double *fdata2, int *size)
{
fftw_plan plan;
double *fdata1,*fdata1p,*fdata2p,
real,imag, fac;
int i, npix,npix2;
/* Convert axis indexing to that of FFTW */
npix = size[0]*size[1];
npix2 = (((size[0]/2) + 1)*2) * size[1];
/* Forward FFT "in place" for data1 */
#ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif
QFFTWMALLOC(fdata1, double, npix2);
plan = fftw_plan_dft_r2c_2d(size[1], size[0], data1,
(fftw_complex *)fdata1, FFTW_ESTIMATE|FFTW_DESTROY_INPUT);
#ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif
fftw_execute(plan);
#ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif
fftw_destroy_plan(plan);
#ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif
/* Actual convolution (Fourier product) */
fac = 1.0/npix;
fdata1p = fdata1;
fdata2p = fdata2;
for (i=npix2/2; i--; fdata2p+=2)
{
real = *fdata1p **fdata2p - *(fdata1p+1)**(fdata2p+1);
imag = *(fdata1p+1)**fdata2p + *fdata1p**(fdata2p+1);
*(fdata1p++) = fac*real;
*(fdata1p++) = fac*imag;
}
/* Reverse FFT */
#ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif
plan = fftw_plan_dft_c2r_2d(size[1], size[0], (fftw_complex *)fdata1,
data1, FFTW_ESTIMATE|FFTW_DESTROY_INPUT);
#ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif
fftw_execute(plan);
#ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif
fftw_destroy_plan(plan);
/* Free the fdata1 scratch array */
QFFTWFREE(fdata1);
#ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif
return;
}
/****** fft_rtf ************************************************************
PROTO double *fft_rtf(double *data, int *size)
PURPOSE Optimized 2-dimensional FFT "in place" using the FFTW library.
INPUT ptr to the image,
ptr to image size vector.
OUTPUT Pointer to the compressed, memory-allocated Fourier transform.
NOTES Input data may end up corrupted.
AUTHOR E. Bertin (IAP)
VERSION 26/03/2007
***/
double *fft_rtf(double *data, int *size)
{
fftw_plan plan;
fftw_complex *fdata;
int npix2;
/* Convert axis indexing to that of FFTW */
npix2 = ((size[0]/2) + 1) * size[1];
/* Forward FFT "in place" for data1 */
#ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif
QFFTWMALLOC(fdata, fftw_complex, npix2);
plan = fftw_plan_dft_r2c_2d(size[1], size[0], data,
fdata, FFTW_ESTIMATE|FFTW_DESTROY_INPUT);
#ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif
fftw_execute(plan);
#ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif
fftw_destroy_plan(plan);
#ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif
return (double *)fdata;
}
|
