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

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

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