;+ 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