/* -*- MaTX -*-
 *  
 *  NAME
 *      dft()      - Discrete Fourier Transform
 *      dft_plot() - Plot Discrete Fourier Transform results
 *
 *  SYNOPSIS
 *      Y = dft(X)
 *         Array Y;
 *         Array X;
 *
 *      Y = dft(X, N)
 *         Array Y;
 *         Array X;
 *         Integer N;
 *
 *       dft_plot(X)
 *         Array X;
 *
 *       dft_plot(X, dt)
 *         Array X;
 *         Real dt;
 *
 *       dft_plot(X, dt, N)
 *         Array X;
 *         Real dt;
 *         Integer N;
 *
 *  DESCRIPTION
 *      Y = dft(X, N) calculates the discrete fourier transfrom of X as
 *      Y(k) = sum_{n=0}^{N-1} X(n) W^{k n}   k = 0,1,...,N-1
 *
 *      W = exp^{-j*2PI/N}
 *
 *      Y(k) = Y(j*k * 2*PI/(N*T))  T := Sampling Interval
 *
 *      dft_plot(X, dt, N) plots the data of discrete fourier transform
 *
 *  SEE ALSO
 *      idft
 */

Func Array dft(X_, N, ...)
	Array X_;
	Integer N;
{
	Integer n, k;
	Complex j;
	Array X, Y, W0, Wn;

	error(nargchk(1, 2, nargs, "dft"));

	X = X_;

	if (nargs == 1) {
		N = Cols(X);
	} else {
		n = Cols(X);
		if (N <= n) {
			X = X(*, 1:N);
		} else if (N > n) {
			X = [X, Z(Rows(X), N - n)];
		}
	}

	j = (0,1);
	W0 = 2*PI*[0:N-1]'/N;
	Wn = exp(-j * W0);

	Y = (Z(X),:);
	for (k = 0; k <= N-1; k++) {
		Y(:,k+1) = Matrix(X) * Matrix(Wn^k);
	}

	return Y;
}

Func void dft_plot(X, dt, N, ...)
	Array X;
	Real dt;	// Sampling interval [s]
	Integer N;
{
	Integer win;
	Array W, Y, Pyy;

	error(nargchk(1, 3, nargs, "dft_plot"));

	if (nargs < 3) {
		N = Cols(X);
	}
	if (nargs == 1) {
		dt = 1.0;
	}

	Y = dft(X, N);

	W = 1.0/dt * [0:(N/2 - 1)]/N; // Frequencies [Hz]
	Pyy = Re(Y * conj(Y))/N;	  // Energy at various frequencies

	win = mgplot_cur_win();
	mgplot_xlabel(win, "Frequency [Hz]");
	mgplot_ylabel(win, "Power");
	mgplot_title(win, "dft plot");
	mgplot(1, W, Pyy(1:N/2));
}