


(*
 * Multi-dimensional Fourier transform.
 * Given a multi-dimensional array
 *	a[[i, j, ...]] with 1 <= i <= n, 1 <= j <= m, ...
 * return an array b (of the same shape) such that b[[s, t, ...]] is
 *	Sum[
 *		a[[i, j, ...]] zeta^((s-1)(i-1)) gamma^((t-1)(j-1)) ...,
 *		{i, n},
 *		{j, m},
 *		...
 *	] / Sqrt[n m ...]
 * where
 *	zeta = Exp[2 Pi I / n], gamma = Exp[2 Pi I / m], ...
 * MultiInverseFourier performs the inverse of this transformation, which
 * simply the same as above except that
 *	zeta = Exp[-2 Pi I / n], gamma = Exp[-2 Pi I / m], ...
 *)

Begin["System`"]

Unprotect[Fourier, InverseFourier]

Begin["Fourier`Private`"]

Fourier`Private`HighDimension

HighDimension[a_]:= (Length[Dimensions[a]]>2)

Fourier[a_?HighDimension] :=
	Block[{
		Fourier`Private`rank = Length[Dimensions[a]],
		Fourier`Private`shift,
		Fourier`Private`res = a
		},
		shift = RotateLeft[Range[rank]];
		Do[
			res = Map[
				Fourier,
				res,
				{rank - 1}
			];
			res = Transpose[res, shift],
			{rank}
		];
		res
	] 

Protect[Fourier]

InverseFourier[a_?HighDimension] :=
	Block[{
		Fourier`Private`rank = Length[Dimensions[a]],
		Fourier`Private`shift,
		Fourier`Private`res = a
		},
		shift = RotateLeft[Range[rank]];
		Do[
			res = Map[
				InverseFourier,
				res,
				{rank - 1}
			];
			res = Transpose[res, shift],
			{rank}
		];
		res
	]

Protect[InverseFourier]

End[]

End[];

