;+ The webpage for this routine is, ; http://www.davidpace.com/physics/graduate-school/farge-wavelet-method.htm ; NAME: ; ; fargeMethod ; ; ; PURPOSE: ; ; Applies the method of M. Farge, K. Schneider, and P. Devynck, ; Phys. Plasmas, 13, p. 042304 (2006), to separate time series ; signals into coherent and noise components. This file contains ; additional routines that may be used for testing and learning ; purposes. ; ; ; EXAMPLE: ; ; Generate test signals with a low noise amplitude, ;IDL> q = sigs( 0.2 ) ; ; Plot the coherent and noise parts of this test signal, ;IDL> plot, q.sig ;IDL> oplot, q.noi ; ; Separate the test signal using the wavelet method, ;IDL> w = wvint( q.sig + q.noi ) ; ; Plot the coherent and noise parts as determined by the routine, ;IDL> plot, w.coh ;IDL> oplot, w.noi ; ; ; MODIFICATION HISTORY: ; ; May 11, 2007 - Written by David Pace ; May 14, 2007 - Removed extra FOR loops. Added wavelet frequency ; calculation (to be used for flatness results). ; ;-------------------------------------------- ;******************************************** ;SIGS ; Generate test signals for testing the routine and learning its ; calling sequence. ; ;Inputs - ; snr - Signal to Noise Ratio: type FLOAT, this value multiplies the ; normalized noise array. The noise array will always be within ; the range of plus/minus SNR. ; ;Outputs - ; STRUCTURE: ; sig - The coherent signal. ; noi - The noise signal. ; ;Notes - ; The testing signal for the wavelet coherent component extraction ; technique is generated by adding the results of this function. For ; example, the input to the extraction function, F_test, may be ; written, ; F_test = sig + noi ; function sigs, snr ;create noise see = 100L noi = randomu( see, 16384, /normal ) noi = noi / max( noi ) noi = snr * noi ;create coherent signal w1 = findgen(2048) / 2e3 w1 = 5. * w1 + 0.2 ;linear ramp increasing w11 = reverse( w1 ) ;linear ramp decreasing w2 = findgen(4096) w2 = ( w2 / 2e3 )^2. - exp( -1. * w2 ) + 1. ;mostly exponential increasing w3 = findgen( 4096 ) / 2e3 w3 = sin( 2. * !PI * 2. * w3 ) + [ w11, w11 ]/5. ;sine on decreasing sawtooth w4 = replicate( 3., 4096 ) ;constant value, step function tmp1 = [ w1, w11, w2, w3, w4 ] ;concatenate into one array tmp1 = tmp1 / max( tmp1 ) ;normalize to maximum value of unity return, { sig : tmp1, \$ ;coherent signal noi : noi } ;noise signal end ;******************************************** ;WVINT ; Applies the Farge method to separate an input signal into coherent ; and noise components. ; ;Inputs - ; data - A one-dimensional array of numbers, FLOAT or DOUBLE type. ; ; daub - KEYWORD, if set, then the Daubechies wavelet basis will be ; used. The default is to use the Coiflet wavelet. ; ;Outputs - ; STRUCTURE: ; raw - Raw Signal for Analysis, contains only that part of the ; input array that is actually analyzed. For instance, if the ; input does not have the correct number of data points (any ; power of two), then RAW will contain only that portion of ; the input data that was actually processed. ; ; coh - Coherent Component, that part of the input signal that is ; considered coherent by this method. ; ; var - Variance, of the wavelet coefficients according to ; iteration. ; ; thresh - Threshold, for determining whether a wavelet coefficient ; is noise. This array provides the threshold for each ; iteration of the routine. ; ; inc - Incoherent Component, that part of the input signal that ; remains after subtracting out the coherent part. ; ; nnoise - Number of wavelet coefficients deemed noise at the final ; iteration. ; ; nold - Same as NNOISE, returned for validation. ; ; sc - Scaling Coefficients, also called Father Coefficients, for ; the wavelet function. Useful if you want to get more detailed ; about the wavelet mathematics. ; ; wv - Wavelet Coefficients, also called Mother Coefficients, for ; the wavelet function. Useful if you want to get more detailed ; about the wavelet basis chosen. ; ; m2 - M-squared terms of the coherent component (see flatness ; calculation). ; ; m4 - M-fourth terms of the coherent component (see flatness ; calculation). ; ; flat - Flatness of the coherent component. ; ; flatnoi - Flatness of noise component. ; ; flattot - Flatness of original signal. ; ;Notes - ; The method requires that the processed data have a total number of ; points that is equal to a power of two. If your input array has a ; number of points that does not satisfy this, then a subsection of ; the input will be extracted for analysis. This subsection will ; begin with the first point of the input and keep as many points as ; possible. ; ; The variable, count1, is the exponent of two that is cycled to ; determine how many data points can be processed. You can change the ; initial value of this variable to effectively change the smallest ; array that can be processed. For example, the default beginning of ; two means that an array of 2^2 = 4 data points can be processed ; successfully. Increasing this initial value (especially if you ; always know the minimum size of the input) will make this part of ; the code run faster. ; ; The flatness and m-power terms are returned as a function of ; frequency. The frequency array is not yet calculate here because we ; are still testing different wavelet bases and the resultant ; frequency array is determined by the basis and any variable ; properties it may have. ; function wvint, data, \$ daub = daub ;constants acq = 100e6 / 64. ;sampling rate (Hz) ;get as many data points as possible, as a power of 2 nel = n_elements( data ) ;number of points given as input n = 2L ;variable will eventually be the number of data points processed count1 = 2. ;exponent, can increase to make computation faster ;determine the largest number of data points that satisfy the 2^n ;requirement while ( n LE nel ) do begin n = 2.^count1 count1++ endwhile j = count1 - 2 ;final exponent value freq = fltarr( j ) ;maximum frequency array (Hz) f = freq ;center frequency array (Hz) n = 2.^j ;total number of points to be analyzed ;reduce input to a power of 2 (number of data points) in = data[ 0 : n - 1 ] ;generate the wavelet depending on user's choice if keyword_set( daub ) then wvlt = wv_fn_daubechies( 24, sc, wv, scoff, wvoff ) else wvlt = wv_fn_coiflet( 5, sc, wv, scoff, wvoff ) nscales = n_elements( sc ) ;total number of wavelet scales waveOff = -1 * n / 2 + 2. ;offset of wavelet during transform ;perform the wavelet transform, centering the wavelet basis over each ;data point out = wv_dwt( in, sc, wv, waveOff, waveOff, /double ) ;flatness calculation of the total signal (calculated before adjusting ;the wavelet scales to remove the noise) m2tot = dblarr( j ) ;generate arrays to hold the results m4tot = m2tot ;also calculate wavelet frequency array within this loop for mm = 0, j - 1 do begin ;calculate m^2 terms m2tot[ mm ] = ( 1. / 2.^mm ) * total( out[ 0 : 2.^mm ]^2. ) ;calculate m^4 terms m4tot[ mm ] = ( 1. / 2.^mm ) * total( out[ 0 : 2.^mm ]^4. ) ;begin frequency work ;maximum frequency for each scale (Hz) freq[ mm ] = 2.^mm * ( acq / n ) endfor ;compute center frequencies f[ 0 ] = freq[ 0 ] / 2. for mm = 1, j - 1 do begin f[ mm ] = ( freq[ mm ] + freq[ mm - 1 ] ) / 2. endfor flattot = m4tot / ( m2tot^2. ) ;flatness of total (input) signal ;BEGIN section to find wavelet coefficients that correspond to ;coherent signal. var = [ total( abs( out )^2. ) / n ] ;compute variance (actually sigma^2) thresh = [ sqrt( 2. * alog(n) * var ) ] ;compute the initial threshold ;define the number of wavelet coefficients considered to be noise nnoise = n_elements( out ) ;all of them considered noise initially nold = 0 ;number of noise coefficients from previous iteration count = 0 ;placeholder for building arrays of threshold and variance ;loop progresses until the number of coefficients deemed noise does ;not change after an iteration while ( nold NE nnoise ) do begin nold = nnoise ;update the number for this iteration ;find coefficients that are below the threshold w = where( abs( out ) LE thresh[ count ] ) nnoise = n_elements( w ) ;update the number of noise coefficients var = [ var, total( abs( out[ w ] )^2. ) / n ] ;update variance array thresh = [ thresh, sqrt( 2. * alog( n ) * var[ count + 1 ] ) ] ;update threshold array count++ endwhile ;reconstruct the signal using the coherent part only out[ w ] = 0. ;all coefficients below the threshold are set to zero ;inverse wavelet transform constructs coherent signal out = wv_dwt( out, sc, wv, waveOff, waveOff, /double, /inverse ) ;noise component (by definition of being additive) noi = in - out ;BEGIN flatness calculation of coherent and noise components m2 = dblarr( j ) ;coherent variables m4 = m2 m2noi = m2 ;noise variables m4noi = m4 for mm = 0, j - 1 do begin ;m^2 terms calculated for both components m2[ mm ] = ( 1. / 2.^mm ) * total( out[ 0 : 2.^mm ]^2. ) m2noi[ mm ] = ( 1. / 2.^mm ) * total( ( noi[ 0 : 2.^mm ] )^2. ) ;m^4 terms calculated for both components m4[ mm ] = ( 1. / 2.^mm ) * total( out[ 0 : 2.^mm ]^4. ) m4noi[ mm ] = ( 1. / 2.^mm ) * total( ( noi[ 0 : 2.^mm ] )^4. ) endfor flatness = m4 / ( m2^2. ) ;flatness of coherent signal flatnoi = m4noi / ( m2noi^2. ) ;flatness of noise signal return, { raw:in, \$ ;raw input signal coh:out, \$ ;coherent part of the signal var:var, \$ ;variance of the wavelet coefficients thresh:thresh, \$ ;threshold for determining whether noise noi:noi, \$ ;noise part of the signal nnoise:nnoise, \$ ;number of wavelet coefficients deemed noise nold:nold, \$ ;same as NNOISE, returned for validation sc:sc, \$ ;scaling coefficients for wavelet function wv:wv, \$ ;wavelet coefficients for wavelet function m2:m2, \$ ;m-squared m4:m4, \$ ;m-fourth flat:flatness, \$;flatness of the coherent part flatnoi:flatnoi, \$ ;flatness of noise part flattot:flattot, \$ ;flatness of original signal freq:freq, \$ ;wavelet frequency (Hz) f:f } ;center frequency (Hz) end