binaries/src/disembl/Tisean_3.0.1/source_c/mem_spec.c

   1 /*
   2  *   This file is part of TISEAN
   3  *
   4  *   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   5  *
   6  *   TISEAN is free software; you can redistribute it and/or modify
   7  *   it under the terms of the GNU General Public License as published by
   8  *   the Free Software Foundation; either version 2 of the License, or
   9  *   (at your option) any later version.
  10  *
  11  *   TISEAN is distributed in the hope that it will be useful,
  12  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  *   GNU General Public License for more details.
  15  *
  16  *   You should have received a copy of the GNU General Public License
  17  *   along with TISEAN; if not, write to the Free Software
  18  *   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  19  */
  20 /*Author: Rainer Hegger */
  21 /*Changes:
  22   Feb 19, 2007: changed meaning of -f flag and added -P flag to be
  23                 consistent with spectrum
  24   Dec 5, 2006: Seg fault when poles > length;
  25  */
  26 #include <stdio.h>
  27 #include <string.h>
  28 #include <stdlib.h>
  29 #include <limits.h>
  30 #include "routines/tsa.h"
  31 #include <math.h>
  32
  33 #define WID_STR "Estimates the power spectrum of the data"
  34
  35 #ifndef M_PI
  36 #define M_PI 3.1415926535897932385E0
  37 #endif
  38
  39 unsigned long poles=128,out=2000;
  40 unsigned long length=ULONG_MAX,exclude=0;
  41 unsigned int column=1;
  42 unsigned int verbosity=0x1;
  43 double samplingrate=1.0;
  44 char *outfile=NULL,stdo=1;
  45 char *infile=NULL;
  46 double *series;
  47
  48 void show_options(char *progname)
  49 {
  50   what_i_do(progname,WID_STR);
  51   fprintf(stderr," Usage: %s [Options]\n",progname);
  52   fprintf(stderr," Options:\n");
  53   fprintf(stderr,"Everything not being a valid option will be interpreted"
  54           " as a possible"
  55           " datafile.\nIf no datafile is given stdin is read. Just - also"
  56           " means stdin\n");
  57   fprintf(stderr,"\t-l length of file [default is whole file]\n");
  58   fprintf(stderr,"\t-x # of lines to be ignored [default is 0]\n");
  59   fprintf(stderr,"\t-c column to read [default is 1]\n");
  60   fprintf(stderr,"\t-p number of poles [default is 128 or file length]\n");
  61   fprintf(stderr,"\t-P number of frequences out [default is 2000]\n");
  62   fprintf(stderr,"\t-f sampling rate in Hz [default is 1]\n");
  63   fprintf(stderr,"\t-o outfile [default is 'datafile'.spec]\n");
  64   fprintf(stderr,"\t-V verbosity level [default is 1]\n\t\t"
  65           "0='only panic messages'\n\t\t"
  66           "1='+ input/output messages'\n\t\t"
  67           "2='+ print the ar coefficients too'\n");
  68   fprintf(stderr,"\t-h show these options\n\n");
  69   exit(0);
  70 }
  71
  72 void scan_options(int argc,char **argv)
  73 {
  74   char *hout;
  75
  76   if ((hout=check_option(argv,argc,'l','u')) != NULL)
  77     sscanf(hout,"%lu",&length);
  78   if ((hout=check_option(argv,argc,'x','u')) != NULL)
  79     sscanf(hout,"%lu",&exclude);
  80   if ((hout=check_option(argv,argc,'c','u')) != NULL)
  81     sscanf(hout,"%u",&column);
  82   if ((hout=check_option(argv,argc,'p','u')) != NULL)
  83     sscanf(hout,"%lu",&poles);
  84   if ((hout=check_option(argv,argc,'P','u')) != NULL)
  85     sscanf(hout,"%lu",&out);
  86   if ((hout=check_option(argv,argc,'f','f')) != NULL)
  87     sscanf(hout,"%lf",&samplingrate);
  88   if ((hout=check_option(argv,argc,'V','u')) != NULL)
  89     sscanf(hout,"%u",&verbosity);
  90   if ((hout=check_option(argv,argc,'o','o')) != NULL) {
  91     stdo=0;
  92     if (strlen(hout) > 0)
  93       outfile=hout;
  94   }
  95 }
  96
  97 double getcoefs(double *coef)
  98 {
  99   long i,j,hp=(long)poles-1;
 100   double ret=0.0,*cov,*help,h1,h2;
 101
 102   check_alloc(cov=(double*)malloc(sizeof(double)*length));
 103   check_alloc(help=(double*)malloc(sizeof(double)*poles));
 104
 105   for (i=0;i<length;i++)
 106     ret += series[i]*series[i];
 107   ret /= length;
 108
 109   for (i=0;i<length;i++)
 110     cov[i]=series[i];
 111   series++;
 112
 113   for (i=0;i<poles;i++) {
 114     h1=h2=0.0;
 115     for (j=0;j<length-i-1;j++) {
 116       h1 += cov[j]*series[j];
 117       h2 += cov[j]*cov[j]+series[j]*series[j];
 118     }
 119     coef[i]=2.0*h1/h2;
 120     ret *= (1.0-coef[i]*coef[i]);
 121     for (j=0;j<i;j++)
 122       coef[j]=help[j]-coef[i]*help[i-1-j];
 123     if (i == hp)
 124       break;
 125     for (j=0;j<=i;j++)
 126       help[j]=coef[j];
 127     for (j=0;j<length-i-1;j++) {
 128       cov[j] -= help[i]*series[j];
 129       series[j]=series[j+1]-help[i]*cov[j+1];
 130     }
 131   }
 132   free(cov);
 133   free(help);
 134
 135   return ret;
 136 }
 137 double powcoef(double dt,double *coef)
 138 {
 139   int i;
 140   double si=0.0,sr=1.0,zr=1.0,zi=0.0,h,omdt,hr,hi;
 141
 142   omdt=2.0*M_PI*dt;
 143   hr=cos(omdt);
 144   hi=sin(omdt);
 145
 146   for (i=0;i<poles;i++) {
 147     h=zr;
 148     zr=zr*hr-zi*hi;
 149     zi=h*hi+zi*hr;
 150     sr -= coef[i]*zr;
 151     si -= coef[i]*zi;
 152   }
 153   return (sr*sr+si*si);
 154 }
 155
 156 int main(int argc,char **argv)
 157 {
 158   char stdi=0;
 159   double fdt,pm,pow_spec,*cof,av,var;
 160   long i;
 161   FILE *fout;
 162
 163   if (scan_help(argc,argv))
 164     show_options(argv[0]);
 165
 166   scan_options(argc,argv);
 167 #ifndef OMIT_WHAT_I_DO
 168   if (verbosity&VER_INPUT)
 169     what_i_do(argv[0],WID_STR);
 170 #endif
 171
 172   infile=search_datafile(argc,argv,&column,verbosity);
 173   if (infile == NULL)
 174     stdi=1;
 175
 176   if (outfile == NULL) {
 177     if (!stdi) {
 178       check_alloc(outfile=(char*)calloc(strlen(infile)+6,(size_t)1));
 179       strcpy(outfile,infile);
 180       strcat(outfile,".spec");
 181     }
 182     else {
 183       check_alloc(outfile=(char*)calloc((size_t)11,(size_t)1));
 184       strcpy(outfile,"stdin.spec");
 185     }
 186   }
 187   if (!stdo)
 188     test_outfile(outfile);
 189
 190   series=(double*)get_series(infile,&length,exclude,column,verbosity);
 191
 192   if (length <= poles) {
 193     fprintf(stderr,"\n\tNo. of poles has to be smaller then the length of the\n"
 194             "\tdata set! Exiting.\n");
 195     exit(MEM_SPEC_TOO_MANY_POLES);
 196   }
 197
 198   variance(series,length,&av,&var);
 199   for (i=0;i<length;i++)
 200     series[i] -= av;
 201
 202   check_alloc(cof=(double*)malloc(sizeof(double)*poles));
 203
 204   pm=getcoefs(cof);
 205
 206   if (!stdo) {
 207     fout=fopen(outfile,"w");
 208     if (verbosity&VER_INPUT)
 209       fprintf(stderr,"Opened %s for writing\n",outfile);
 210     if (verbosity&VER_USR1) {
 211       fprintf(fout,"#sigma^2=%e\n",pm);
 212       for (i=0;i<poles;i++)
 213         fprintf(fout,"#%ld %e\n",i+1,cof[i]);
 214     }
 215     for(i=0;i<out;i++) {
 216       fdt=i/(2.0*out);
 217       pow_spec=powcoef(fdt,cof);
 218       fprintf(fout,"%e %e\n",fdt*samplingrate,
 219               pm/pow_spec/sqrt((double)length));
 220       fflush(fout);
 221     }
 222     fclose(fout);
 223   }
 224   else {
 225     if (verbosity&VER_INPUT)
 226       fprintf(stderr,"Writing to stdout\n");
 227     if (verbosity&VER_USR1) {
 228       fprintf(stdout,"#sigma^2=%e\n",pm);
 229       for (i=0;i<poles;i++)
 230         fprintf(stdout,"#%ld %e\n",i+1,cof[i]);
 231     }
 232     for(i=0;i<out;i++) {
 233       fdt=i/(2.0*out);
 234       pow_spec=powcoef(fdt,cof);
 235       fprintf(stdout,"%e %e\n",fdt*samplingrate,
 236               pm/pow_spec/*/sqrt((double)length)*/);
 237     }
 238   }
 239
 240   return 0;
 241 }
 242