home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / numana01.zip / SRC / TESTS / FFTTEST.MOD < prev    next >
Text File  |  1996-07-31  |  13KB  |  428 lines
  1. MODULE FFTtest;
  2.  
  3.         (********************************************************)
  4.         (*                                                      *)
  5.         (*              Test of module Fourier                  *)
  6.         (*                                                      *)
  7.         (*  Programmer:         P. Moylan                       *)
  8.         (*  Last edited:        31 July 1996                    *)
  9.         (*  Status:             Working                         *)
  10.         (*                                                      *)
  11.         (*      Observation: although the FFT is clearly a      *)
  12.         (*      lot faster than using the direct definition     *)
  13.         (*      of a DFT (i.e. procedure SlowFT), its           *)
  14.         (*      accuracy improvement is smaller than I would    *)
  15.         (*      have expected.  This is probably because I'm    *)
  16.         (*      using iterative methods to calculate the        *)
  17.         (*      complex exponentials.                           *)
  18.         (*                                                      *)
  19.         (********************************************************)
  20.  
  21. FROM Fourier IMPORT
  22.     (* proc *)  FFT3, FFT8, SlowFT, RealFFT;
  23.  
  24. FROM MiscM2 IMPORT
  25.     (* const*)  PI,
  26.     (* proc *)  SelectWindow, PressAnyKey, WriteString, WriteLn,
  27.                 WriteCard, WriteLongReal, Sin, Cos;
  28.  
  29. FROM LongComplexMath IMPORT
  30.     (* proc *)  abs;
  31.  
  32. IMPORT Cx;
  33.  
  34. (*
  35. FROM Windows IMPORT
  36.     (* type *)  Window, Colour, FrameType, DividerType,
  37.     (* proc *)  OpenWindow, CloseWindow;
  38. *)
  39.  
  40. (************************************************************************)
  41.  
  42. TYPE
  43.     CxArray = ARRAY [0..1023] OF LONGCOMPLEX;
  44.  
  45. VAR
  46.     (* To avoid stack overflow problems, we use static storage for our  *)
  47.     (* big arrays.  (Heap allocation would also have been a viable      *)
  48.     (* option.)                                                         *)
  49.  
  50.     TestData: CxArray;
  51.     Results: CxArray;
  52.     SpareCopy: CxArray;
  53.     RealData: ARRAY [0..1023] OF LONGREAL;
  54.  
  55. (************************************************************************)
  56. (*                      PLAUSIBILITY TESTS                              *)
  57. (************************************************************************)
  58.  
  59. PROCEDURE TransformAndWrite (N: CARDINAL);
  60.  
  61.     (* Compares results of slow and fast methods. *)
  62.  
  63.     VAR j: CARDINAL;
  64.  
  65.     BEGIN
  66.         WriteString ("     Slow result             FFT result");
  67.         WriteString ("      Magnitude of difference");
  68.         WriteLn;
  69.         SlowFT (TRUE, N, TestData, Results);
  70.         FFT8 (TRUE, N, TestData);
  71.         FOR j := 0 TO N-1 DO
  72.             Cx.Write (Results[j], 8);
  73.             WriteString ("     ");
  74.             Cx.Write (TestData[j], 8);
  75.             WriteString ("     ");
  76.             WriteLongReal (abs(Results[j] - TestData[j]), 8);
  77.             WriteLn;
  78.         END (*FOR*);
  79.         PressAnyKey;
  80.     END TransformAndWrite;
  81.  
  82. (************************************************************************)
  83.  
  84. PROCEDURE PlausibilityTest;
  85.  
  86.     (* FFT of some sample data arrays. *)
  87.  
  88.     CONST testsize = 16;
  89.  
  90.     VAR j: [0..testsize-1];  theta: LONGREAL;
  91.         (*w: Window;*)
  92.  
  93.     BEGIN
  94.         (*
  95.         OpenWindow (w, black, brown, 0, 24, 0, 79, simpleframe, nodivider);
  96.         SelectWindow (w);
  97.         *)
  98.         WriteString ("COMPARISON OF FFT WITH DIRECT TRANSFORM");
  99.         WriteLn;
  100.  
  101.         (* Test on a constant function. *)
  102.  
  103.         WriteString ("Constant function, 1 data points.");
  104.         WriteLn;
  105.         FOR j := 0 TO testsize-1 DO
  106.             TestData[j] := CMPLX (1.0, 0.0);
  107.         END (*FOR*);
  108.         TransformAndWrite (1);
  109.         WriteString ("Constant function, 2 data points.");
  110.         WriteLn;
  111.         FOR j := 0 TO testsize-1 DO
  112.             TestData[j] := CMPLX (1.0, 0.0);
  113.         END (*FOR*);
  114.         TransformAndWrite (2);
  115.  
  116.         WriteString ("Constant function, 4 data points.");
  117.         WriteLn;
  118.         FOR j := 0 TO testsize-1 DO
  119.             TestData[j] := CMPLX (1.0, 0.0);
  120.         END (*FOR*);
  121.         TransformAndWrite (4);
  122.  
  123.         (* Test on a sine function. *)
  124.  
  125.         WriteString ("Cosine function, ");
  126.         WriteCard (testsize);
  127.         WriteString (" data points.");
  128.         WriteLn;
  129.         FOR j := 0 TO testsize-1 DO
  130.             theta := PI*VAL(LONGREAL,j)/VAL(LONGREAL,testsize DIV 2);
  131.             TestData[j] := CMPLX (Cos(theta), 0.0);
  132.         END (*FOR*);
  133.         TransformAndWrite (testsize);
  134.  
  135.         (* Test on a ramp function. *)
  136.  
  137.         WriteString ("Ramp function.");
  138.         WriteLn;
  139.         FOR j := 0 TO testsize-1 DO
  140.             TestData[j] := CMPLX (VAL(LONGREAL,j), 0.0);
  141.         END (*FOR*);
  142.         TransformAndWrite (testsize);
  143.  
  144.         (*CloseWindow (w);*)
  145.  
  146.     END PlausibilityTest;
  147.  
  148. (************************************************************************)
  149. (*                      TEST USING THE INVERSE TRANSFORM                *)
  150. (************************************************************************)
  151.  
  152. TYPE TransformProc = PROCEDURE (BOOLEAN, CARDINAL, VAR ARRAY OF LONGCOMPLEX);
  153.  
  154. PROCEDURE OneComparison (N: CARDINAL;  Transform, InvTransform: TransformProc;
  155.                                         title: ARRAY OF CHAR);
  156.  
  157.     (* Does a transform with "Transform", then an inverse transform     *)
  158.     (* with "InvTransform", and summarises the errors.                  *)
  159.  
  160.     VAR j: CARDINAL;  error, max, sum: LONGREAL;
  161.  
  162.     BEGIN
  163.         (* Take a copy of the data. *)
  164.  
  165.         FOR j := 0 TO N-1 DO
  166.             Results[j] := TestData[j];
  167.         END (*FOR*);
  168.  
  169.         (* Do a direct and then an inverse transform. *)
  170.  
  171.         Transform (TRUE, N, Results);
  172.         InvTransform (FALSE, N, Results);
  173.  
  174.         (* Work out the errors. *)
  175.  
  176.         max := 0.0;  sum := 0.0;
  177.         FOR j := 0 TO N-1 DO
  178.             error := abs (Results[j] - TestData[j]);
  179.             IF error > max THEN max := error END(*IF*);
  180.             sum := sum + error;
  181.         END (*FOR*);
  182.  
  183.         WriteString (title);
  184.         FOR j := 0 TO 16 - HIGH(title) DO WriteString (" ") END(*FOR*);
  185.         WriteLongReal (max, 8);
  186.         FOR j := 0 TO 7 DO WriteString (" ") END(*FOR*);
  187.         WriteLongReal (sum, 8);
  188.         WriteLn;
  189.  
  190.     END OneComparison;
  191.  
  192. (************************************************************************)
  193.  
  194. PROCEDURE DoubleComparison (N: CARDINAL);
  195.  
  196.     (* Compares results of the available methods. *)
  197.  
  198.     (* Conclusion so far: FFT3 seems to be a little more accurate than  *)
  199.     (* FFT8 (although I suspect that FFT8 is faster).  I'd like to      *)
  200.     (* check out whether further variants are worth following up.       *)
  201.  
  202.     VAR j: CARDINAL;  max, sum, error: LONGREAL;
  203.  
  204.     BEGIN
  205.         WriteString ("   Method        max error");
  206.         WriteString ("     sum of errors");
  207.         WriteLn;
  208.  
  209.         OneComparison (N, FFT3, FFT3, "FFT3/FFT3");
  210.         OneComparison (N, FFT3, FFT8, "FFT3/FFT8");
  211.         OneComparison (N, FFT8, FFT3, "FFT8/FFT3");
  212.         OneComparison (N, FFT8, FFT8, "FFT8/FFT8");
  213.  
  214.         (* For the slow FT, the calling conventions are different       *)
  215.         (* because we don't have an in-place algorithm.                 *)
  216.  
  217.         (* Take a copy of the data. *)
  218.  
  219.         FOR j := 0 TO N-1 DO
  220.             SpareCopy[j] := TestData[j];
  221.         END (*FOR*);
  222.  
  223.         (* Do a direct and then an inverse transform. *)
  224.  
  225.         SlowFT (TRUE, N, TestData, Results);
  226.         SlowFT (FALSE, N, Results, TestData);
  227.  
  228.         (* Work out the errors. *)
  229.  
  230.         max := 0.0;  sum := 0.0;
  231.         FOR j := 0 TO N-1 DO
  232.             error := abs (SpareCopy[j] - TestData[j]);
  233.             IF error > max THEN max := error END(*IF*);
  234.             sum := sum + error;
  235.         END (*FOR*);
  236.  
  237.         WriteString ("Slow             ");
  238.         WriteLongReal (max, 8);
  239.         FOR j := 0 TO 7 DO WriteString (" ") END(*FOR*);
  240.         WriteLongReal (sum, 8);
  241.         WriteLn;
  242.  
  243.         PressAnyKey;
  244.  
  245.     END DoubleComparison;
  246.  
  247. (************************************************************************)
  248.  
  249. PROCEDURE DoubleTransformTest;
  250.  
  251.     (* FFT followed by inverse FFT. *)
  252.  
  253.     CONST testsize = 128;       (* I had to cut this down, the slow     *)
  254.                                 (* transform was taking too long.       *)
  255.         k = testsize DIV 2;
  256.  
  257.     VAR j: [0..testsize-1];  theta: LONGREAL;
  258.         (*w: Window;*)
  259.  
  260.     BEGIN
  261.         (*
  262.         OpenWindow (w, black, green, 0, 24, 0, 79, simpleframe, nodivider);
  263.         SelectWindow (w);
  264.         *)
  265.         WriteString ("DIRECT TRANSFORM FOLLOWED BY INVERSE TRANSFORM");
  266.         WriteLn;
  267.  
  268.         (* Test on a constant function. *)
  269.  
  270.         WriteLn;
  271.         WriteString ("Constant function, ");
  272.         WriteCard (testsize);
  273.         WriteString (" data points.");
  274.         WriteLn;
  275.         FOR j := 0 TO testsize-1 DO
  276.             TestData[j] := CMPLX (1.0, 0.0);
  277.         END (*FOR*);
  278.         DoubleComparison (testsize);
  279.  
  280.         (* Test on a sine function. *)
  281.  
  282.         WriteLn;
  283.         WriteString ("Cosine function, ");
  284.         WriteCard (testsize);
  285.         WriteString (" data points.");
  286.         WriteLn;
  287.         FOR j := 0 TO testsize-1 DO
  288.             theta := PI*VAL(LONGREAL,j)/VAL(LONGREAL,testsize DIV 2);
  289.             TestData[j] := CMPLX (Cos(theta), 0.0);
  290.         END (*FOR*);
  291.         DoubleComparison (testsize);
  292.  
  293.         (* Test on some other simple functions. *)
  294.  
  295.         WriteLn;
  296.         WriteString ("Ramp function.");
  297.         WriteLn;
  298.         FOR j := 0 TO testsize-1 DO
  299.             TestData[j] := CMPLX (VAL(LONGREAL,j), 0.0);
  300.         END (*FOR*);
  301.         DoubleComparison (testsize);
  302.  
  303.         WriteLn;
  304.         WriteString ("Square wave.");
  305.         WriteLn;
  306.         FOR j := 0 TO k-1 DO
  307.             TestData[j] := CMPLX (1.0, 0.0);
  308.             TestData[j+k] := CMPLX (-1.0, 0.0);
  309.         END (*FOR*);
  310.         DoubleComparison (testsize);
  311.  
  312.         (*CloseWindow (w);*)
  313.  
  314.     END DoubleTransformTest;
  315.  
  316. (************************************************************************)
  317. (*                      TEST OF REAL FFT                                *)
  318. (************************************************************************)
  319.  
  320. PROCEDURE DoRealTest (N: CARDINAL);
  321.  
  322.     (* Performs transform and writes results. *)
  323.  
  324.     VAR j: CARDINAL;
  325.  
  326.     BEGIN
  327.         WriteString ("Real part         Imag part");
  328.         WriteLn;
  329.         RealFFT (TRUE, N, RealData);
  330.         WriteLongReal (RealData[0], 8);
  331.         WriteString ("        ");
  332.         WriteLongReal (0.0, 8);
  333.         WriteLn;
  334.         FOR j := 1 TO (N DIV 2)-1 DO
  335.             WriteLongReal (RealData[2*j], 8);
  336.             WriteString ("        ");
  337.             WriteLongReal (RealData[2*j + 1], 8);
  338.             WriteLn;
  339.         END (*FOR*);
  340.         WriteLongReal (RealData[1], 8);
  341.         WriteString ("        ");
  342.         WriteLongReal (0.0, 8);
  343.         WriteLn;
  344.         PressAnyKey;
  345.     END DoRealTest;
  346.  
  347. (************************************************************************)
  348.  
  349. PROCEDURE RealTest;
  350.  
  351.     (* FFT of some sample data arrays. *)
  352.  
  353.     CONST testsize = 16;
  354.  
  355.     VAR j: [0..testsize-1];  theta: LONGREAL;
  356.         (*w: Window;*)
  357.  
  358.     BEGIN
  359.         (*
  360.         OpenWindow (w, black, brown, 0, 24, 0, 79, simpleframe, nodivider);
  361.         SelectWindow (w);
  362.         *)
  363.         WriteString ("TEST OF REAL FFT");
  364.         WriteLn;
  365.  
  366.         (* Test on a constant function. *)
  367.  
  368.         WriteString ("Constant function, 2 data points.");
  369.         WriteLn;
  370.         FOR j := 0 TO testsize-1 DO
  371.             RealData[j] := 1.0;
  372.         END (*FOR*);
  373.         DoRealTest (2);
  374.  
  375.         WriteString ("Constant function, 4 data points.");
  376.         WriteLn;
  377.         FOR j := 0 TO testsize-1 DO
  378.             RealData[j] := 1.0;
  379.         END (*FOR*);
  380.         DoRealTest (4);
  381.  
  382.         (* Test on a sine function. *)
  383.  
  384.         WriteString ("Cosine function, ");
  385.         WriteCard (testsize);
  386.         WriteString (" data points.");
  387.         WriteLn;
  388.         FOR j := 0 TO testsize-1 DO
  389.             theta := PI*VAL(LONGREAL,j)/VAL(LONGREAL,testsize DIV 2);
  390.             RealData[j] := Cos(theta);
  391.         END (*FOR*);
  392.         DoRealTest (testsize);
  393.  
  394.         (* Test on a ramp function. *)
  395.  
  396.         WriteString ("Ramp function.");
  397.         WriteLn;
  398.         FOR j := 0 TO testsize-1 DO
  399.             RealData[j] := VAL(LONGREAL,j);
  400.         END (*FOR*);
  401.         DoRealTest (testsize);
  402.  
  403.         (* Do an inverse transform on this last result. *)
  404.  
  405.         WriteString ("Inverse transform of this last result.");
  406.         WriteLn;
  407.         RealFFT (FALSE, testsize, RealData);
  408.         FOR j := 0 TO testsize-1 DO
  409.             WriteLongReal (RealData[j], 12);
  410.             WriteLn;
  411.         END (*FOR*);
  412.         PressAnyKey;
  413.  
  414.         (*CloseWindow (w);*)
  415.  
  416.     END RealTest;
  417.  
  418. (************************************************************************)
  419. (*                              MAIN PROGRAM                            *)
  420. (************************************************************************)
  421.  
  422. BEGIN
  423.     PlausibilityTest;
  424.     DoubleTransformTest;
  425.     RealTest;
  426. END FFTtest.
  427. 
  428.