psd.c

最大似然估计算法

#include "timeseries.h"

void psd(double *P, double *f, double *d, int nfft, double fs) {

	int j, k, i, lwork, info;

	double KMU;

	double *array, *window;
	double *work, *U, *S, *VT;

	array    = (double *) calloc((size_t)(nfft*2),     sizeof(double));
	window   = (double *) calloc((size_t)(nfft),       sizeof(double));
	for (j = 0; j < nfft; j++) {
		array[j*2]  = d[j] * 0.5 * (1.0 - cos(2.0 * M_PI * (double) (j+1) / (double)(nfft+1) ) );
		window[j] = 0.5 * (1.0 - cos(2.0 * M_PI * (double) (j+1) / (double)(nfft+1) ) );
	}

	i = 1;
	lwork = 3 + nfft;
	work = (double *) calloc((size_t)lwork, sizeof(double));
	U    = (double *) calloc((size_t)i,     sizeof(double));
	S    = (double *) calloc((size_t)i,     sizeof(double));
	VT   = (double *) calloc((size_t)i,     sizeof(double));

	dgesvd_("N","N",&nfft,&i,window,&nfft,S,U,&i,VT,&i,work,&lwork,&info);

	KMU = S[0] * S[0];

	fft1d(array,nfft,1);

	for (k = 0, j = 2; k < nfft/2; j+= 2, k++) {
		f[k] = (double) (k+1) * fs / (double) nfft;
		P[k] = (array[j]*array[j] + array[j+1]*array[j+1]) / KMU;
	}
	free(array);
	free(window);
	free(work);
	free(U);
	free(S);
	free(VT);

}