four1.c

最大似然估计算法

#include "timeseries.h"

#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;
	}

}