create_power_spectrum.c

最大似然估计算法

#include "timeseries.h"

void create_power_spectrum(noise_model nm, double *f, double fs, int n_psd,double *Power, int is_unit) {

int i, j, k;

double D, alpha, ss;
double a, b, fl, fh, pole, flfh2, fl2, fh2, f2;

ss =  is_unit ? 1.0  : nm.sigma[0] * nm.sigma[0];

switch (nm.model) {
	case 'p':
		alpha = nm.pvec[0];
		D = 2.0 * pow(2.0 * M_PI,alpha) * pow(sec_per_year,alpha/2.0);

		for (j = 0; j < n_psd; j++) Power[j] = D * ss * pow(f[j],alpha) / pow(fs,alpha/2.0 + 1.0);
		break;
	case 'f':
		/* Maybe use the discrete version here? */
		a = pow(nm.pvec[0] / sec_per_year,2.0);
		b = 4.0 * M_PI * M_PI;
		for (j = 0; j < n_psd; j++) Power[j] = 2.0 * ss / sec_per_year / (a + b * pow(f[j],2.0));
		break;
	case 'b':
		fl = nm.pvec[0] / nm.pvec[1];
                fh = nm.pvec[0] * nm.pvec[1];
		pole = nm.pvec[2];
		flfh2 = pow(fl+fh,2.0);
		fl2   = pow(fl,2.0);
		fh2   = pow(fh,2.0);
		for (j = 0; j < n_psd; j++) {
			f2 = f[j]*f[j];
			Power[j] = ss * pow(f2 * flfh2 / (fl2 + f2) / (fh2 + f2),pole);
		}
		break;
	case 'v':
		/* Note variable white noise is indistinguishable from white in a power spectrum */
	case 'w':
		for (j = 0; j < n_psd; j++) Power[j] = 2.0 * ss / fs;
		break;
	case 'g':
		/* We now have the equation for this */
		D = 2.0 * pow(sec_per_year,nm.pvec[1]/2.0) / pow(fs,nm.pvec[1]/2.0 + 1.0);
		a = pow(nm.pvec[0] / sec_per_year,2.0);
		b = 4.0 * M_PI * M_PI;
		for (j = 0; j < n_psd; j++) Power[j] = D * ss * pow(a + b * pow(f[j],2.0),nm.pvec[1]/2.0);
		break;
	case 's':
		fprintf(stderr, " This model does not have a spectral equivalent that\n");
		fprintf(stderr, " is different from white noise. Since it only makes sense\n");
		fprintf(stderr, " when accompanied by white noise then these two cant be\n");
		fprintf(stderr, " distinguished from each other in the spectral simplex.\n");
		exit(EXIT_FAILURE);
		break;
}

}