IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / 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