realft.c

最大似然估计算法

#include "timeseries.h"

void realft(double *data, int n, int isign) {


	int i, i1, i2, i3, i4, n2p3;

	double c1, c2, h1i, h1r, h2i, h2r;
	double theta, wi, wpi, wpr, wr, wtemp;

	c1 = 0.5;
	theta = M_PI / ((double)n / 2.0);

	if (isign == 1) {
		c2 = -0.5;
		four1(data,n/2,1);
	} else {
		c2 = 0.5;
		theta = -theta;
	}
	wpr = -2.0 * pow(sin(0.5*theta),2.0);
	wpi = sin(theta);
	wr = 1.0 + wpr;
	wi = wpi;
	n2p3 = n + 3;
	for (i = 2; i <= n/4; i++) {
		i1 = 2*i-1;
		i2 = i1+1;
		i3 = n2p3 - i2;
		i4 = i3 + 1;

		h1r =  c1 * (data[i1-1]+data[i3-1]);
		h1i =  c1 * (data[i2-1]-data[i4-1]);
		h2r = -c2 * (data[i2-1]+data[i4-1]);
		h2i =  c2 * (data[i1-1]-data[i3-1]);

		data[i1-1] =  h1r + wr * h2r - wi * h2i;
		data[i2-1] =  h1i + wr * h2i + wi * h2r;
		data[i3-1] =  h1r - wr * h2r + wi * h2i;
		data[i4-1] = -h1i + wr * h2i + wi * h2r;
		wtemp = wr;
		wr = wr * wpr - wi * wpi + wr;
		wi = wi * wpr + wtemp * wpi + wi;
	}
	if (isign == 1) {
		h1r = data[0];
		data[0] = h1r + data[1];
		data[1] = h1r - data[1];
	} else {
		h1r = data[0];
		data[0] = c1 * (h1r + data[1]);
		data[1] = c1 * (h1r - data[1]);
		four1(data,n/2,-1);
	}
}