/* -*- MaTX -*-
 *  
 *  NAME
 *      idft()      - Inverse discrete Fourier transform
 *      idft_plot() - plot the data of inverse discrete Fourier transform
 *
 *  SYNOPSIS
 *      X = idft(Y, N)
 *         Array X;
 *         Array Y;
 *
 *      X = idft(Y, N)
 *         Array X;
 *         Array Y;
 *         Integer N;
 *
 *  DESCRIPTION
 *      Y = idft(X, N) calculates the inverse discrete fourier transfrom of X
 *      idft_plot(X, N) 
 *	
 *      X(n) = 1/N sum_{k=0}^{N-1} Y(k) W^{- k n}   n = 0,1,...,N-1
 *      W = exp^{-j*2PI/N}
 *
 *      Y(k) = Y(j*k * 2*PI/(N*dt)),  dt := Sampling Interval
 *
 *   SEE ALSO
 *      dft
 *
 */

Func Array idft(Y_, N, ...)
	Array Y_;
	Integer N;
{
	Integer n, k;
	Real Nr;
	Complex j;
	Array X, Y, W0, Wn;

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

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

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

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

	return X/N;
}

Func void idft_plot(Y, N, ...)
	Array Y;
	Integer N;
{
	Integer win;
	Array X;

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

	if (nargs == 1) {
		X = idft(Y);
	} else if (nargs == 2) {
		X = idft(Y, N);
	}

	win = mgplot_cur_win();
	mgplot_xlabel(win, "time");
	mgplot_ylabel(win, "");
	mgplot_title(win, "idft plot");
	mgplot(win, Re(X));
}