dst.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_DST_H__
39 #define __INCLUDE_DST_H__
40 
41 
42 #ifndef __INCLUDE_MIST_H__
43 #include "../mist.h"
44 #endif
45 
46 #ifndef __INCLUDE_FFT_UTIL_H__
47 #include "fft_util.h"
48 #endif
49 
50 #if defined( __MIST_WINDOWS__ ) && __MIST_WINDOWS__ > 0
51 #define USE_CDFT_WINTHREADS
52 #define USE_FFT2D_WINTHREADS
53 #define USE_FFT3D_WINTHREADS
54 #else
55 #define USE_CDFT_PTHREADS
56 #define USE_FFT2D_PTHREADS
57 #define USE_FFT3D_PTHREADS
58 #endif
59 
60 #include "fftsg.h"
61 #include "fftsg2d.h"
62 #include "fftsg3d.h"
63 
64 // mist名前空間の始まり
65 _MIST_BEGIN
66 
67 
68 
71 
79 
80 
92 template < class T1, class T2, class Allocator1, class Allocator2 >
93 bool dst( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out )
94 {
95  if( !__fft_util__::size_check( ( unsigned int ) in.size( ) ) )
96  {
97  return( false );
98  }
99 
100  typedef typename Allocator1::size_type size_type;
101  size_type i;
102  __fft_util__::FFT_MEMORY1 mem;
103 
104  if( !__fft_util__::allocate_memory( mem,
105  in.size( ),
106  static_cast< size_t >( std::sqrt( static_cast< double >( in.size( ) / 2 ) ) + 3 ),
107  in.size( ) * 5 / 4 ) )
108  {
109  __fft_util__::deallocate_memory( mem );
110  return( false );
111  }
112 
113  double *data = mem.data;
114  double *w = mem.w;
115  int *ip = mem.ip;
116 
117  for( i = 0 ; i < in.size( ) ; i++ )
118  {
119  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in[ i ] ) );
120  data[ i ] = c.real( );
121  }
122 
123  ip[ 0 ] = 0;
124 
125  ooura_fft::ddst( static_cast< int >( in.size( ) ), -1, data, ip, w );
126 
127 
128  out.resize( in.size( ) );
129 
130  out[ 0 ] = __fft_util__::convert_complex< T2 >::convert_from( data[ in.size( ) - 1 ], 0.0 );
131 
132  for( i = 1 ; i < out.size( ) ; i++ )
133  {
134  out[ i ] = __fft_util__::convert_complex< T2 >::convert_from( data[ i ], 0.0 );
135  }
136 
137  __fft_util__::deallocate_memory( mem );
138 
139  return( true );
140 }
141 
142 
143 
154 template < class T1, class T2, class Allocator1, class Allocator2 >
155 bool idst( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out )
156 {
157  if( !__fft_util__::size_check( ( unsigned int ) in.size( ) ) )
158  {
159  return( false );
160  }
161 
162  typedef typename Allocator1::size_type size_type;
163  size_type i;
164  __fft_util__::FFT_MEMORY1 mem;
165 
166  if( !__fft_util__::allocate_memory( mem,
167  in.size( ),
168  static_cast< size_t >( std::sqrt( static_cast< double >( in.size( ) / 2 ) ) + 3 ),
169  in.size( ) * 5 / 4 ) )
170  {
171  __fft_util__::deallocate_memory( mem );
172  return( false );
173  }
174 
175  double *data = mem.data;
176  double *w = mem.w;
177  int *ip = mem.ip;
178 
179  for( i = 1 ; i < in.size( ) ; i++ )
180  {
181  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in[ i ] ) );
182  data[ i ] = c.real( );
183  }
184 
185  data[ 0 ] = data[ in.size( ) - 1 ];
186 
187  ip[ 0 ] = 0;
188  data[ 0 ] *= 0.5;
189 
190  ooura_fft::ddst( static_cast< int >( in.size( ) ), 1, data, ip, w );
191 
192 
193  out.resize( in.size( ) );
194 
195  for( i = 0 ; i < out.size( ) ; i++ )
196  {
197  out[ i ] = __fft_util__::convert_complex< T2 >::convert_from( data[ i ] * 2.0 / out.size( ), 0.0 );
198  }
199 
200  __fft_util__::deallocate_memory( mem );
201 
202  return( true );
203 }
204 
205 
206 
218 template < class T1, class T2, class Allocator1, class Allocator2 >
219 bool dst( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out )
220 {
221  if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) || !__fft_util__::size_check( ( unsigned int ) in.height( ) ) )
222  {
223  return( false );
224  }
225 
226  typedef typename Allocator1::size_type size_type;
227  size_type i, j, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
228  __fft_util__::FFT_MEMORY2 mem;
229 
230  if( !__fft_util__::allocate_memory( mem,
231  in.width( ),
232  in.height( ),
233  4 * in.width( ) * FFT2D_MAX_THREADS,
234  static_cast< size_t >( std::sqrt( static_cast< double >( size / 2 ) ) + 3 ),
235  size * 3 / 2 ) )
236  {
237  __fft_util__::deallocate_memory( mem );
238  return( false );
239  }
240 
241  double **data = mem.data;
242  double *t = mem.t;
243  double *w = mem.w;
244  int *ip = mem.ip;
245 
246  for( i = 0 ; i < in.width( ) ; i++ )
247  {
248  for( j = 0 ; j < in.height( ) ; j++ )
249  {
250  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j ) ) );
251  data[ i ][ j ] = c.real( );
252  }
253  }
254 
255  ip[ 0 ] = 0;
256 
257  ooura_fft::ddst2d( static_cast< int >( in.width( ) ), static_cast< int >( in.height( ) ), -1, data, t, ip, w );
258 
259 
260  out.resize( in.width( ), in.height( ) );
261 
262  for( i = 0 ; i < out.width( ) ; i++ )
263  {
264  for( j = 0 ; j < out.height( ) ; j++ )
265  {
266  out( i, j ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ], 0.0 );
267  }
268  }
269 
270  __fft_util__::deallocate_memory( mem );
271 
272  return( true );
273 }
274 
275 
276 
287 template < class T1, class T2, class Allocator1, class Allocator2 >
288 bool idst( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out )
289 {
290  if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) || !__fft_util__::size_check( ( unsigned int ) in.height( ) ) )
291  {
292  return( false );
293  }
294 
295  typedef typename Allocator1::size_type size_type;
296  size_type i, j, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
297  __fft_util__::FFT_MEMORY2 mem;
298 
299  if( !__fft_util__::allocate_memory( mem,
300  in.width( ),
301  in.height( ),
302  4 * in.width( ) * FFT2D_MAX_THREADS,
303  static_cast< size_t >( std::sqrt( static_cast< double >( size / 2 ) ) + 3 ),
304  size * 3 / 2 ) )
305  {
306  __fft_util__::deallocate_memory( mem );
307  return( false );
308  }
309 
310  double **data = mem.data;
311  double *t = mem.t;
312  double *w = mem.w;
313  int *ip = mem.ip;
314 
315  for( i = 0 ; i < in.width( ) ; i++ )
316  {
317  for( j = 0 ; j < in.height( ) ; j++ )
318  {
319  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j ) ) );
320  data[ i ][ j ] = c.real( );
321  }
322  }
323 
324  ip[ 0 ] = 0;
325 
326  for( i = 0 ; i < in.width( ) ; i++ )
327  {
328  data[ i ][ 0 ] *= 0.5;
329  }
330 
331  for( j = 0 ; j < in.height( ) ; j++ )
332  {
333  data[ 0 ][ j ] *= 0.5;
334  }
335 
336  ooura_fft::ddst2d( static_cast< int >( in.width( ) ), static_cast< int >( in.height( ) ), 1, data, t, ip, w );
337 
338 
339  out.resize( in.width( ), in.height( ) );
340 
341  for( i = 0 ; i < out.width( ) ; i++ )
342  {
343  for( j = 0 ; j < out.height( ) ; j++ )
344  {
345  out( i, j ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ] * 4.0 / in.size( ), 0.0 );
346  }
347  }
348 
349  __fft_util__::deallocate_memory( mem );
350 
351  return( true );
352 }
353 
354 
355 
367 template < class T1, class T2, class Allocator1, class Allocator2 >
368 bool dst( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
369 {
370  if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) ||
371  !__fft_util__::size_check( ( unsigned int ) in.height( ) ) ||
372  !__fft_util__::size_check( ( unsigned int ) in.depth( ) ) )
373  {
374  return( false );
375  }
376 
377  typedef typename Allocator1::size_type size_type;
378  size_type i, j, k, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
379  __fft_util__::FFT_MEMORY3 mem;
380 
381  if( !__fft_util__::allocate_memory( mem,
382  in.width( ),
383  in.height( ),
384  in.depth( ),
385  4 * size * FFT3D_MAX_THREADS,
386  static_cast< size_t >( std::sqrt( static_cast< double >( size / 2 ) ) + 3 ),
387  size * 3 / 2 ) )
388  {
389  __fft_util__::deallocate_memory( mem );
390  return( false );
391  }
392 
393  double ***data = mem.data;
394  double *t = mem.t;
395  double *w = mem.w;
396  int *ip = mem.ip;
397 
398  for( i = 0 ; i < in.width( ) ; i++ )
399  {
400  for( j = 0 ; j < in.height( ) ; j++ )
401  {
402  for( k = 0 ; k < in.depth( ) ; k++ )
403  {
404  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j, k ) ) );
405  data[ i ][ j ][ k ] = c.real( );
406  }
407  }
408  }
409 
410  ip[ 0 ] = 0;
411 
412  ooura_fft::ddst3d( static_cast< int >( in.width( ) ), static_cast< int >( in.height( ) ), static_cast< int >( in.depth( ) ), -1, data, t, ip, w );
413 
414 
415  out.resize( in.width( ), in.height( ), in.depth( ) );
416 
417  for( i = 0 ; i < out.width( ) ; i++ )
418  {
419  for( j = 0 ; j < out.height( ) ; j++ )
420  {
421  for( k = 0 ; k < out.depth( ) ; k++ )
422  {
423  out( i, j, k ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ][ k ], 0.0 );
424  }
425  }
426  }
427 
428  __fft_util__::deallocate_memory( mem );
429 
430  return( true );
431 }
432 
433 
434 
445 template < class T1, class T2, class Allocator1, class Allocator2 >
446 bool idst( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
447 {
448  if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) ||
449  !__fft_util__::size_check( ( unsigned int ) in.height( ) ) ||
450  !__fft_util__::size_check( ( unsigned int ) in.depth( ) ) )
451  {
452  return( false );
453  }
454 
455  typedef typename Allocator1::size_type size_type;
456  size_type i, j, k, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
457  __fft_util__::FFT_MEMORY3 mem;
458 
459  if( !__fft_util__::allocate_memory( mem,
460  in.width( ),
461  in.height( ),
462  in.depth( ),
463  4 * size * FFT3D_MAX_THREADS,
464  static_cast< size_t >( std::sqrt( static_cast< double >( size / 2 ) ) + 3 ),
465  size * 3 / 2 ) )
466  {
467  __fft_util__::deallocate_memory( mem );
468  return( false );
469  }
470 
471  double ***data = mem.data;
472  double *t = mem.t;
473  double *w = mem.w;
474  int *ip = mem.ip;
475 
476  for( i = 0 ; i < in.width( ) ; i++ )
477  {
478  for( j = 0 ; j < in.height( ) ; j++ )
479  {
480  for( k = 0 ; k < in.depth( ) ; k++ )
481  {
482  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j, k ) ) );
483  data[ i ][ j ][ k ] = c.real( );
484  }
485  }
486  }
487 
488  ip[ 0 ] = 0;
489 
490 
491  for( i = 0 ; i < in.width( ) ; i++ )
492  {
493  for( j = 0 ; j < in.height( ) ; j++ )
494  {
495  data[ i ][ j ][ 0 ] *= 0.5;
496  }
497 
498  for( k = 0 ; k < in.depth( ) ; k++ )
499  {
500  data[ i ][ 0 ][ k ] *= 0.5;
501  }
502  }
503 
504  for( j = 0 ; j < in.height( ) ; j++ )
505  {
506  for( k = 0 ; k < in.depth( ) ; k++ )
507  {
508  data[ 0 ][ j ][ k ] *= 0.5;
509  }
510  }
511 
512  ooura_fft::ddst3d( static_cast< int >( in.width( ) ), static_cast< int >( in.height( ) ), static_cast< int >( in.depth( ) ), 1, data, t, ip, w );
513 
514  out.resize( in.width( ), in.height( ), in.depth( ) );
515 
516  for( i = 0 ; i < out.width( ) ; i++ )
517  {
518  for( j = 0 ; j < out.height( ) ; j++ )
519  {
520  for( k = 0 ; k < out.depth( ) ; k++ )
521  {
522  out( i, j, k ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ][ k ] * 8.0 / in.size( ), 0.0 );
523  }
524  }
525  }
526 
527  __fft_util__::deallocate_memory( mem );
528 
529  return( true );
530 }
531 
532 
534 // DST グループの終わり
535 
537 // Fourier グループの終わり
538 
539 
540 
541 
542 _MIST_END
543 
544 #endif
545 
546 
Generated on Wed Nov 12 2014 19:44:14 for MIST by doxygen 1.8.1.2