##..........................................................................## ## ## ## Script functions in Lastwave language for computing the multifractal and ## ## singularity spectra of binary data ## ## ## ##..........................................................................## if ([var exist DISPLAY] == 0) { DISPLAY = [new &int] DISPLAY = 1 } if ([var exist BINARY] == 0) { BINARY = [new &int] BINARY = 0 } if ([var exist QLIST] == 0) { QLIST=[new &listv] QLIST = {-4. -3.6 -3.2 -3. -2.8 -2.6 -2.4 -2.2 -2. -1.8 -1.6 -1.4 -1.2 -1.1 -1. -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1. 1.05 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2. 2.3 2.6 3. 3.5 4. 5. 6. 7. 8.} } ## ## Get the WT and the extrema representation of a given signal ## setproc WTransExtCompute {{&signal s} {&num aMin 1.5} {&num nOct 8.} {&num nVoice 10.} \ {&string wavelet 'g2'}} \ "{{{[] [=1.5] [=8] [=10] [='g2']} \ {WTransExtCompute computes the wavelet transform and the extrema representation of the 1d-signal . Results are returned in a wtrans variable .}}}" { global DISPLAY w=[new &wtrans] # Copy the original signal in a wtrans structure 0w = s # copy s 0w echo Compute the wavelet transform ## Compute WT of the cantor using the th derivative of the gaussian cwtd w aMin nOct nVoice wavelet -e 0 :: >> ## ## cwtd [=objCur] [='g2'] [-b = ## mir] [-c] [-m] [-e = -1] ## Performs the continuous wavelet transform of the signal A[0][0] in ## from scale using octaves, voices per octave, using ## the nth derivative of the Gaussian function (=g where n is ## in [0..4] are such that g0, -g1, g2, -g3, g4 are the gaussian function ## and its successive derivatives) as analyzing wavelet or the morlet wavelet ## (which is a complex wavelet). ## note: the commmand cwtd returns a time elapsed in seconds: who cares... that's ## the reason of the redirection :: >> echo Compute the extrema representation from the wavelet decomposition ## compute the extrema representation extrema w w.extrep # -i -C ## ## extrema [=objCur] [=] [-icC] [-e ] ## Computes the extrema representation associated to the wavelet transform ## and puts it in variable . It returns the number of extrema ## found. Let us note that this command leads to unstable results at large ## sacles (i.e., too many extrema). In order to make the computation stable, ## by default, the 'extrema' command chain the extrema (using the 'chain' ## command) and delete the extrema (using the 'chaindelete' command) which ## have not been chained. ## -i : The extremum positions are computed without using interpolation and ## thus always correspond to a sample of the original signal ## -c : Causal effects are taken into account (The extrema are computed only ## between the 'firstp' and 'lastp' indexes of each wavelet coefficient signal). ## -C : Disable computation of chains and deletion of extrema which do not ## belong to the chains. ## -e : Sets the below which no maxima will be considered. ## If we want to display the extrema representation if (DISPLAY == 1) {disp Extrema w.extrep -pos 700 60 -size 300 300} return w } ## ## Compute the PF and the different multifractal spectra associated to a given WT ## setproc PFCompute {{&wtrans w} \ {&num aMin 1.5} {&num nOct 6.} {&num nVoice 8.} \ {&listv qlist {-200.}} \ {&num logaMin -200.} {&num logaMax -200.}} \ "{{{[] [=1.5] [=6] [=8.] [=QLIST] [=log2()] [=+]} \ {PFCompute computes the (singularity and multifractal) spectra associated to the WTMM representation stored in the variable . The returned result is a list { }, where is a partition function, and and are two signals storing the multifractal and the singularity spectra resp.}}}" { # note : -200 => absurd value global DISPLAY pf=[new &PF] if (qlist.size == 1 && qlist[0] == -200) { qlist = QLIST } echo Compute the partition function on the maxima lines i.e. the WTMM method ## compute the partition functions pf wtmm pf w.extrep qlist ## ## pf wtmm [] ## Computes the partition functions using the extrema representation. It is ## the so-called WTMM method. It returns a new partition function if none ## were specified on the command line or itself. ## if we want to display ALL the computed 't' partition function # echo display tau(q) for all the scales if (DISPLAY==1) { tauq = [pf get t pf] disp PF {tauq} -pos 381 395 -size 300 300 \ -..fv1 -axisFont 'fixed-10-plain' -yLabel "log2 Z(q,a)" -xLabel "log2 a" } if (logaMin == -200.) {logaMin=log2(aMin)} if (logaMax == -200.) {logaMax=logaMin+nOct} ## The linear regression of the previous partition function let us access the ## tau(q) spectrum. ## compute and display the multifractal exponents tauq = [myTauqSpectrum pf logaMin logaMax] if (DISPLAY == 1) {disp Tau_q tauq -pos 700 395 -size 300 300 \ -..fv1 -yLabel "tau(q)" -xLabel "q" -bound -7 9 -7 5 \ -..1 -curve 'o' -3 -fg 'blue'} ## The linear regression of the previous partition functions let us access the ## D(h) spectrum. ## compute and display the singularity spectrum dh = [mySingSpectrum pf logaMin logaMax] if (DISPLAY == 1) {disp D_h dh -pos 60 60 -size 300 300 -x 0 1 -y 0 1.5 \ -..fv1 -title "D(h)" -..title -font 'fixed-10-plain' \ -..1 -curve 'o' 3 -fg 'red'} return {pf tauq dh} } # # # setproc AnalyzeSeriesLastWTMM {{&string storedir} {&string savedir} {&string file} {&int N 1} \ {&num aMin 1.5} {&num nOct 8.} {&num nVoice 10.} \ {&string wavelet 'g2'} \ {&listv qlist {-200}} \ {&num logaMin -200} {&num logaMax -200}} \ "{{{[] [] [] [=1] [=1.5] [=8] [=10] [='g2'] [=QLIST] [=log2()] [=+]} \ {AnalyzeSeriesWTMM computes the (singularity and multifractal) spectra associated to the WTMM representation of a serie of signals stored in the text or binary files with generic name -N???.dat in directory . The returned result is a list { }, where is a partition function, and and are two signals storing the multifractal and the singularity spectra resp. Saves the results in files Tauq_lw_ and Dh_lw_ of directory .}}}" { # note : -200 => absurd value global DISPLAY global BINARY global QLIST w=[new &wtrans] pf=[new &PF] if (qlist.size == 1 && qlist[0] == -200) { qlist = QLIST } qldisp = [sprintf toto "%.1f" qlist[0]] foreach i 1:qlist.length-1 {qldisp += [sprintf toto " %.1f " qlist[i]]} # Loop on the trials foreach i 0:N-1 { if (BINARY == 1) { sprintf infile "%s/%s-N%03d.fdat" storedir file i read 0w infile -r } else { sprintf infile "%s/%s-N%03d.txt" storedir file i read 0w infile } ## ## read ( | ) [[] ] [-f ] ## [-s ] [-r [('little' | 'big')]] [-b] ## Reads a signal from a or a . The file must be in the LW ## format (i.e., created with 'write') or in the raw format (option '-r'). ## If option '-r' is not set then you can specify the column number for the ## x-values and the column number for the y-values (first column is 1). If ## none are specified then this command tries to read the first two columns ## as x and y or, if there is only one colum, directly the y values. If a ## stream is specicified and if option '-r' is not set then, at the end of ## the command the stream is positionned at the end of the signal (even if ## just a few points are read). The options are ## -s : It reads only values ## -r : The data (just the y\'s) are in raw format (binary floats) with ## no header. ## compute WT of the signal using the th derivative of the gaussian cwtd w aMin nOct nVoice wavelet :: >> ## compute the extrema representation: see above extrema w # w.extrep ## Use the scale adaptive method to chain the maxima lines ## Replace all the extrema ordinates by the max ordinate on the extrema line ## (finer scale). It takes into account a factor for computing the max value chainmax w.extrep 1 ## Replaces the wavelet trasnform value of an extremum by the maximum value ## along its chain. This routine allows to perform the scale adaptative ## version of the WTMM method. The coefficients are first multiply by ## aMin*(2^(exponent*(o-1+v/nvoice))). This is useful in the case some maxima ## chains do not decay when going to small scales. ## compute the partition functions if (i == 0) { pft = [pf wtmm w.extrep qlist] } else { pft = [pf add pft [pf wtmm w.extrep qlist]] } } copy pft pf printf "\n" ## Display ALL the computed 'd' and 'h' partition function ## display h(q) and D(q) if (DISPLAY == 1) {disp Spectrum {[pf get h pf]} {[pf get d pf]} -pos 460 60 -size 300 600 \ -..fv1 -title "H(q=$qldisp)" \ -..fv2 -title "D(q=$qldisp)" \ -..title -font 'fixed-10-plain'} if (logaMin == -200.) {logaMin=log2(aMin)} if (logaMax == -200.) {logaMax=logaMin+nOct} ## Other methods tauq = [myTauqSpectrum pf logaMin logaMax] if (DISPLAY == 1) {disp Tau_q tauq -pos 700 395 -size 300 300 \ -..fv1 -yLabel "tau(q)" -xLabel "q" -bound -7 9 -7 5 \ -..1 -curve 'o' -3 -fg 'blue'} dh = [mySingSpectrum pf logaMin logaMax] ## Display D(h) if (DISPLAY == 1) {disp D_h dh -pos 60 60 -size 300 300 -x 0 1 -y 0 1.5 \ -..fv1 -title "D(h)" -..title -font 'fixed-10-plain' \ -..1 -curve 'o' 3 -fg 'red'} ## write ( | ) [('xy' | 'yx' | 'x' | 'y']) [-h] [-b] ## [-r [('little' | 'big')]] ## Writes a signal into a or a . It writes the file either ## in a raw format (using '-r' option) or (by default) in the LastWave format ## with ascii coding (unless '-b' is specified in which case binary codes are used). ## In the LastWave format, there are different modes : ## - 'xy' : First column is X and second is Y (in ascii mode), or the X-values ## are stored first then the Y-values (in binary mode). ## - 'yx' : First column is Y and second is X (in ascii mode), or the Y-values ## are stored first then the X-values (in binary mode). ## - 'y' : A single column made of Y ## - 'x' : A single column made of X ## By default it uses 'xy' mode for xy-signals and 'y' mode for y-signals. ## The options are : ## -h : No Header (ONLY in ascii mode). ## -b : Binary writing (instead of ascii). ## -r : The data (just the y\'s) are written in raw format (binary floats) ## sprintf taufile "%s/Tauq_lw_%s" savedir file write tauq taufile 'xy' -h sprintf dfile "%s/Dh_lw_%s" savedir file write dh dfile 'xy' -h return {pf tauq dh} } ## ## Script functions to automatically get the singularity - graph (h,D(h)) - and ## the multifractal - graph (q,tau(q)) - spectra from a PF structure. ## ## Note : mySingSpectrum and myTauqSpectrum are nothing else than the functions ## singSpectrum and tauqSpectrum of the file wtmm1d.pkg. They have just been copied ## here for convenience ## ## ## Get the D(h) spectrum from a PF structure. ## setproc mySingSpectrum {{&PF pf} {&float log2aMin} {&float log2aMax}} \ "{{{ } \ {Computes the singularity spectrum D(h) from the partition function structure . It reads the partition function values for both h(q) and D(q) for all the available q's then performs automatic fits of the partition functions between the log scale values and in order to obtain the h(q) and D(q) values. The D(h) function is stored in the (xy) signal . The fits are done automatically between and .}}}" \ { ## Init the result signal as an XY empty signal result = XY(Zero(pf.qnumber),Zero) ## Init the loop i = 0 ## main Loop foreach q pf.qlist { ## Get the h(q) value and sets it in result (x array) result.X[i] = [stats fit [pf get h pf q] -x log2aMin log2aMax][0] ## Get the D(q) value and sets it in result (y array) result.Y[i] = [stats fit [pf get d pf q] -x log2aMin log2aMax][0] i+=1 } ## sort according to x-array sort result return result } ## ## Get the tau(q) spectrum from a PF structure. ## setproc myTauqSpectrum {{&PF pf} {&float log2aMin} {&float log2aMax}} \ "{{{ } \ {Computes the tau(q) spectrum from the partition function structure . It reads the partition function values for tau(q) for all the available q's then performs automatic fits of the partition functions between the log scale values and in order to obtain the tau(q) values which are stored in . The fits are done automatically between and .}}}" \ { result = XY(Zero(pf.qnumber),Zero) ## Init the loop i = 0 ## main Loop foreach q pf.qlist { ## Set the q result.X[i] = q ## Get the T(q) value and sets it in result (y array) result.Y[i] = [stats fit [pf get t pf q] -x log2aMin log2aMax][0] i+=1 } return result }