sing.c

时间序列工具

/* certain changes (mainly global names) by Andreas Schmitz (1997) */
/* certain changes (mainly cleanup) by Thomas Schreiber (1998)     */

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

	Javier Soley, Ph. D,   FJSOLEY @UCRVM2.BITNET
	Escuela de F韘ica y Centro de Investigaciones Geof韘icas
	Universidad de Costa Rica

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

/*	Computes the DISCRETE FOURIER TRANSFORM of very long data series.
 *
 *   This functions are translations from the fortran program in
 *
 *   R. C. Singleton, An algorithm for computing the mixed radix fast
 *   Fourier transform 
 *
 *	IEEE Trans. Audio Electroacoust., vol. AU-17, pp. 93-10, June 1969.
 *   Some features are:
 *
 *		1-) Accepts an order of transform that can be factored not only
 *		    in prime factors such 2 and 4, but also including odd factors
 *		    as 3, 5, 7, 11, etc.
 *		2-) Generates sines and cosines recursively and includes
 *		    corrections for truncation errors.
 *		3-) The original subroutine accepts multivariate data. This 
 *			translation does not implement that option (because I
 *			do not needed right now).
 *
 *	Singleton wrote his subroutine in Fortran and in such a way that it 
 *	could be ported allmost directly to assembly language. I transcribed
 *   it to C with little effort to make it structured. So I apologize to
 *   all those C purists out there!!!!!!!!
 *
 */

				/*  Version 2.0 March/30/92 */
/* Includes */

#include <stdlib.h>
#include <stddef.h>
#include <stdio.h>
#include <math.h>

/* Defines */
#ifndef M_PI
   #define M_PI   3.14159265358979323846       /* should have been in math.h */
#endif
#ifndef M_PI_2
   #define M_PI_2 1.57079632679489661923       /* should have been in math.h */
#endif
#define	  TWO_PI	((double)2.0 * M_PI)
#define   MAXF  20000
#define   MAXP  20000 

/* Globals */

long nn_gl,m_gl,flag_gl,jf_gl,jc_gl,kspan_gl,ks_gl, kt_gl,nt_gl,kk_gl,i_gl;
double c72_gl, s72_gl, s120_gl, cd_gl, sd_gl, rad_gl, radf_gl;
double at_gl[MAXF], bt_gl[MAXF];
long nfac_gl[MAXF];   
int inc_gl;
long np_gl[MAXP];


/* The functions */

void radix_2(double *a, double *b)
{ 
  long k1, k2;
  double ak, bk, c1, s1;
  kspan_gl >>= 1;
  k1 = kspan_gl +2;
  
  do {
    do  {
      k2 = kk_gl + kspan_gl;
      ak = a[k2-1];
      bk = b[k2-1];
      a[k2-1] = a[kk_gl-1] -ak;
      b[k2-1] = b[kk_gl-1] -bk;
      a[kk_gl-1] += ak;
      b[kk_gl-1] += bk;
      kk_gl = k2 + kspan_gl;
    } while ( kk_gl <= nn_gl);
    kk_gl = kk_gl - nn_gl;
  } while ( kk_gl <= jc_gl);
  
  if ( kk_gl > kspan_gl) 
    flag_gl = 1; 
  else
    {
      do {
	c1 = 1.0 - cd_gl;
	s1 = sd_gl;
	do {
	  do  {
	    do {
	      k2 = kk_gl + kspan_gl;
	      ak = a[kk_gl-1]- a[k2-1];
	      bk = b[kk_gl-1]- b[k2-1];
	      a[kk_gl-1] += a[k2-1];
	      b[kk_gl-1] += b[k2-1];
	      a[k2-1]  = c1*ak - s1*bk;
	      b[k2-1]  = s1*ak + c1*bk;
	      kk_gl = k2 + kspan_gl;
	    } while ( kk_gl < nt_gl );
	    k2 = kk_gl - nt_gl;
	    c1 = -c1;
	    kk_gl = k1 - k2;
	  } while ( kk_gl > k2 );
	  ak = c1- (cd_gl*c1+sd_gl*s1);
	  s1 = (sd_gl*c1-cd_gl*s1) +s1;
	  
	  /***** Compensate for truncation errors   *****/
	  
	  c1 = 0.5/(ak*ak+s1*s1)+0.5;
	  s1 *= c1;
	  c1 *= ak;
	  kk_gl += jc_gl;
	} while ( kk_gl < k2);
	k1 = k1 + inc_gl + inc_gl;
	kk_gl = (k1- kspan_gl) /2 + jc_gl;
      } while ( kk_gl <= jc_gl + jc_gl );
    }
}

void radix_4(int isn, double *a, double *b)
{
  long	k1, k2, k3;
  double  akp, akm, ajm, ajp, bkm, bkp, bjm, bjp;
  double  c1, s1, c2, s2, c3, s3;

  kspan_gl /= 4;
 cuatro_1: 
  c1 = 1.0;
  s1 = 0;
  do {
    do {
      do {
	k1  =  kk_gl + kspan_gl;
	k2  =  k1 + kspan_gl;
	k3  =  k2 + kspan_gl;
	akp =  a[kk_gl-1] + a[k2-1];
	akm =  a[kk_gl-1] - a[k2-1];
	ajp =  a[ k1-1] + a[k3-1];
	ajm =  a[ k1-1] - a[k3-1];
	a[kk_gl-1] = akp + ajp;
	ajp = akp - ajp;
	bkp = b[kk_gl-1] + b[k2-1];
	bkm = b[kk_gl-1] - b[k2-1];
	bjp = b[k1-1] + b[k3-1];
	bjm = b[k1-1] - b[k3-1];
	b[kk_gl-1] = bkp + bjp;
	bjp = bkp - bjp;
	if ( isn < 0) goto cuatro_5;
	akp = akm - bjm;
	akm = akm + bjm;
	bkp = bkm + ajm;
	bkm = bkm - ajm;
	if (s1 == 0.0) goto cuatro_6;
      cuatro_3:		a[ k1-1] = akp*c1 - bkp*s1;
	b[ k1-1] = akp*s1 + bkp*c1;
	a[ k2-1] = ajp*c2 - bjp*s2;
	b[ k2-1] = ajp*s2 + bjp*c2;
	a[ k3-1] = akm*c3 - bkm*s3;
	b[ k3-1] = akm*s3 + bkm*c3;
	kk_gl = k3 + kspan_gl;
      }  while ( kk_gl <= nt_gl);
      
    cuatro_4: 
      c2 = c1 - (cd_gl*c1 + sd_gl*s1);
      s1 = (sd_gl*c1 - cd_gl*s1) + s1;
      
      /***** Compensate for truncation errors *****/
      
      c1 = 0.5 / (c2*c2 + s1*s1) +0.5;
      s1 = c1 * s1;
      c1 = c1 * c2;
      c2 = c1*c1 - s1*s1;
      s2 = 2.0 * c1 *s1;
      c3 = c2*c1 - s2*s1;
      s3 = c2*s1 + s2*c1;
      kk_gl = kk_gl -nt_gl + jc_gl;
    } while ( kk_gl <= kspan_gl);
    
    kk_gl = kk_gl - kspan_gl + inc_gl;
    if ( kk_gl <= jc_gl)  goto cuatro_1;
    if ( kspan_gl == jc_gl)  flag_gl =1;
    s1 = s1 + 0.0;
    return;
  cuatro_5:
    akp = akm + bjm;
    akm = akm - bjm;
    bkp = bkm - ajm;
    bkm = bkm + ajm;
    if (s1 != 0.0) goto cuatro_3;
  cuatro_6:
    a[k1-1] = akp;
    b[k1-1] = bkp;
    b[k2-1] = bjp;
    a[k2-1] = ajp;
    a[k3-1] = akm;
    b[k3-1] = bkm;
    kk_gl = k3 + kspan_gl;
  } while ( kk_gl <= nt_gl);
  goto cuatro_4;
}

/* Find prime factors of n */

void fac_des(long  n)
{
  long k, j, jj;
  k  = n;
  m_gl = 0;
  while ( k-(k / 16)*16 == 0 ) {
    m_gl++;
    nfac_gl[m_gl-1] = 4;
    k /= 16;
  } 
  j  = 3;
  jj = 9;
  do {
    while (k % jj == 0) {
      m_gl++;
      nfac_gl[m_gl-1] = j;
      k /= jj;
    }
    j += 2;
    jj = j * j;
  } while ( jj <= k);
  if (k <= 4) {
    kt_gl = m_gl;
    nfac_gl[m_gl] = k;
    if (k != 1) m_gl++;
  }
  else {
    if (k-(k / 4)*4 == 0) {
      m_gl++;
      nfac_gl[m_gl-1] = 2;
      k /= 4;
    }
    kt_gl = m_gl;
    j = 2;
    do {
      if (k % j == 0 ) {
	m_gl++;
	nfac_gl[m_gl-1] = j;
	k /= j;
      }
      j = ((j+1)/ 2)*2 + 1;
    } while ( j <= k);
  }
  if (kt_gl != 0) {
    j = kt_gl;
    do {
      m_gl++;
      nfac_gl[m_gl-1] = nfac_gl[j-1];
      j--;
    } while ( j != 0);
  }
}    

/* Permute the results to normal order  */

void permute(long n, double *a, double *b)
{
  long  k, j, k1, k2, k3, kspnn, maxf;
  
  double ak, bk;
  long  ii, jj;
  maxf = MAXF;
  np_gl[0] = ks_gl;
  if (kt_gl != 0) {
    k = kt_gl +kt_gl +1;
    if (m_gl < k) k--;
    j = 1;
    np_gl[k] = jc_gl;
    do {
      np_gl[j]   = np_gl[j-1] / nfac_gl[j-1];
      np_gl[k-1] = np_gl[k]   * nfac_gl[j-1];
      j++;
      k--;
    } while (j < k);
    k3 = np_gl[k];   
    kspan_gl = np_gl[1];
    kk_gl = jc_gl+1;
    k2 = kspan_gl + 1;  
    j = 1;
    
    do {
    	do {
    	  ak      = a[kk_gl-1];
    	  a[kk_gl-1] = a[k2-1];
    	  a[k2-1] = ak;
    	  bk      = b[kk_gl-1];
    	  b[kk_gl-1] = b[k2-1];
    	  b[k2-1] = bk;
    	  kk_gl     += inc_gl;
    	  k2     += kspan_gl;
    	} while ( k2 < ks_gl);
    ocho_30:  		do {
    	k2 -= np_gl[j-1];
    	j++;
    	k2 += np_gl[j];
    } while (k2 > np_gl[j-1]);
    	j = 1;
    ocho_40:  		j = j + 0;
    } while (kk_gl < k2);
    kk_gl += inc_gl;
    k2 += kspan_gl;
    if (k2 < ks_gl)  goto ocho_40;
    if (kk_gl < ks_gl)  goto ocho_30;
    jc_gl = k3;
  }
  
  if ( (2*kt_gl +1) < m_gl) {
    kspnn = np_gl[kt_gl];
    /* Permutation of square-free factors of n */
    j = m_gl - kt_gl;
    nfac_gl[j] = 1;
    do {
      nfac_gl[j-1] *= nfac_gl[j];
      j--;
    } while (j != kt_gl);
    kt_gl++;
    nn_gl = nfac_gl[kt_gl-1] -1;
    if (nn_gl > MAXP) {
      printf("product of square free factors exceeds allowed limit\n");
      exit(2);
    }
    jj =0;
    j=0;
    goto nueve_06;
  nueve_02:
    jj -= k2;
    k2 = kk_gl;
    k++;
    kk_gl = nfac_gl[k-1];
    do {
      jj += kk_gl;
      if ( jj >= k2 ) goto nueve_02;
      np_gl[j-1] = jj;
    nueve_06:
      k2 = nfac_gl[kt_gl-1];
      k = kt_gl+1;
      kk_gl = nfac_gl[k-1];
      j++;
    } while (j <= nn_gl);
    /* determine the permutation cycles  of length greater then 1 */
    j =0;
    goto nueve_14;
    do {
      do {		
	k = kk_gl;
	kk_gl = np_gl[k-1];
	np_gl[k-1] = -kk_gl;
      } while ( kk_gl != j);
      k3 = kk_gl;
    nueve_14:
      do {
	j++;
	kk_gl = np_gl[j-1];
      } while (kk_gl <0);
    } while ( kk_gl != j);
    np_gl[j-1] = -j;
    if (j != nn_gl) goto nueve_14;
    maxf *= inc_gl;				
    /* Reorder a and b following the permutation cycles */
    goto nueve_50;
    do {
      do {
	do { j--;} while (np_gl[j-1] <0);
	jj = jc_gl;
	do {