multifractal.c

计算时间序列的多分形程序

/* 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 ([e-mail]).  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);
}