create_spectrum.c

最大似然估计算法

#include "timeseries.h"

int create_spectrum(double *P, double *f, time_series ts, data_kernel dk, options op) {  

	int j, k, l, idt;
	int nfft, len, jmax, nout;
	int n_segments;

	int *ind;

	double mle, prob;

	double *data, *fit, *cov_fit, *data_hat, *dy, *t;

	data_hat = (double *) calloc((size_t) dk.n_data,          sizeof(double));

	if (op.speed < 3) {
		dy       = (double *) calloc((size_t) dk.n_data,          sizeof(double));
		for (j = 0; j < dk.n_data; j++) dy[j] = 1.0;
		fit      = (double *) calloc((size_t) dk.n_par,           sizeof(double));
		cov_fit  = (double *) calloc((size_t)(dk.n_par*dk.n_par), sizeof(double));
        	linefit_fast(dk.A, dk.d, dy, dk.n_data, dk.n_par, fit, cov_fit, data_hat);
		free(fit);
		free(cov_fit);
		free(dy);
	} else {
		for (j = 0; j < dk.n_data; j++) data_hat[j] = 0.0;
	}

	/*********************************************************/
	/* Find if there is a segment in the time series that is */
	/* continuous and long enough to be used in the fft psd  */
	/*********************************************************/

	nfft = (int)pow(2.0,floor( log((double)ts.n_data) / log(2.0) ));

	ind = (int *) calloc((size_t) (ts.n_full-ts.n_data), sizeof(int));

	ind[0] = 0;
	for (k = 1, j = 1; j < ts.n_data; j++) {
		idt = ts.index[j]-ts.index[j-1];
		if (idt > 1) ind[k++] = j; 
	}
	ind[k++] = ts.n_data;
	n_segments = k;

	l = -1;
	for (j = 0; j < n_segments-1; j++) {
		len = ind[j+1]-ind[j];
		if (len >= nfft) { 
			l = j;
			fprintf(op.fpout, "%d to %d a continuous sequence (%.8f to %.8f)\n", ind[j], ind[j+1]-1, ts.t[ind[j]], ts.t[ind[j+1]-1]);
		}
	}

	/***********************************/
	/* Now estimate the power spectrum */
	/***********************************/

	if (l > 0) {
		data = (double *) calloc((size_t) nfft, sizeof(double));
		for (j = 0; j < nfft; j++) data[j] = dk.d[ind[l]+j] - data_hat[ind[l]+j];

		psd(P,f,data,nfft,ts.fs);
		for (j = 0; j < nfft/2; j++) P[j] *= (2.0 / ts.fs);
		free(data);
		free(data_hat);
		nout = nfft/2;
	} else {
		t = (double *) calloc((size_t) dk.n_data,     sizeof(double));

		for (j = 0; j < dk.n_data; j++) {
			t[j] = ts.t[j] - ts.t[0];
			t[j] *= sec_per_year;
			data_hat[j] = dk.d[j] - data_hat[j];
		}
		scargle(t,data_hat,dk.n_data,1.0,1.0,f,P,dk.n_data,&nout,&jmax,&prob);
		for (j = 0; j < nout; j++) P[j] *= (2.0 / ts.fs);
		free(t);
		free(data_hat);
	}
	free(ind);
	return(nout);
}