transform_vectors.c

最大似然估计算法

#include "timeseries.h"

void realft(double *, int , int);
void four1(double *, int, int);

void power_law_tv(double *tv, int n, double degree) {

	int i, j;

	double *b;

	b = (double *)calloc((size_t)n, sizeof(double));

        b[0] = 1.0;

        for (j = 0; j < n-1; j++)  b[j+1] = ((double) j - degree) / ((double) (j+1));

        tv[0] = b[0];
        for (j = 1; j < n; j++) tv[j] = tv[j-1] * b[j];

	free(b);
}

void fogm_tv(double *tv, int n, double beta, double dt) {

	int j;

	double sec_per_year, fs, Beta, a;

        sec_per_year = 365.24219*24.0*3600.0;

	Beta = beta / sec_per_year;

	fs = 1.0 / dt / sec_per_year;

	a = exp(-Beta / fs) * sqrt(dt) + exp(-fs/Beta) /  Beta * sqrt(fs/sec_per_year);

	for (j = 0; j < n; j++) {
		tv[j] = a * exp(-Beta * sec_per_year * (double) j * dt);
	}

}


void band_pass_tv(double *tv, int n, double flow, double fhigh, int npole) {

/*****************************************************************/
/*    band_pass_cov                                              */
/*                                                               */
/*    Compute the band-pass transfromation vector                */
/*    adapted from John Langbein's code                          */
/*                                                               */
/*    sigma          standard deviation of noise (meters)        */
/*    flow                                                       */
/*    fhigh                                                      */
/*    npole                                                      */
/*                                                               */
/* Simon  Williams, [e-mail], 25-Jan-2002 c-version       */
/*****************************************************************/

	int i, j, k, kk;
	int len, nf, imax, kmax;

	double fll, fhh, fo;
	double h1r, h1i, h2r, h2i, h3r, h3i, h4r, h4i;
	double h1h2r, h1h2i, h3h4r, h3h4i;
	double hcr, hci, re, im, re2, im2;
	double calc_amp, a, f, dsum, temp;

	double *H;

	if (fhigh < flow) {  
        	temp  = fhigh;
        	fhigh = flow;
        	flow  = temp;
	}

	fll = flow  / 365.24219;
	fhh = fhigh / 365.24219;
	
	fo = (fll + fhh) / 2.0;

	a = log( (double)n) / log(2.0);
	a = ceil(a);

	len = (int) pow(2.0, a);
	nf = len / 2;

	/* fprintf(stdout, "In band_pass n = %d, len = %d, nf = %d\n", n, len, nf); */
	H  = (double *) calloc((size_t) len, sizeof(double));

	H[0] = 1.0;
	for (i = 1; i < nf; i++) {
		f = (1.0 / (double)len) * (double)i;
		h1r = 1.0;
		h1i = 0.0;
		h2r = 0.0;
		h2i = f / fll;
		h3r = 1.0;
		h3i = f / fll;
		h4r = 1.0;
		h4i = f / fhh;

		h1h2r = (h1r * h2r - h1i * h2i);
		h1h2i = (h1r * h2i + h1i * h2r);

		h3h4r = (h3r * h4r - h3i * h4i);
		h3h4i = (h3r * h4i + h3i * h4r);

		hcr = (h1h2r * h3h4r + h1h2i * h3h4i) / (h3h4r*h3h4r + h3h4i * h3h4i);
		hci = (h3h4r * h1h2i - h1h2r * h3h4i) / (h3h4r*h3h4r + h3h4i * h3h4i);
		re = hcr;
		im = hci;

		if (npole > 1) {
			re2 = hcr*re - hci*im;
			im2 = hci*re + hcr*im;
			re = re2;
			im = im2;
		} 
		if (npole > 2) {
			re2 = hcr*re - hci*im;
			im2 = hci*re + hcr*im;
			re = re2;
			im = im2;
		}
		if (npole > 3) {
			re2 = hcr*re - hci*im;
			im2 = hci*re + hcr*im;
			re = re2;
			im = im2;
		}

		H[i*2]   = re;
		H[i*2+1] = -im;
			
	}
	H[1] = sqrt(re*re + im*im);

	calc_amp = (fo/fll) / (sqrt(1.0 + (fo/fll)*(fo/fll)) * sqrt(1.0 + (fo/fhh)*(fo/fhh)));
	calc_amp = pow(calc_amp,(double)npole);

	realft(H,len,-1);
	
	for (i = 0; i < n; i++) {
		tv[i] = (1.0 / calc_amp) * ( H[i] * 2.0 / (double)len) * sqrt(182.621095);
	}

	free(H);

}

void realft(double *data, int n, int isign) {


	int i, i1, i2, i3, i4, n2p3;

	double c1, c2, h1i, h1r, h2i, h2r;
	double theta, wi, wpi, wpr, wr, wtemp;

	c1 = 0.5;
	theta = M_PI / ((double)n / 2.0);

	if (isign == 1) {
		c2 = -0.5;
		four1(data,n/2,1);
	} else {
		c2 = 0.5;
		theta = -theta;
	}
	wpr = -2.0 * pow(sin(0.5*theta),2.0);
	wpi = sin(theta);
	wr = 1.0 + wpr;
	wi = wpi;
	n2p3 = n + 3;
	for (i = 2; i <= n/4; i++) {
		i1 = 2*i-1;
		i2 = i1+1;
		i3 = n2p3 - i2;
		i4 = i3 + 1;

		h1r =  c1 * (data[i1-1]+data[i3-1]);
		h1i =  c1 * (data[i2-1]-data[i4-1]);
		h2r = -c2 * (data[i2-1]+data[i4-1]);
		h2i =  c2 * (data[i1-1]-data[i3-1]);

		data[i1-1] =  h1r + wr * h2r - wi * h2i;
		data[i2-1] =  h1i + wr * h2i + wi * h2r;
		data[i3-1] =  h1r - wr * h2r + wi * h2i;
		data[i4-1] = -h1i + wr * h2i + wi * h2r;
		wtemp = wr;
		wr = wr * wpr - wi * wpi + wr;
		wi = wi * wpr + wtemp * wpi + wi;
	}
	if (isign == 1) {
		h1r = data[0];
		data[0] = h1r + data[1];
		data[1] = h1r - data[1];
	} else {
		h1r = data[0];
		data[0] = c1 * (h1r + data[1]);
		data[1] = c1 * (h1r - data[1]);
		four1(data,n/2,-1);
	}
}

#define SWAP(a,b) tempr = (a); (a) = (b); (b) = tempr

void four1(double *data, int nn, int isign) {

	int n, mmax, m, j, istep, i;
	double wtemp, wr, wpr, wpi, wi, theta;
	double tempr, tempi;

	n = 2 * nn;

	j = 1;

	for (i = 1; i < n; i+= 2) {
		if (j > i) {
			SWAP(data[j-1],data[i-1]);
			SWAP(data[j],data[i]);
		}
		m = n / 2;
		while ( m >= 2 && j > m) {
			j -= m;
			m /= 2;
		}
		j += m;
	}

	mmax = 2;
	while ( n > mmax) {
		istep = 2 * mmax;
		theta = (2.0 * M_PI) / (double)(isign * mmax);
		wtemp = sin(0.5*theta);
		wpr = -2.0 * wtemp * wtemp;
		wpi = sin(theta);
		wr = 1.0;
		wi = 0.0;
		for(m = 1; m < mmax; m+=2) {
			for (i = m; i <= n; i+= istep) {
				j=i+mmax;
				tempr = wr * data[j-1]-wi*data[j];
				tempi = wr * data[j] + wi*data[j-1]; 
				data[j-1]  = data[i-1] - tempr;
				data[j]    = data[i]   - tempi;
				data[i-1] += tempr;
				data[i]   += tempi;
			}
			wr = (wtemp=wr)*wpr-wi*wpi+wr;
			wi = wi*wpr+wtemp*wpi+wi;
		}
		mmax = istep;
	}

}