home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / PRELOAD.PAK / FOURIER.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  1.4 KB  |  78 lines
  1.  
  2.  
  3.  
  4. (*
  5.  * Multi-dimensional Fourier transform.
  6.  * Given a multi-dimensional array
  7.  *    a[[i, j, ...]] with 1 <= i <= n, 1 <= j <= m, ...
  8.  * return an array b (of the same shape) such that b[[s, t, ...]] is
  9.  *    Sum[
  10.  *        a[[i, j, ...]] zeta^((s-1)(i-1)) gamma^((t-1)(j-1)) ...,
  11.  *        {i, n},
  12.  *        {j, m},
  13.  *        ...
  14.  *    ] / Sqrt[n m ...]
  15.  * where
  16.  *    zeta = Exp[2 Pi I / n], gamma = Exp[2 Pi I / m], ...
  17.  * MultiInverseFourier performs the inverse of this transformation, which
  18.  * simply the same as above except that
  19.  *    zeta = Exp[-2 Pi I / n], gamma = Exp[-2 Pi I / m], ...
  20.  *)
  21.  
  22. Begin["System`"]
  23.  
  24. Unprotect[Fourier, InverseFourier]
  25.  
  26. Begin["Fourier`Private`"]
  27.  
  28. Fourier`Private`HighDimension
  29.  
  30. HighDimension[a_]:= (Length[Dimensions[a]]>2)
  31.  
  32. Fourier[a_?HighDimension] :=
  33.     Block[{
  34.         Fourier`Private`rank = Length[Dimensions[a]],
  35.         Fourier`Private`shift,
  36.         Fourier`Private`res = a
  37.         },
  38.         shift = RotateLeft[Range[rank]];
  39.         Do[
  40.             res = Map[
  41.                 Fourier,
  42.                 res,
  43.                 {rank - 1}
  44.             ];
  45.             res = Transpose[res, shift],
  46.             {rank}
  47.         ];
  48.         res
  49.     ] 
  50.  
  51. Protect[Fourier]
  52.  
  53. InverseFourier[a_?HighDimension] :=
  54.     Block[{
  55.         Fourier`Private`rank = Length[Dimensions[a]],
  56.         Fourier`Private`shift,
  57.         Fourier`Private`res = a
  58.         },
  59.         shift = RotateLeft[Range[rank]];
  60.         Do[
  61.             res = Map[
  62.                 InverseFourier,
  63.                 res,
  64.                 {rank - 1}
  65.             ];
  66.             res = Transpose[res, shift],
  67.             {rank}
  68.         ];
  69.         res
  70.     ]
  71.  
  72. Protect[InverseFourier]
  73.  
  74. End[]
  75.  
  76. End[];
  77.  
  78.