/* 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 #include #include #include #include #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); }