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

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

hht.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
33
34 #include <vector>
35
36 #ifndef __INCLUDE_HHT_H__
37 #define __INCLUDE_HHT_H___
38
39 #ifndef __INCLUDE_MIST_H__
40 #include "mist.h"
41 #endif
42
43
44 #ifndef __INCLUDE_FFT_H__
45 #include "fft/fft.h"
46 #endif
47
48 // mist名前空間の始まり
49 _MIST_BEGIN
50
58
59
60 template < class T >
61 bool emd_envelope(const array< T > &x, const array< T > &y, size_t xsize, array< T > &yout)
62 {
63 T p, sig;
64
65 if( x.size( ) < xsize || y.size( ) < xsize )
66 {
67 return ( false );
68 }
69
70 array< T > u( xsize );
71 array< T > y2( xsize );
72
73 y2[ 0 ] = u[ 0 ] = 0.0;
74
75
76 for( size_t i = 1 ; i < xsize - 1 ; i++ )
77 {
78 sig = ( x[ i ] - x[ i - 1 ] ) / ( x[ i + 1 ] - x[ i - 1 ] );
79 p = sig * y2[ i - 1 ] + 2.0;
80 y2[ i ] = ( sig - 1.0 ) / p;
81 u[ i ] = ( y[ i + 1 ] - y[ i ] ) / ( x[ i + 1 ] - x[ i ] ) - ( y[ i ] - y[ i - 1 ] ) / ( x[ i ] - x[ i - 1 ] );
82 u[ i ] = ( 6.0 * u[ i ] / ( x[ i + 1 ] - x[ i - 1 ] ) - sig * u[ i - 1 ] ) / p;
83 }
84
85 T qn = 0.0, un = 0.0;
86
87 y2[ xsize - 1 ] = ( un - qn * u[ xsize - 2 ] ) / ( qn * y2[ xsize - 2 ] + 1.0 );
88
89 for( size_t i = xsize - 1 ; i > 0 ; i-- )
90 {
91 y2[ i - 1 ] = y2[ i - 1 ] * y2[ i ] + u[ i - 1 ];
92 }
93
94 size_t k = 0;
95 T h, b, a;
96
97 for( size_t i = 0 ; i < yout.size( ) ; i++ )
98 {
99 if( x[ k + 1 ] < static_cast< T >( i ) )
100 {
101 k++;
102 }
103 h = x[ k + 1 ] - x[ k ];
104 a = ( x[ k + 1 ] - static_cast< T >( i ) ) / h;
105 b = ( static_cast< T >( i ) - x[ k ] ) / h;
106 yout[ i ] = a * y[ k ] + b * y[ k + 1 ] + ( ( a * a * a - a ) * y2[ k ] + ( b * b * b - b ) * y2[ k + 1 ] ) * ( h * h ) / 6.0;
107 }
108
109 return ( true );
110 }
111
112
125 template < class T >
126 bool hilbert( const array< T > &in, array< std::complex< T > > &out, bool mirror = false )
127 {
128
129 double f_index = std::log( static_cast< double >( in.size( ) ) ) / std::log( 2.0 );
130 if( mirror )
131 {
132 f_index += 1.0;
133 }
134 array< std::complex< T > > as( static_cast< size_t >( std::pow( 2.0, std::floor( f_index ) + 1.0 ) ) );
135
136 for( size_t i = 0 ; i < in.size( ) ; i ++ )
137 {
138 as[ i ] = in[ i ];
139 }
140
141 for( size_t i = in.size( ) ; i < as.size( ) ; i ++ )
142 {
143 as[ i ] = 0.0;
144 }
145
146 // 端がなめらかになるように折り返す場合
147 if( mirror )
148 {
149 size_t j;
150
151 j = in.size( ) - 3;
152 while( j > 0 && ( in[ in.size( ) - 1 ] - in[ j ] ) * ( in[ in.size( ) - 1 ] - in[ j + 1 ] ) > 0 )
153 {
154 j --;
155 }
156 if( j < in.size( ) / 4 )
157 {
158 j = in.size( ) - 2;
159 }
160
161 for( size_t i = in.size( ) ; i < ( in.size( ) + as.size( ) ) / 2 ; i ++ )
162 {
163 if( j >= 0 )
164 {
165 as[ i ] = std::complex< T >( in[ j -- ], as[ i ].imag( ) );
166 }
167 }
168
169 j = 2;
170 while( j < in.size( ) - 1 && ( in[ 0 ] - in[ j ] ) * ( in [ 0 ] - in[ j - 1 ] ) > 0 )
171 {
172 j ++;
173 }
174 if( j > in.size( ) * 3 / 4 )
175 {
176 j = 1;
177 }
178
179 for( size_t i = as.size( ) - 1 ; i >= ( in.size( ) + as.size( ) ) / 2 ; i -- )
180 {
181 if( j < in.size( ) )
182 {
183 as[ i ] = std::complex< T >( in[ j ++ ], as[ i ].imag( ) );
184 }
185 }
186 }
187
188 _fft( as, as );
189
190 // 正の周波数成分を2倍に
191 for( size_t i = 1 ; i < as.size( ) / 2 ; i ++ )
192 {
193 as[ i ] *= 2.0;
194 }
195
196 // 負の周波数成分を0に
197 for( size_t i = as.size( ) / 2 ; i < as.size( ) ; i ++ )
198 {
199 as( i ) = 0.0;
200 }
201
202 _ifft( as, as );
203
204 out.resize( in.size( ) );
205
206 for( size_t i = 0 ; i < in.size( ) ; i++ )
207 {
208 out[ i ] = std::complex< T >( in[ i ], as[ i ].imag( ) );
209 }
210
211 return ( true );
212 }
213
233 template < class T1, class T2 >
234 size_t emd( const array< T1 > &in, array< T2 > &imf, double sd = 0.3 )
235 {
236 std::vector< T2 > imf_v;
237
238 // 初期化
239 T2 r( in.size( ) );
240 std::copy( in.begin( ), in.end( ), r.begin( ) );
241 size_t minmax_old = in.size( );
242
243 while( true )
244 {
245 T2 h, h_new;
246 size_t minmax;
247
248 h = h_new = r;
249 while( true )
250 {
251 size_t max_i = 0, min_i = 0;
252 minmax = 0;
253 array< T1 > max_x( in.size( ) );
254 array< T1 > max_y( in.size( ) );
255 array< T1 > min_x( in.size( ) );
256 array< T1 > min_y( in.size( ) );
257 max_x[ 0 ] = min_x[ 0 ] = 0.0;
258 max_x[ 1 ] = min_x[ 1 ] = static_cast< T1 >( in.size( ) - 1 );
259
260 // 極値を抽出
261 for( size_t i = 1 ; i < in.size( ) - 1 ; i ++ )
262 {
263 if( h[ i ] > h[ i - 1 ] && h[ i ] > h[ i + 1 ] )
264 {
265 max_i ++;
266 max_x[ max_i ] = static_cast< T1 >( i );
267 max_y[ max_i ] = h[ i ];
268 minmax ++;
269 }
270 if( h[ i ] < h[ i - 1 ] && h[ i ] < h[ i + 1 ] )
271 {
272 min_i ++;
273 min_x[ min_i ] = static_cast< T1 >( i );
274 min_y[ min_i ] = h[ i ];
275 minmax ++;
276 }
277 }
278
279 if( minmax < 1 )
280 {
281 break;
282 }
283
284 // 左端を調節
285 if( max_x[ 1 ] < min_x[ 1 ] )
286 {
287 if( min_i > 0 && min_y[ 1 ] < h[ 0 ] )
288 {
289 min_x[ 0 ] = max_x[ 1 ] * 2.0 - min_x[ 1 ];
290 min_y[ 0 ] = min_y[ 1 ];
291 }
292 else
293 {
294 min_x[ 0 ] = 0.0;
295 min_y[ 0 ] = h[ 0 ];
296 }
297 max_x[ 0 ] = min_x[ 0 ] * 2.0 - max_x[ 1 ];
298 max_y[ 0 ] = max_y[ 1 ];
299 }
300 else
301 {
302 if( max_i > 0 && max_y[ 1 ] > h[ 0 ] )
303 {
304 max_x[ 0 ] = min_x[ 1 ] * 2.0 - max_x[ 1 ];
305 max_y[ 0 ] = max_y[ 1 ];
306 }
307 else
308 {
309 max_x[ 0 ] = 0.0;
310 max_y[ 0 ] = h[ 0 ];
311 }
312 min_x[ 0 ] = max_x[ 0 ] * 2.0 - min_x[ 1 ];
313 min_y[ 0 ] = min_y[ 1 ];
314 }
315 if( max_x[ 0 ] > 0 )
316 {
317 max_x[ 0 ] = 0.0;
318 }
319 if( min_x[ 0 ] > 0 )
320 {
321 min_x[ 0 ] = 0.0;
322 }
323
324 // 右端を調節
325 if( max_x[ max_i ] > min_x[ min_i ] )
326 {
327 if( min_i > 0 && min_y[ min_i ] < h[ in.size( ) - 1 ] )
328 {
329 min_x[ min_i + 1 ] = max_x[ max_i ] * 2.0 - min_x[ min_i ];
330 min_y[ min_i + 1 ] = min_y[ min_i ];
331 }
332 else
333 {
334 min_x[ min_i + 1 ] = static_cast< T1 >( in.size( ) - 1 );
335 min_y[ min_i + 1 ] = h[ in.size( ) - 1 ];
336 }
337 max_x[ max_i + 1 ] = min_x[ min_i + 1 ] * 2.0 - max_x[ max_i ];
338 max_y[ max_i + 1 ] = max_y[ max_i ];
339 }
340 else
341 {
342 if( max_i > 0 && max_y[ max_i ] > h[ in.size( ) - 1 ] )
343 {
344 max_x[ max_i + 1 ] = min_x[ min_i ] * 2.0 - max_x[ max_i ];
345 max_y[ max_i + 1 ] = max_y[ max_i ];
346 }
347 else
348 {
349 max_x[ max_i + 1 ] = static_cast< T1 >( in.size( ) - 1 );
350 max_y[ max_i + 1 ] = h[ in.size( ) - 1 ];
351 }
352 min_x[ min_i + 1 ] = max_x[ max_i + 1 ] * 2.0 - min_x[ min_i ];
353 min_y[ min_i + 1 ] = min_y[ min_i ];
354 }
355 if( max_x[ max_i + 1 ] < static_cast< T1 >( in.size( ) - 1 ) )
356 {
357 max_x[ max_i + 1 ] = static_cast< T1 >( in.size( ) - 1 );
358 }
359 if( min_x[ min_i + 1 ] < static_cast< T1 >( in.size( ) - 1 ) )
360 {
361 min_x[ min_i + 1 ] = static_cast< T1 >( in.size( ) - 1 );
362 }
363 array< T1 > max_e( in.size( ) );
364 array< T1 > min_e( in.size( ) );
365 array< T1 > mean_e( in.size( ) );
366
367 emd_envelope( max_x, max_y, max_i + 2, max_e );
368 emd_envelope( min_x, min_y, min_i + 2, min_e );
369
370 T1 sd1 = 0.0;
371
372 // 包短線の中間の線を抽出し，シフト処理
373 for( size_t i = 0 ; i < in.size( ) ; i ++ )
374 {
375 mean_e[ i ] = ( max_e[ i ] + min_e[ i ] ) / 2.0;
376 h_new[ i ] = h[ i ] - mean_e[ i ];
377 sd1 += ( h[ i ] - h_new[ i ] ) * ( h[ i ] - h_new[ i ] ) / ( h[ i ] * h[ i ] + 0.00001 );
378 }
379
380 // IMFかどうかの判定
381 if( static_cast< double >( sd1 ) < sd )
382 {
383 break;
384 }
385
386 h = h_new;
387 }
388
389 if( minmax < 1 || minmax > minmax_old )
390 {
391 break;
392 }
393 minmax_old = minmax;
394
395 for( size_t i = 0 ; i < in.size( ) ; i ++ )
396 {
397 r[ i ] -= h[ i ];
398 }
399 imf_v.push_back( h );
400
401 }
402
403 imf.resize( imf_v.size( ) + 1 );
404 for( size_t i = 0 ; i < imf_v.size( ) ; i ++ )
405 {
406 imf[ i ] = imf_v[ i ];
407 }
408 imf( imf_v.size( ) ) = r;
409
410 return ( imf_v.size( ) );
411 }
412
414 // グループの終わり
415
416 // mist名前空間の終わり
417 _MIST_END
418
419
420 #endif // __INCLUDE_HHT_H__
421