/* file: multifractal.c Y. Ashkenazy 8 February 2004
Last revised: 14 October 2004 (GBM, IH)
-------------------------------------------------------------------------------
multifractal: Calculate the multifractal partition functions of a time series
Copyright (C) 2004 Yossi Ashkenazy
This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place, Suite 330, Boston, MA 02111-1307 USA.
You may contact the author by e-mail (ashkena@bgu.ac.il). For updates to this
software, please visit PhysioNet (http://www.physionet.org/).
_______________________________________________________________________________
*/
#include <stdio.h>
#include <math.h>
#include <values.h>
#include <stdlib.h>
#include <string.h>
#define SQR(x) ( (x) * (x) )
#define Min(x,y) ( ((x) < (y)) ? (x) : (y) )
#define Max(x,y) ( ((x) > (y)) ? (x) : (y) )
#define SIGN(x) (((x)>0.0)?(1):(-1))
#define SIGN1(x) (((x)<0.0)?(-1):(1))
#define EPS 1e-4
#define EPS1 1e-6
#define MAX_CHAR_IN_LINE 4096
#define STRL 80
#define PI M_PI
#define Pi2 2.0*PI
#define NumOfArgc 5 /* Number of arguments called by the executable file */
#define TimesWS 6 /* The factor of convolution range; i.e., wavelet
box =-TimesWS*WS..TimesWS*WS */
#define Ratio (1.0/8.0) /* Defines the ratio between maximum scale (wavelet
box) and signal length */
#define DQ 0.2 /* Resolution of moments q */
/* SUBROUTINE PrintColor
Prints the color cascade in a ppm format (can be read by xv or ImageMagick).
Subroutine for constructing color coded 3D wavelet decomposition:
x axis is time; y axis is wavelet scale; z axis is wavelet coefficient
(bright colors represent large values).
*/
PrintColor(double *vec, long i, long nx, long ny, double max, double min)
{
long l, j, bin, color1, color2, color3, N, N1, N2, N3, N4, N5, N6;
double delta;
N = 255;
N1 = N;
N2 = 2*N;
N3 = 3*N;
N4 = 4*N;
N5 = 5*N;
N6 = 6*N;
delta = (max - min)/(double)(N6-2);
printf("P3\n# CREATOR: multifractal\n%d %d\n255\n", nx, ny);
j = l = 0;
while (l < i) {
bin = (long)(((vec[l] - min)/delta));
if (bin < N1) {
color1 = N;
color2 = bin;
color3 = 0;
}
else if (bin < N2) {
color1 = N2-bin;
color2 = N;
color3 = 0;
}
else if (bin < N3) {
color1 = 0;
color2 = N;
color3 = bin-N2;
}
else if (bin < N4) {
color1 = 0;
color2 = N4-bin;
color3 = N;
}
else if (bin < N5) {
color1 = bin-N4;
color2 = 0;
color3 = N;
}
else {
color1 = N;
color2 = 0;
color3 = N6-bin;
}
printf("%3d %3d %3d ", color1, color2, color3);
l++;
if (++j >= 5) {
printf("\n");
j = 0;
}
}
}
/* SUBROUTINE F
Calculates continuous Gaussian wavelet functions (zero to 7th derivative).
*/
double F(long n, double t)
{
switch (n) {
default:
case 0: return exp(-0.5*SQR(t));
case 1: return -t*exp(-0.5*SQR(t));
case 2: return (-1+t*t)*exp(-0.5*SQR(t));
case 3: return t*exp(-0.5*SQR(t))*(3-t*t);
case 4: return exp(-0.5*SQR(t))*(pow(t,4.0)-6*t*t+3);
case 5: return -t*exp(-0.5*SQR(t))*(pow(t,4.0)-10*t*t+15);
case 6: return exp(-0.5*SQR(t))*(pow(t,6.0)-15*pow(t,4.0)+45*t*t+15);
case 7: return -t*exp(-0.5*SQR(t))*(pow(t,6.0)-21*pow(t,4.0)+105*t*t-105);
}
}
/* SUBROUTINE PrepareVecF
Vector storing the values of the wavelet function; done for numerical
efficiency
*/
void PrepareVecF( double ws, double vec[], long Derivative )
{
long i, tempi = (long)((double)TimesWS*ws);
for (i = 0; i <= 2*tempi; i++)
vec[i] = F((long)Derivative,((double)i-(double)(tempi))/ws);
}
/* SUBROUTINE MF1
Find the maximum for the scale window ws, by performing the following steps:
1) Wavelet convolution of the signal for increasing wavelet scale.
2) Locate the local maxima of the absolute value of wavelet
coefficient as a function of time for each wavelet scale.
3) Check whether a local maximum at a given wavelet scale is located
close to a maximum at a smaller scale - if yes connect both maxima,
otherwise cancel it. Generate maxima lines.
4) Check that the number of maxima at larger scales is less or equal
to that at a smaller scale.
5) Track maxima lines for increasing wavelet scale by choosing at each
scale value the supremum between all previous values at smaller
scales.
*/
double MF1(double *vec, long n, double min_q, double max_q, double dq,
double min_ws, double max_ws, double dws, long Derivative, long code)
{
long nq,i,i1,i2,j,jj,sign,sign1,tempi_max,tempi,tempi1,tempi2,n_max,icolor;
double *s, s1, min, max, temp, temp1, temp2, q;
double *MaxVec, *VecWL, *VecF, *VecColor, *pt, ws;
long *VecWLX, *MaxVecX, *ptl;
nq = (long)((max_q-min_q)/dq)+1;
MaxVec = (double *)calloc((size_t)(n+1),(size_t)(sizeof(double)));
MaxVecX = (long *)calloc((size_t)(n+1),(size_t)(sizeof(long)));
tempi_max = (long)((double)TimesWS*max_ws);
s = (double *) calloc( (size_t)(nq+1),sizeof(double) );
VecF = (double *) calloc( (size_t)(2*tempi_max+1),sizeof(double) );
VecWL = (double *) calloc( (size_t)(n+1),sizeof(double) );
VecWLX = (long *) calloc( (size_t)(n+1),sizeof(long) );
if (code == 1) {
printf("# %f\n", min_q);
fflush(NULL);
}
else if (code == 3) {
for (ws = min_ws, tempi1 = 0; ws <= max_ws; ws *= dws)
tempi1++;
icolor = 1;
max = -1e20;
min = 1e20;
VecColor = (double *)calloc((size_t)((n-2*tempi_max-1)*tempi1+1),
sizeof(double));
}
jj=0;
for (ws = min_ws; ws <= max_ws; ws *= dws) {
tempi = (long)((double)TimesWS*ws);
for (i = 0; i < (2*tempi+1); i++)
VecF[i] = 0.0;
for (i = 0; i < (n+1); i++)
VecWL[i] = 0.0;
for (i = 0; i < (n+1); i++)
VecWLX[i] = 0;
PrepareVecF(ws, VecF, Derivative);
/* do the convolution */
for (i = tempi_max; i < (n-tempi_max-1); i++) {
s1 = 0.0;
for (j = i-tempi; j <= (i+tempi); j++)
s1 += vec[j] * VecF[j-i+tempi];
VecWL[i] = fabs(s1/ws);
if (code == 3) {
VecColor[icolor] = VecWL[i];
icolor++;
if (max < VecWL[i])
max = VecWL[i];
if (min >VecWL[i])
min = VecWL[i];
}
}
/* Find the local maxima of the wavelet coefficient for each scale */
n_max = 0;
sign = SIGN(VecWL[tempi_max+1] - VecWL[tempi_max]);
temp1=0.0;
for (j = 0, i = tempi_max+2; i < (n-tempi_max-1); i++) {
if ((fabs(VecWL[i] - VecWL[i-1]) > 0.0) &&
((sign == 1) && ((SIGN(VecWL[i] - VecWL[i-1])) == (-1)))) {
temp = VecWL[i-1];
if (code == 2) { /* Print the maxima lines (code 2) */
printf("%d %g\n", i, log(ws)/log(10.0));
fflush(NULL);
}
n_max++;
VecWL[j]=temp;
VecWLX[j]=i-1;
temp1=temp;
j++;
}
sign=SIGN(VecWL[i]-VecWL[i-1]);
}
/* Tracking the maxima lines: test for supremum */
if (ws > min_ws) {
for (i1 = i2 = i = 0; i < nq;i++)
s[i]=0.0;
while (((i1-1) < n_max) && ((i2-1) < jj)) {
if ((VecWLX[i1] - MaxVecX[i2]) <= (MaxVecX[i2+1] - VecWLX[i1]))
VecWL[i1] = Max(VecWL[i1],MaxVec[i2]);
else
VecWL[i1] = Max(VecWL[i1],MaxVec[i2+1]);
for (i = 0; i < nq; i++)
s[i] += pow(VecWL[i1], min_q+(double)i*dq);
i1++;
i2++;
while ((i2 < jj) && (VecWLX[i1] >= MaxVecX[i2]))
i2++;
i2--;
}
if (code == 1) { /* Print the partition function */
printf("%g ", log(ws)/log(10.0));
for (i = 0; i < nq; i++)
printf("%g ", log(s[i])/log(10.0));
printf("\n");
fflush(NULL);
}
}
jj = n_max;
pt = MaxVec;
MaxVec = VecWL;
VecWL = pt;
ptl = MaxVecX;
MaxVecX = VecWLX;
VecWLX = ptl;
}
free(VecF);
free(VecWL);
free(VecWLX);
free(MaxVec);
free(MaxVecX);
/* Print the wavelet cascade (code 3) */
if (code == 3) {
PrintColor(VecColor,icolor,n-2*tempi_max-1,tempi1,max,min);
free(VecColor);
}
}
main(int argc, char *argv[])
{
long i, j;
long N, Derivative = 3, code = 1;
double *VecX, temp, tempx, tempy, ws, dw, MinScale, MaxScale,
q, q_min, q_max, dq;
FILE *fp;
if (argc < NumOfArgc) {
printf("Usage: %s INPUT N QMIN QMAX DW MODE >OUTPUT\n", argv[0]);
printf(" INPUT name of file containing the input time series\n");
printf(" N number of points (lines) in INPUT\n");
printf(" QMIN minimum MF order\n");
printf(" QMAX maximum MF order\n");
printf(" DW order of the Gaussian derivative wavelet (0-7)\n");
printf(" MODE the type of output to be produced, one of:\n");
printf(" 1: partition functions (text)\n");
printf(" 2: maxima lines (text)\n");
printf(" 3: wavelet cascade (PPM image)\n\n");
printf(" INPUT is a text file containing two columns of numbers; the\n");
printf(" first is ignored, and the second contains the data values.\n");
exit(0);
}
/* N is the series length */
N = atol(argv[2]);
/* VecX is the vector storing the data */
VecX = (double *) malloc((size_t)(N+1) * sizeof(double));
for (i = 0; i < N+1; i++)
VecX[i] = 0.0;
/* Open data file for reading two column data file. The second column
should be the data that are analyzed. */
fp = fopen( argv[1], "r");
for (i = 0; (i < N) && (fscanf(fp, "%lf%lf", &temp, &(VecX[i])) == 2); i++)
;
fclose(fp);
N = i;
/* MaxScale is the maximum scale analyzed */
MaxScale = Ratio*0.5*(((double)N)-1.0)/(double)TimesWS;
/* MinScale is the minimum scale analyzed */
MinScale = 2.0; /* 0.5*M_E */
/* q_min is the minimum moment q for which the partition function is
calculated */
q_min = atof(argv[3]);
/* q_max is the maximum moment q for which the partition function is
calculated */
q_max = atof(argv[4]);
/* resolution of moments in steps dq for which the partition function is
calculated */
dq = DQ;
if (argc > NumOfArgc) {
Derivative = atoi(argv[5]);
if (argc > (NumOfArgc+1))
code = atoi(argv[6]);
}
/* scale resolution for which the partition function is calculated */
dw = pow(2.0,0.05);
MF1(VecX, N, q_min, q_max, dq, MinScale, MaxScale, dw, Derivative, code);
free(VecX);
}