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 
Generated on Wed Nov 12 2014 19:44:15 for MIST by doxygen 1.8.1.2