band_pass_tv.c

最大似然估计算法

#include "timeseries.h"

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);

}