share/numeric/rtest_rfft.mac

   1 /* Test real_fft routines */
   2
   3 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<=tol) then true else abse);
   4 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<=tol) then true else abse);
   5
   6 (compare_real_fft(tol, x) :=
   7   block([rft: real_fft(x),
   8          ft: fft(x),
   9          maperror: false,
  10          maxdiff],
  11     maxdiff: lreduce('max, abs(rft - ft)),
  12     closeto(maxdiff, tol)),
  13  0);
  14 0;
  15
  16 (compare_inverse_real_fft(tol, x) :=
  17   block([rft: inverse_real_fft(real_fft(x)),
  18          maxdiff],
  19     maxdiff: lreduce('max, abs(x - rft)),
  20     closeto(maxdiff, tol)),
  21  0);
  22 0;
  23
  24
  25 /* verify that real_fft does not modify its argument */
  26
  27 (foo0 : [1, 2, 4, 8], 0);
  28 0;
  29
  30 (kill (foo), foo : foo0, real_fft (foo), foo);
  31 ''foo0;
  32
  33 (foo : make_array (any, 4), fillarray (foo, foo0), real_fft (foo), listarray (foo));
  34 ''foo0;
  35
  36 (kill (foo), array (foo, 3), fillarray (foo, foo0), real_fft (foo), listarray (foo));
  37 ''foo0;
  38
  39 /* Verify basic computation */
  40
  41 real_fft([]);
  42 [];
  43
  44 real_fft([123]);
  45 [123.0];
  46
  47 real_fft([1,2]);
  48 [1.5, -0.5];
  49
  50 errcatch(real_fft([1,2,3]));
  51 [];
  52
  53 compare_real_fft(0, makelist(k, k, 1, 4));
  54 true;
  55
  56 compare_real_fft(1.2413e-16, makelist(k, k, 1, 16));
  57 true;
  58
  59 compare_real_fft(9.8712e-14, makelist(k, k, 1, 1024));
  60 true;
  61
  62 compare_real_fft(3.8703e-13, makelist(k, k, 1, 4096));
  63 true;
  64
  65 inverse_real_fft([]);
  66 [];
  67
  68 inverse_real_fft([123]);
  69 [123.0];
  70
  71 inverse_real_fft([1.5 ,-0.5]);
  72 [1.0, 2.0];
  73
  74 compare_inverse_real_fft(0, makelist(k, k, 1, 4));
  75 true;
  76
  77 compare_inverse_real_fft(1.7764e-15, makelist(k, k, 1, 16));
  78 true;
  79
  80 compare_inverse_real_fft(3.9791e-13, makelist(k, k, 1, 1024));
  81 true;
  82
  83 compare_inverse_real_fft(2.2738e-12, makelist(k, k, 1, 4096));
  84 true;
  85
  86 /* verify that bf_real_fft does not modify its argument */
  87
  88 (foo0 : [1, 2, 4, 8], 0);
  89 0;
  90
  91 (kill (foo), foo : foo0, bf_real_fft (foo), foo);
  92 ''foo0;
  93
  94 (foo : make_array (any, 4), fillarray (foo, foo0), bf_real_fft (foo), listarray (foo));
  95 ''foo0;
  96
  97 (kill (foo), array (foo, 3), fillarray (foo, foo0), bf_real_fft (foo), listarray (foo));
  98 ''foo0;
  99
 100 /* Verify basic computation */
 101
 102 bf_real_fft([]);
 103 [];
 104
 105 bf_real_fft([123]);
 106 [123b0];
 107
 108 bf_real_fft([1,2]);
 109 [1.5b0, -0.5b0];
 110
 111 errcatch(bf_real_fft([1,2,3]));
 112 [];
 113
 114 (compare_bf_real_fft(tol, x) :=
 115   block([rft: bf_real_fft(x),
 116          ft: bf_fft(x),
 117          maperror: false,
 118          maxdiff],
 119     maxdiff: lreduce('max, abs(rft - ft)),
 120     closeto(maxdiff, tol)),
 121  0);
 122 0;
 123
 124 (compare_bf_inverse_real_fft(tol, x) :=
 125   block([rft: bf_inverse_real_fft(bf_real_fft(x)),
 126          maxdiff],
 127     maxdiff: lreduce('max, abs(x - rft)),
 128     closeto(maxdiff, tol)),
 129  0);
 130 0;
 131
 132 compare_bf_real_fft(0, makelist(k, k, 1, 4));
 133 true;
 134
 135 compare_bf_real_fft(2.7756b-17, makelist(k, k, 1, 16));
 136 true;
 137
 138 compare_bf_real_fft(5.7034b-15, makelist(k, k, 1, 1024));
 139 true;
 140
 141 compare_bf_real_fft(2.4145b-14, makelist(k, k, 1, 4096));
 142 true;
 143
 144 bf_inverse_real_fft([]);
 145 [];
 146
 147 bf_inverse_real_fft([1.5 ,-0.5]);
 148 [1b0, 2b0];
 149
 150 compare_bf_inverse_real_fft(0, makelist(k, k, 1, 4));
 151 true;
 152
 153 compare_bf_inverse_real_fft(3.3088b-24, makelist(k, k, 1, 16));
 154 true;
 155
 156 compare_bf_inverse_real_fft(8.4704b-22, makelist(k, k, 1, 1024));
 157 true;
 158
 159 compare_bf_inverse_real_fft(5.0822b-21, makelist(k, k, 1, 4096));
 160 true;
 161
 162 (errcatch(real_fft([1,2,3])), error);
 163 ["real_fft: size of input must be a power of 2, not ~M",3];
 164
 165 (errcatch(inverse_real_fft([1,2,3,4])), error);
 166 ["inverse_real_fft: input length must be one more than a power of two, not ~M",4];
 167
 168 (errcatch(bf_real_fft([1,2,3])), error);
 169 ["bf_fft: size of input must be a power of 2, not ~M",3];
 170
 171 (errcatch(bf_inverse_real_fft([1,2,3,4])), error);
 172 ["bf_inverse_real_fft: input length must be one more than a power of two, not ~M",4];
 173
 174