// fftfun.C

/******************************************************************************
 *
 *  MiXViews - an X window system based sound & data editor/processor
 *
 *  Copyright (c) 1993, 1994 Regents of the University of California
 *
 *  Author:     Douglas Scott
 *  Date:       December 13, 1994
 *
 *  Permission to use, copy and modify this software and its documentation
 *  for research and/or educational purposes and without fee is hereby granted,
 *  provided that the above copyright notice appear in all copies and that
 *  both that copyright notice and this permission notice appear in
 *  supporting documentation. The author reserves the right to distribute this
 *  software and its documentation.  The University of California and the author
 *  make no representations about the suitability of this software for any 
 *  purpose, and in no event shall University of California be liable for any
 *  damage, loss of data, or profits resulting from its use.
 *  It is provided "as is" without express or implied warranty.
 *
 ******************************************************************************/


#ifdef __GNUG__
#pragma implementation
#endif

#include <math.h>
#include "fftfun.h"
#include "fftdata.h"
#include "request.h"
#include "requester.h"
#include "windowtable.h"

class FFTRequester : public TitledRequester {
public:
	FFTRequester(FFT_Function *);
protected:
	redefined void configureRequest(Request *);
	redefined boolean confirmValues();
private:
	FFT_Function* client;
	int foffset;
	double frate;
};

FFTRequester::FFTRequester(FFT_Function* f)
	: TitledRequester("Fast-Fourier Transform Analysis of Selected Region:"),
	  client(f), foffset(100), frate(0.0) {}

void
FFTRequester::configureRequest(Request* request) {
	request->appendValue("Frame offset (samples):",
			    &foffset, Range(1, client->target()->length() / 2),
			    false);
	request->appendValue("Frame rate (frames/second):", &frate,
	                     PositiveNumbers);
	request->appendChoice("Number of FFT points:",
	                      "|64|128|256|512|1024|2048|",
	                      &client->npoints);
}

boolean
FFTRequester::confirmValues() {
	// { 1, 2, 4, .. } => { 64, 128, 256, ... }
	client->npoints <<= 6;
	return client->setBaseValues(client->npoints * 2, frate, foffset);
}

//********

FFT_Function::FFT_Function(Data* data, int points, int frameoffset)
		: ArrayFunction(data, points*2, frameoffset),
		  xreal(nil), ximag(nil), window(nil), npoints(points) {
	initialize();
}

// FIX ME:  put default values into ArrayFunction args

FFT_Function::FFT_Function(Data* data) : ArrayFunction(data),
	xreal(nil), ximag(nil), window(nil), npoints(P_512) {}

FFT_Function::~FFT_Function() {
	delete [] xreal;
	delete [] ximag;
	Resource::unref(window);
}

Requester *
FFT_Function::createRequester() {
	return new FFTRequester(this);
}

void
FFT_Function::initialize() {
	window = new HanningWindow(framesize());
	window->ref();
	xreal = new double[framesize()];
	ximag = new double[framesize()];
	Super::initialize();
	setAnalysis(
	    new FFTData(analysisLength(), npoints, framerate(), sampRate())
	);
}

int
FFT_Function::operator () (double* input, Data* frame) {
	int fsize = framesize();
	int i;
	for(i=0; i < fsize; i++) {
		*(xreal+i) = *input++;
		*(ximag+i) = 0.0;
	}
	window->applyTo(xreal);			// window the input
	fft(fsize, xreal, ximag);		// do the fft
	// calculate magnitude
	int nvals = fsize / 2;			// only convert and write first half
	for(i=0; i < nvals; i++)
		frame->set(hypot(xreal[i], ximag[i]), 0, i);
	return true;
}

// private static functions for doing transform

void
FFT_Function::fft(int nvalues, double* xreal, double* ximag) {
	static double logTwo = M_LN2;		/* log(2.0) */
	int l = int(log( (double)(abs(nvalues)) ) / logTwo + 0.001);
	int n = int(pow(2.0, double(l)) + 0.5);

	fftalgorithm( (int)(nvalues/abs(nvalues)),l,n,xreal,ximag);

	descramble( ((nvalues>=0)?n:1),l,n,xreal,ximag);
}

void
FFT_Function::fftalgorithm(int wexpsign, int l, int n, double* xr, double* xi) {
	static double PI_x_2 = 2.0 * M_PI;

/* printf("fft: wes,l,n = %d,%d,%d\n",wexpsign,l,n);
*/	int ia = n/2;
	int ib = 1;

	for (int i = 0; i < l; ++i) {
/* printf("level %d\n",i);
*/		int ic = 0;
		int id = ia;

		for (int k = 0; k < ib; ++k) {
			int ie = scramble( l,n,(int)(ic/ia) );
/* printf("ic = %d, ia = %d\n",ic,ia);
printf("W power = %d\n",ie);
*/			double wexp = -wexpsign*PI_x_2*(double)(ie)/(double)(n);
			double zr = cos(wexp);
			double zi = sin(wexp);

			for (int m = ic; m < id; ++m) {
				double ar = xr[m];
				double ai = xi[m];

				double br = xr[m+ia] * zr - xi[m+ia] * zi;
				double bi = xr[m+ia] * zi + xi[m+ia] * zr;

/* printf("x(%d) = x(%d) + W^(%d) * x(%d)\n",m,m,ie,m+ia);
*/				xr[m] = ar + br;
				xi[m] = ai + bi;

/* printf("x(%d) = x(%d) - W^(%d) * x(%d)\n",m+ia,m,ie,m+ia);
*/				xr[m+ia] = ar - br;
				xi[m+ia] = ai - bi;
			}

			ic += 2*ia;
			id += 2*ia;
		}

		ia /= 2;
		ib *= 2;
	}

}

void
FFT_Function::descramble(int scalediv, int l, int n, double* xreal, double* ximag) {
	for (int m1 = 0; m1 < n; ++m1) {
		int m2 = scramble(l,n,m1);
		if (m1 == m2) {
			xreal[m1] /= scalediv;
			ximag[m1] /= scalediv;
		}
		else if (m1 > m2) {
			double swap = xreal[m1]/scalediv;
			xreal[m1] = xreal[m2]/scalediv;
			xreal[m2] = swap;

			swap = ximag[m1]/scalediv;
			ximag[m1] = ximag[m2]/scalediv;
			ximag[m2] = swap;
		}
	}
}

int
FFT_Function::scramble(int l, int n, int min) {
	int mout = 0;
	int mplus = 1;
	int mminus = n/2;
	for (int i = 0; i < l; ++i) {
		if (min >= mminus) {
			mout += mplus;
			min -= mminus;
		}
		mplus *= 2;
		mminus /= 2;
	}
	return mout;
}