- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

fft.h
説明を見る。
1 //
2 // Copyright (c) 2003-2011, MIST Project, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28
37
38 #ifndef __INCLUDE_FFT_H__
39 #define __INCLUDE_FFT_H__
40
41
42 #ifndef __INCLUDE_MIST_H__
43 #include "../mist.h"
44 #endif
45
46 #ifndef __INCLUDE_MIST_THREAD__
47 #include "../thread.h"
48 #endif
49
50 #ifndef __INCLUDE_FFT_UTIL_H__
51 #include "fft_util.h"
52 #endif
53
54
55 #if defined( __MIST_WINDOWS__ ) && __MIST_WINDOWS__ > 0
56 #define USE_CDFT_WINTHREADS
57 #define USE_FFT2D_WINTHREADS
58 #define USE_FFT3D_WINTHREADS
59 #else
60 #define USE_CDFT_PTHREADS
61 #define USE_FFT2D_PTHREADS
62 #define USE_FFT3D_PTHREADS
63 #endif
64
65 #include "fftsg.h"
66 #include "fftsg2d.h"
67 #include "fftsg3d.h"
68
69 // mist名前空間の始まり
70 _MIST_BEGIN
71
72
75
83
84
85
86 template < class T1, class T2, class Allocator1, class Allocator2 >
87 bool _fft( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out )
88 {
89 if( !__fft_util__::size_check( ( unsigned int ) in.size( ) ) )
90 {
91 return( false );
92 }
93
94 // 1次元高速フーリエ変換
95 typedef typename Allocator1::size_type size_type;
96 size_type i;
97 __fft_util__::FFT_MEMORY1 mem;
98
99 if( !__fft_util__::allocate_memory( mem,
100 in.size( ) * 2,
101 static_cast< size_t >( std::sqrt( static_cast< double >( in.size( ) ) ) + 3 ),
102 in.size( ) / 2 ) )
103 {
104 __fft_util__::deallocate_memory( mem );
105 return( false );
106 }
107
108 double *data = mem.data;
109 double *w = mem.w;
110 int *ip = mem.ip;
111
112 for( i = 0 ; i < in.size( ) ; i++ )
113 {
114 std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in[ i ] ) );
115 data[ i * 2 ] = c.real( );
116 data[ i * 2 + 1 ] = c.imag( );
117 }
118
119 ip[ 0 ] = 0;
120
121 ooura_fft::cdft( static_cast< int >( in.size( ) * 2 ), -1, data, ip, w );
122
123 out.resize( in.size( ) );
124
125 for( i = 0 ; i < out.size( ) ; i++ )
126 {
127 out[ i ] = __fft_util__::convert_complex< T2 >::convert_from( data[ 2 * i ], data[ 2 * i + 1 ] );
128 }
129
130 __fft_util__::deallocate_memory( mem );
131
132 return( true );
133 }
134
135
136
137 template < class T1, class T2, class Allocator1, class Allocator2 >
138 bool _ifft( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out )
139 {
140 if( !__fft_util__::size_check( ( unsigned int ) in.size( ) ) )
141 {
142 return( false );
143 }
144
145 // 1次元高速フーリエ変換
146 typedef typename Allocator1::size_type size_type;
147 size_type i;
148 __fft_util__::FFT_MEMORY1 mem;
149
150 if( !__fft_util__::allocate_memory( mem,
151 in.size( ) * 2,
152 static_cast< size_t >( std::sqrt( static_cast< double >( in.size( ) ) ) + 3 ),
153 in.size( ) / 2 ) )
154 {
155 __fft_util__::deallocate_memory( mem );
156 return( false );
157 }
158
159 double *data = mem.data;
160 double *w = mem.w;
161 int *ip = mem.ip;
162
163 for( i = 0 ; i < in.size( ) ; i++ )
164 {
165 std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in[ i ] ) );
166 data[ i * 2 ] = c.real( );
167 data[ i * 2 + 1 ] = c.imag( );
168 }
169
170 ip[ 0 ] = 0;
171
172 ooura_fft::cdft( static_cast< int >(in.size( ) * 2), 1, data, ip, w );
173
174 out.resize( in.size( ) );
175
176 double __value__ = 1.0 / in.size( );
177 for( i = 0 ; i < out.size( ) ; i++ )
178 {
179 out[ i ] = __fft_util__::convert_complex< T2 >::convert_from( data[ 2 * i ] * __value__, data[ 2 * i + 1 ] * __value__ );
180 }
181
182 __fft_util__::deallocate_memory( mem );
183
184 return( true );
185 }
186
187
199 template < class T1, class T2, class Allocator1, class Allocator2 >
200 bool fft( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out )
201 {
202 return _fft( in, out );
203 }
204
205
206
217 template < class T1, class T2, class Allocator1, class Allocator2 >
218 bool ifft( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out )
219 {
220 return _ifft( in, out );
221 }
222
223
224
236 template < class T1, class T2, class Allocator1, class Allocator2 >
237 bool fft( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out )
238 {
239 if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) || !__fft_util__::size_check( ( unsigned int ) in.height( ) ) )
240 {
241 return( false );
242 }
243
244 typedef typename Allocator1::size_type size_type;
245 size_type i, j, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
246 __fft_util__::FFT_MEMORY2 mem;
247
248 if( !__fft_util__::allocate_memory( mem,
249 in.width( ),
250 in.height( ) * 2,
251 8 * in.width( ) * FFT2D_MAX_THREADS,
252 static_cast< size_t >( std::sqrt( static_cast< double >( size ) ) + 3 ),
253 size / 2 ) )
254 {
255 __fft_util__::deallocate_memory( mem );
256 return( false );
257 }
258
259 double **data = mem.data;
260 double *t = mem.t;
261 double *w = mem.w;
262 int *ip = mem.ip;
263
264 for( i = 0 ; i < in.width( ) ; i++ )
265 {
266 for( j = 0 ; j < in.height( ) ; j++ )
267 {
268 std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j ) ) );
269 data[ i ][ 2 * j ] = c.real( );
270 data[ i ][ 2 * j + 1 ] = c.imag( );
271 }
272 }
273
274 ip[ 0 ] = 0;
275
276 ooura_fft::cdft2d( static_cast< int >( in.width( ) ), static_cast< int >( in.height( ) * 2 ), -1, data, t, ip, w );
277
278
279 out.resize( in.width( ), in.height( ) );
280
281 for( i = 0 ; i < out.width( ) ; i++ )
282 {
283 for( j = 0 ; j < out.height( ) ; j++ )
284 {
285 out( i, j ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ 2 * j ], data[ i ][ 2 * j + 1 ] );
286 }
287 }
288
289 __fft_util__::deallocate_memory( mem );
290
291 return( true );
292 }
293
294
295
306 template < class T1, class T2, class Allocator1, class Allocator2 >
307 bool ifft( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out )
308 {
309 if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) || !__fft_util__::size_check( ( unsigned int ) in.height( ) ) )
310 {
311 return( false );
312 }
313
314 typedef typename Allocator1::size_type size_type;
315 size_type i, j, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
316 __fft_util__::FFT_MEMORY2 mem;
317
318 if( !__fft_util__::allocate_memory( mem,
319 in.width( ),
320 in.height( ) * 2,
321 8 * in.width( ) * FFT2D_MAX_THREADS,
322 static_cast< size_t >( std::sqrt( static_cast< double >( size ) ) + 3 ),
323 size / 2 ) )
324 {
325 __fft_util__::deallocate_memory( mem );
326 return( false );
327 }
328
329 double **data = mem.data;
330 double *t = mem.t;
331 double *w = mem.w;
332 int *ip = mem.ip;
333
334 for( i = 0 ; i < in.width( ) ; i++ )
335 {
336 for( j = 0 ; j < in.height( ) ; j++ )
337 {
338 std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j ) ) );
339 data[ i ][ 2 * j ] = c.real( );
340 data[ i ][ 2 * j + 1 ] = c.imag( );
341 }
342 }
343
344 ip[ 0 ] = 0;
345
346 ooura_fft::cdft2d( static_cast< int >(in.width( )), static_cast< int >(in.height( ) * 2), 1, data, t, ip, w );
347
348 out.resize( in.width( ), in.height( ) );
349
350 double __value__ = 1.0 / in.size( );
351 for( i = 0 ; i < out.width( ) ; i++ )
352 {
353 for( j = 0 ; j < out.height( ) ; j++ )
354 {
355 out( i, j ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ 2 * j ] * __value__, data[ i ][ 2 * j + 1 ] * __value__ );
356 }
357 }
358
359 __fft_util__::deallocate_memory( mem );
360
361 return( true );
362 }
363
364
365
377 template < class T1, class T2, class Allocator1, class Allocator2 >
378 bool fft( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
379 {
380 if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) ||
381 !__fft_util__::size_check( ( unsigned int ) in.height( ) ) ||
382 !__fft_util__::size_check( ( unsigned int ) in.depth( ) ) )
383 {
384 return( false );
385 }
386
387 typedef typename Allocator1::size_type size_type;
388 size_type i, j, k, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
389 __fft_util__::FFT_MEMORY3 mem;
390
391 if( !__fft_util__::allocate_memory( mem,
392 in.width( ),
393 in.height( ),
394 in.depth( ) * 2,
395 8 * size * FFT3D_MAX_THREADS,
396 static_cast< size_t >( std::sqrt( static_cast< double >( size ) ) + 3 ),
397 size / 2 ) )
398 {
399 __fft_util__::deallocate_memory( mem );
400 return( false );
401 }
402
403 double ***data = mem.data;
404 double *t = mem.t;
405 double *w = mem.w;
406 int *ip = mem.ip;
407
408 for( i = 0 ; i < in.width( ) ; i++ )
409 {
410 for( j = 0 ; j < in.height( ) ; j++ )
411 {
412 for( k = 0 ; k < in.depth( ) ; k++ )
413 {
414 std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j, k ) ) );
415 data[ i ][ j ][ 2 * k ] = c.real( );
416 data[ i ][ j ][ 2 * k + 1 ] = c.imag( );
417 }
418 }
419 }
420
421 ip[ 0 ] = 0;
422
423 ooura_fft::cdft3d( static_cast< int >( in.width( ) ), static_cast< int >( in.height( ) ), static_cast< int >( in.depth( ) * 2 ), -1, data, t, ip, w );
424
425 out.resize( in.width( ), in.height( ), in.depth( ) );
426
427 for( i = 0 ; i < out.width( ) ; i++ )
428 {
429 for( j = 0 ; j < out.height( ) ; j++ )
430 {
431 for( k = 0 ; k < out.depth( ) ; k++ )
432 {
433 out( i, j, k ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ][ 2 * k ], data[ i ][ j ][ 2 * k + 1 ] );
434 }
435 }
436 }
437
438 __fft_util__::deallocate_memory( mem );
439
440 return( true );
441 }
442
443
444
445
456 template < class T1, class T2, class Allocator1, class Allocator2 >
457 bool ifft( const array3< T1, Allocator1 > &in, array3< T2 , Allocator2 > &out )
458 {
459 if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) ||
460 !__fft_util__::size_check( ( unsigned int ) in.height( ) ) ||
461 !__fft_util__::size_check( ( unsigned int ) in.depth( ) ) )
462 {
463 return( false );
464 }
465
466 typedef typename Allocator1::size_type size_type;
467 size_type i, j, k, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
468 __fft_util__::FFT_MEMORY3 mem;
469
470 if( !__fft_util__::allocate_memory( mem,
471 in.width( ),
472 in.height( ),
473 in.depth( ) * 2,
474 8 * size * FFT3D_MAX_THREADS,
475 static_cast< size_t >( std::sqrt( static_cast< double >( size ) ) + 3 ),
476 size / 2 ) )
477 {
478 __fft_util__::deallocate_memory( mem );
479 return( false );
480 }
481
482 double ***data = mem.data;
483 double *t = mem.t;
484 double *w = mem.w;
485 int *ip = mem.ip;
486
487 for( i = 0 ; i < in.width( ) ; i++ )
488 {
489 for( j = 0 ; j < in.height( ) ; j++ )
490 {
491 for( k = 0 ; k < in.depth( ) ; k++ )
492 {
493 std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j, k ) ) );
494 data[ i ][ j ][ 2 * k ] = c.real( );
495 data[ i ][ j ][ 2 * k + 1 ] = c.imag( );
496 }
497 }
498 }
499
500 ip[ 0 ] = 0;
501
502 ooura_fft::cdft3d( static_cast< int >( in.width( ) ), static_cast< int >( in.height( ) ), static_cast< int >( in.depth( ) * 2 ), 1, data, t, ip, w );
503
504 out.resize( in.width( ), in.height( ), in.depth( ) );
505
506 double __value__ = 1.0 / in.size( );
507 for( i = 0 ; i < out.width( ) ; i++ )
508 {
509 for( j = 0 ; j < out.height( ) ; j++ )
510 {
511 for( k = 0 ; k < out.depth( ) ; k++ )
512 {
513 out( i, j, k ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ][ 2 * k ] * __value__, data[ i ][ j ][ 2 * k + 1 ] * __value__ );
514 }
515 }
516 }
517
518 __fft_util__::deallocate_memory( mem );
519
520 return( true );
521 }
522
523
525 // FFT グループの終わり
526
528 // Fourier グループの終わり
529
530
531 _MIST_END
532
533 #endif
534
535