#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <complex.h>
#include <math.h>
#include <fftw3.h>

/*
 * fourierfit.c:
 * Fits a series of sines/cosines onto a series of character values from a
 * string.
 *
 * Copyright (c) 2008, Jan Krueger <jast@heapsort.de>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * Compile like this:
 * cc -o fourierfit -lm -lfftw3 fourierfit.c
 */

/*
 * Precision of output function (in digits after decimal point).
 * Yes, it's a string. printf much?
 */
#define PRECISION "3"
/* Names of cos/sin functions to use in output */
#define COSNAME "cos"
#define SINNAME "sin"
/* Name of function variable */
#define FNAME "t"

/* Don't change below here etc. */

#define PREC_FORMAT "%."PRECISION"f"
const char *cos_format = " %c "PREC_FORMAT"*"COSNAME"("PREC_FORMAT"*"FNAME
	")";
const char *sin_format = " %c "PREC_FORMAT"*"SINNAME"("PREC_FORMAT"*"FNAME
	")";

void format_term(const char *f, double d, int p, int N)
{
	printf(f, (d >= 0 ? '+' : '-'),
		fabs(d), /* scale of function */
		(double)p * M_PI * 2.0 / (double)N /* scale of parameter */
	);
}

int main(int argc, char **argv)
{
	int N;
	int odd;
	int i;
	double *in;
	fftw_complex *out;
	fftw_plan p;
	if (argc != 2) {
		puts("Syntax: fourierfit <string>");
		exit(1);
	}
	N = strlen(argv[1]);
	odd = N % 2;

	in = fftw_malloc(sizeof (double) * N);
	out = fftw_malloc(sizeof (fftw_complex) * N);

	p = fftw_plan_dft_r2c_1d(N, in, out,
		FFTW_ESTIMATE | FFTW_DESTROY_INPUT);

	/* Fill in array */
	for (i = 0; i < N; i++) {
		in[i] = (double) argv[1][i];
	}

	fftw_execute(p);
	fftw_destroy_plan(p);
	fftw_free(in);

	/* "DC" */
	printf(PREC_FORMAT, creal(out[0])/N);

	for (i = 1; i <= (N-1)/2; i++) {
		format_term(cos_format, 2.0*creal(out[i])/N, i, N);
		format_term(sin_format, -2.0*cimag(out[i])/N, i, N);
	}
	/* Nyquist component */
	if (!odd) {
		format_term(cos_format, creal(out[N/2])/N, N/2, N);
	}
	printf("\n");
	fftw_free(out);
	exit(0);
}