;+ The webpage for this routine is, 
; http://www.davidpace.com/2015/11/12/farge-method/
; 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