chaos_analyser.c

混沌分析工具

/*
 * The whole analyser : SVD reduction
                        LSVD analysis
			Time domain plot
			Synthesiser (for speech only):
			      Normaliser
			      start/mid/end morphing
			      Across vowel morphing
			Poicare Sections
			Lyapunov exponents
			
 
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <X11/X.h>
#include <X11/Xlib.h>
#include <X11/Xos.h> 
#include <xview/cms.h>
#include <xview/xview.h>
#include <xview/canvas.h>
#include <xview/panel.h>
#include <xview/xv_xrect.h>
#include <xview/notify.h>
/*#include <X11/bitmaps/xlogo64>*/
#include <xview/svrimage.h>
#include <xview/icon.h>
#include "Recipes/read_double.c"
#include "Recipes/gasdev.c"
#include "Recipes/ran1.c"
#include "Recipes/avevar.c"
#include "Recipes/rotated.h"
#include "Recipes/rotated.c"
#include "Recipes/fig_drawing.c"

/* Color indices */
#define WHITE           0
#define RED             1
#define GREEN           2
#define BLUE            3
#define GREY            4
#define BLACK           5
#define NUM_COLORS      6

#define PI 3.14159
#define PIbyTWO PI/2.0
#define MAX_rows 51000
#define LABEL_WIDTH 160

#include "lor_icon.xbm"



/*******  general global variables *****************************************/
void    event_proc(), repaint_proc(), redraw_proc();
void    newfile(),event_proc(),set_angles(),readfile(),
        save_to_file(),get_print_name(),print_attr(),print_mut();
Frame   frame,
        mes_frame,
        options_frame,
        print_frame;
Display *dpy;
GC      gc,gc_dash;
Window  canvas_win;
Notify_value     step();
struct itimerval timer;
Canvas           canvas;
unsigned long *pixel_table; /* pixel values for colors */
Pixmap        xlogo;        /* the xlogo */
Panel_item M_SLIDER,
           SVD_SLIDER,
           POS_SLIDER,
           NUM_POINTS,
           FR_BUT,
           SVD_ONOFF,
           SPEED,
           PRINT_MESSAGE,
           AXES_ON_OFF;
char    kbd_msg[128], ptr_msg[128], but_msg[128],FILENAME[50],
        PRINT,PRINT_NAME[50];
double  **DATA,
        SIZE=0.5,
        MAX_VAL=0.000001,
        ALPHA=0.7854,BETA=0,GAMMA=5.6723,
        T1,T2,T3,T4,T5,T6,T7,T8,T9;
int     LENGTH=150,
        STEP=1,
        X_SET=250,
        Y_SET=250,
        GHOST_LENGTH=200,
        M=10,
        MULT=100,
        NO_POINTS=5000, /* number of points in main file */
        COUNTER=1,
        SVD_COLS=10,
        SVD=1,
        FREEZE=-1,
        SPEECH=0,
        COUNT_LAST=1,
        OPTIONS=-1,
        AXES=1;

/*******   Local Singular Value Global Variables ***************************/
void    lsvd_paint_proc(),print_lsvd();
Frame   lsvd_frame;
Window  canvas_win_lsvd;
Canvas  canvas_lsvd;
int     LSVD_NUM=20,
        LSVD_LENGTH=1000,
        LSVD_POS;
double  LSVD_MULT=1,
        LSVD_ST_RAD=100,
        svd[50][7],
        lsvd_rad[50];


/******   Time Domain Global Variables *************************************/
void    time_paint_proc();
Frame   time_frame;
Window  canvas_win_time;
Canvas  canvas_time;
double  TMULT=1,
        TSIZE=1;


/******   Poincare section global variables *******************************/
void poin_paint_proc();
Frame  poin_frame;
Canvas canvas_poin;
Window canvas_win_poin;
double P_SIZE;
int    POIN_ON=1,
       X_SET_P=150,
       Y_SET_P=150;

/******* mutual information section ***************************************/
void mut_paint_proc();
Frame mut_frame;
Window canvas_win_mut;
Canvas canvas_mut;
double MUT_SIZE=1,
       MUT_RAD=100,
       mut[200],
       mut_delay[200];
int    MUT_DELAY=100,
       MUT_NO_POINTS=1000,
       MUT_NUM=100;

/*******   Lyapunov expoent Global Variables ***************************/
void    lya_paint_proc(),lya_setup_proc();
Frame   lya_frame;
Window  canvas_win_lya;
Canvas  canvas_lya;
Panel_item SVD_M_SLIDER,
           LYA_EMBEDDING,
           LYA_START_STOP;
int     B_VECT=20,            /* number of vectors in B         */
        B_AVE=1,             /* number of vectors in R         */
        REINIT=5,            /* steps between reinitialisation */
        GLOBAL_DIM=3,        /* global embedding dimension     */
        LOC_DIM=3,           /* embedding dimension            */
        EV_POS,
        EVOLVE=1000,
        LYA_STOP=1,
        LYA_ZERO=150,
        SVD_WIN=20;           /* svd window length              */
double  **LYA_EXP,            /* cummulative sum of exponents   */
        TAU=1,
        LYA_XSCALE=1,
        LYA_YSCALE=10;
char    EMBED='t',         /* type of embedding (s/t/a)      */
        DIAG,              /* diagnostics flag (y/n)         */
        CENTRE,            /* difference or constant centre  */
        LOCAL;             /* local tangent reduction ?      */

#include "lya_section.c"

/***********************************************************************

                  Local Singular Value procedures

***********************************************************************/
void
lsvd_print(item,value)
     Panel_item item;
     int        value;
{
  xv_set(print_frame,XV_SHOW,TRUE,NULL);
  
  PRINT='s';
}

void localsvd(E,matrix_array,num_row,num_col)
     double **matrix_array,*E;
     int   num_row,num_col;

{
double **matrix_transpose;
double **corr,**tmp1;
int   i,j,nrot;


/************ initialise all matrices  ****************************/
matrix_transpose=dmatrix(1,num_col,1,num_row);
corr=dmatrix(1,num_col,1,num_col);
tmp1=dmatrix(1,num_col,1,num_col);

/*************  calculate C matrix ***********************************/
transpose(matrix_array,matrix_transpose,num_row,num_col);
multiply_matrix(tmp1,matrix_transpose,matrix_array,num_col,num_row,num_row,num_col);
jacobi(tmp1,num_col,E,corr,&nrot);
eigsrt(E,corr,num_col);

free_dmatrix(tmp1,1,num_col,1,num_col);
free_dmatrix(matrix_transpose,1,num_col,1,num_row);
free_dmatrix(corr,1,num_col,1,num_col);


}

void obtain_R(X,R,num,row,col,radius,m,pos)
     double **X,*R,*radius;
     int   *num,row,col,m,pos;

{
  int i,j,count=0;
  double distance;

  while(count<=col){
    (*radius) *=1.25;
    count=0;
    i=0;
    while (i<=row-10){
      i++;
      distance=0.0;
      j=0;
      while (distance <= *radius && j<col){
	j++;
	distance=fabs(X[pos][j]-X[i][j]);
      }
      if (distance<= *radius &&  i != pos){
	count++;
	R[i]=1;
      }
      else R[i]=0;
    }
    *num=count;
  }
}


void obtain_B(X,R,B,row,col,k)
     double *R,**B,**X;
     int   k,row,col;
     
{
  int i,j,count;
  double distance,test;

  count=0;
  for (i=1;i<=row;i++){
    if (R[i]==1){
      count++;
      for (j=1;j<=col;j++){
	B[count][j]=X[LSVD_POS][j]-X[i][j];
      }
    }
  }
}


void readfile_lsvd(dataarray,window,m,num)
     double  **dataarray;
     int    window,m,*num;

{
  double databit,buffer[100];
  FILE *ioptr;
  int i,l,j,check,sect=1;
  char data[256],space;

  ioptr = fopen(FILENAME,"r");

  if (ioptr == (FILE *)NULL) {
    printf("\nCannot open file !");
    exit(1);
  }
 
  for(l=1;l<=m* window;l++){
    read_double(ioptr,&databit);
    buffer[l]=databit;
  }
  j=0;
  for (i=1;i<COUNTER-LSVD_LENGTH;i++){
    j++;
    read_double(ioptr,&databit);
  }
  
  i=0;
  while (read_double(ioptr,&databit)==0 && j<COUNTER+LSVD_LENGTH){
    i++;
    j++;
    if (j==COUNTER)
      LSVD_POS=i;
    for (l=0;l< window;l++)
      dataarray[i][l+1]=buffer[l* m +1];
    for (l=1;l<=m* window -1;l++)
      buffer[l]=buffer[l+1];
    buffer[m* window]=databit;
  }
   *num=i;
  fclose(ioptr);
}


void loc_svd()
{
  char  reduced;
  double **X,**B,*E,*R;
  double radius,mult=1.25,start_radius;
  int    window=6,loop,j,B_elements=20,i;
  int    evolve,num;
  FILE   *optr;


/************* initialise matrices, vectors and variables *****************/
  X=dmatrix(1,2*LSVD_LENGTH+10,1,window);
  B=dmatrix(1,2*LSVD_LENGTH+10,1,window);
  E=dvector(1,window);
 
  R=dvector(1,2*LSVD_LENGTH+10);
  num=2*LSVD_LENGTH;
  readfile_lsvd(X,window,M,&num); 

  radius=LSVD_ST_RAD;
   
/****** main loop for evolving around the attractor *******************/ 
  for(evolve=1;evolve<=LSVD_NUM;evolve++){
 
    obtain_R(X,R,&B_elements,num,window,&radius,M,LSVD_POS);

    obtain_B(X,R,B,num,window,M);  

    localsvd(E,B,B_elements,window);

    for (i=1;i<=window;i++)
      svd[evolve][i]=log(fabs((1/sqrt(B_elements)) * sqrt(fabs(E[i]))));

    lsvd_rad[evolve]=log(radius);
  }
  free_dmatrix(X,1,2*LSVD_LENGTH+10,1,window);
  free_dmatrix(B,1,2*LSVD_LENGTH+10,1,window);
  free_dvector(E,1,window);

}


void
lsvd_on_off(item,value)
     Panel_item item;
     int value;
{
  if (value==1)
    xv_set(lsvd_frame,XV_SHOW,FALSE,NULL);
  else{
    xv_set(lsvd_frame,XV_SHOW,TRUE,NULL);
   }
}


lsvd_num_adjust(item,event)
     Panel_item item;
     Event *event;
{
  LSVD_NUM=xv_get(item,PANEL_VALUE);
}

lsvd_st_rad(item,event)
     Panel_item item;
     Event      *event;
{
  char *test[50],str[50];
 
  strcpy(str,(char *)xv_get(item,PANEL_VALUE));

  LSVD_ST_RAD=strtod(str,test);
}

lsvd_length(item,event)
     Panel_item item;
     Event *event;
{
  LSVD_LENGTH=xv_get(item,PANEL_VALUE);
}

  
lsvd_accept(item,value)
     Panel_item item;
     int value;
{
  Xv_Window window2=canvas_paint_window(canvas_lsvd);

  loc_svd();
  lsvd_paint_proc((Canvas)NULL,window2,
        dpy,
        canvas_win_lsvd, (Xv_xrectlist *) NULL); 
}


void
lsvd_paint_proc(canvas_lsvd,window, dpy, xwin, xrects)
Canvas        canvas_lsvd;           /* Ignored */
Xv_Window    window;
Display      *dpy;
Window        xwin;
Xv_xrectlist *xrects;           /* Ignored */

{
 int i,j,start,length;
 double MAX_rad,min_rad,MAX_svd,min_svd,ymult,xmult;
 char msg[50];
 

 MAX_rad=lsvd_rad[LSVD_NUM];
 min_rad=lsvd_rad[1];
 MAX_svd=svd[LSVD_NUM][1];
 min_svd=svd[1][6];
 xmult=200/(MAX_rad-min_rad);
 ymult=200/(MAX_svd-min_svd);
 XSetForeground(dpy, gc, pixel_table[BLACK]); 
  
 XClearWindow(dpy, xwin);
 XDrawLine(dpy,xwin,gc,40,260,290,260);
 XDrawLine(dpy,xwin,gc,40,260,40,10);
 sprintf(msg,"%3.1lf",min_rad);
 XDrawString(dpy,xwin,gc,60,280,msg,strlen(msg));
 sprintf(msg,"%3.1lf",MAX_rad);
 XDrawString(dpy,xwin,gc,260,280,msg,strlen(msg)); 
 sprintf(msg,"%3.1lf",MAX_svd);
 XDrawString(dpy,xwin,gc,10,20,msg,strlen(msg));
 sprintf(msg,"%3.1lf",min_svd);
 XDrawString(dpy,xwin,gc,10,220,msg,strlen(msg));
 sprintf(msg,"log (Radius)");
 XDrawString(dpy,xwin,gc,120,280,msg,strlen(msg));
sprintf(msg,"SVD");
 XDrawString(dpy,xwin,gc,10,120,msg,strlen(msg));
 
 for (i=1;i<=6;i++){
    for (j=1;j<LSVD_NUM;j++){
      XDrawLine(dpy,xwin,gc,60+(lsvd_rad[j]-min_rad)*xmult,220-(svd[j][i]-min_svd)*ymult,60+(lsvd_rad[j+1]-min_rad)*xmult,220-(svd[j+1][i]-min_svd)*ymult);
    }
  }
}

void
print_lsvd(filename)
     char filename[100];
{
 int i,j,start,length,p,interval;
 double MAX_rad,min_rad,MAX_svd,min_svd,ymult,xmult;
 char msg[100],fig_file[100];
 FILE *optr;


 MAX_rad=lsvd_rad[LSVD_NUM];
 min_rad=lsvd_rad[1];
 MAX_svd=svd[LSVD_NUM][1];
 min_svd=svd[1][6];
 
 xmult=4000/(MAX_rad-min_rad);
 ymult=4000/(MAX_svd-min_svd);
 sprintf(fig_file,"%s.gnu_setup",filename);
 optr = fopen(fig_file,"w");

 fprintf(optr,"set xlabel \"ln (Radius)\"\n");
 fprintf(optr,"set ylabel \"ln (Singular value / sqrt(N))\"\n");
 fprintf(optr,"set term post port 24\n");
 fprintf(optr,"set output \"%s.ps\"\n",filename);
 fprintf(optr,"plot \"%s.gnu_data\" notitle w l\n",filename);
 fprintf(optr,"set term x11\n");
 fclose(optr);
 sprintf(fig_file,"%s.gnu_data",filename);
 optr = fopen(fig_file,"w");
 
 for (i=1;i<=6;i++){
   for (j=1;j<=LSVD_NUM;j++){
     fprintf(optr,"%lf %lf\n",lsvd_rad[j],svd[j][i]);
   }
   fprintf(optr,"\n");
 }
 fclose(optr);

 sprintf(msg,"gnuplot < %s.gnu_setup",filename);
 system(msg);
 sprintf(msg,"rm %s.gnu_setup\n rm %s.gnu_data",filename,filename);
 system(msg);
}

/***********************************************************************
  SINGULAR VALUE DECOMPOSITION

************************************************************************/
void
svd_window(item,value)
     Panel_item item;
     int value;
{
  SVD_COLS=value;
}

void
svd_on_off(item,value)
     Panel_item item;
     int value;
{
  if (value==1){
    EMBED='t';
    xv_set(SVD_M_SLIDER,PANEL_LABEL_STRING,"Time Delay",NULL);