25 #define REPORT { static ExeCounter ExeCount(__LINE__,19); ++ExeCount; } 35 long n4 = n * 4;
int sector = (int)floor( (
Real)n4 / (
Real)d + 0.5 );
37 if (sector < 0) {
REPORT sector = 3 - (3 - sector) % 4; }
38 else {
REPORT sector %= 4; }
43 case 0:
REPORT c = cos(ratio); s = sin(ratio);
break;
44 case 1:
REPORT c = -sin(ratio); s = cos(ratio);
break;
45 case 2:
REPORT c = -cos(ratio); s = -sin(ratio);
break;
46 case 3:
REPORT c = sin(ratio); s = -cos(ratio);
break;
56 const int gamma = after * before;
const int delta = now * after;
61 const int m = A.
Nrows() - gamma;
63 for (
int j = 0; j < now; j++)
66 Real* x1 = x;
Real* y1 = y; x += after; y += after;
67 for (
int ia = 0; ia < after; ia++)
71 cossin(-(j*after+ia), delta, r_arg, i_arg);
80 Real* a2 = m + a1;
Real* b2 = m + b1; a1 += after; b1 += after;
81 Real r_value = *a2;
Real i_value = *b2;
82 *x2 = r_value * r_arg - i_value * i_arg + *(a2-gamma);
83 *y2 = r_value * i_arg + i_value * r_arg + *(b2-gamma);
85 x2 += delta; y2 += delta;
94 Real* a2 = m + a1;
Real* b2 = m + b1; a1 += after; b1 += after;
95 Real r_value = *a2;
Real i_value = *b2;
96 int in = now-1;
while (in--)
101 a2 -= gamma; b2 -= gamma;
Real temp = r_value;
102 r_value = r_value * r_arg - i_value * i_arg + *a2;
103 i_value = temp * i_arg + i_value * r_arg + *b2;
105 *x2 = r_value; *y2 = i_value;
107 x2 += delta; y2 += delta;
127 const Real n = X.
Nrows(); X /= n; Y /= (-n);
135 const int n = U.
Nrows();
136 const int n2 = n / 2;
141 while (i--) { *a++ = *u++; *b++ = *u++; }
146 a = A.Store(); b = B.
Store();
149 Real* xn = x + n2;
Real* yn = y + n2;
151 *x++ = *a + *b; *y++ = 0.0;
152 *xn-- = *a++ - *b++; *yn-- = 0.0;
154 int j = -1; i = n2/2;
158 Real am = *a - *an;
Real ap = *a++ + *an--;
159 Real bm = *b - *bn;
Real bp = *b++ + *bn--;
160 Real samcbp = s * am + c * bp;
Real sbpcam = s * bp - c * am;
161 *x++ = 0.5 * ( ap + samcbp); *y++ = 0.5 * ( bm + sbpcam);
162 *xn-- = 0.5 * ( ap - samcbp); *yn-- = 0.5 * (-bm + sbpcam);
171 const int n21 = A.
Nrows();
172 if (n21 != B.
Nrows() || n21 == 0)
174 const int n2 = n21 - 1;
const int n = 2 * n2;
int i = n2 - 1;
178 Real* an = a + n2;
Real* bn = b + n2;
183 *x++ = hn * (*a + *an); *y++ = - hn * (*a - *an);
184 a++; an--; b++; bn--;
185 int j = -1; i = n2/2;
189 Real am = *a - *an;
Real ap = *a++ + *an--;
190 Real bm = *b - *bn;
Real bp = *b++ + *bn--;
191 Real samcbp = s * am - c * bp;
Real sbpcam = s * bp + c * am;
192 *x++ = hn * ( ap + samcbp); *y++ = - hn * ( bm + sbpcam);
193 *xn-- = hn * ( ap - samcbp); *yn-- = - hn * (-bm + sbpcam);
198 while (i--) { *u++ = *x++; *u++ = - *y++; }
208 const int n = U.
Nrows();
209 if (n != V.
Nrows() || n == 0)
211 if (n == 1) {
REPORT X = U; Y = V;
return; }
224 const int nextmx = 8;
225 int prime[8] = { 2,3,5,7,11,13,17,19 };
226 int after = 1;
int before = n;
int next = 0;
bool inzee =
true;
233 if (next < nextmx) {
REPORT now = prime[next]; }
234 b1 = before / now;
if (b1 * now == before) {
REPORT break; }
239 if (inzee) {
REPORT fftstep(A, B, X, Y, after, now, before); }
242 inzee = !inzee; after *= now;
258 const int n = U.
Nrows();
259 const int n2 = n / 2;
const int n4 = n * 4;
265 while (i--) { *a++ = *u++; *(--b) = *u++; }
277 *(++v) = xi * c + yi * s; *(--w) = xi * s - yi * c;
286 const int n = V.
Nrows();
287 const int n2 = n / 2;
const int n4 = n * 4;
const int n21 = n2 + 1;
294 int i = n2;
int k = 0;
299 *(++x) = vi * c + wi * s; *(++y) = vi * s - wi * c;
305 while (i--) { *u++ = *a++; *u++ = *(--b); }
313 const int n = U.
Nrows();
314 const int n2 = n / 2;
const int n4 = n * 4;
320 while (i--) { *a++ = *u++; *(--b) = -(*u++); }
332 *v++ = xi * s - yi * c; *(--w) = xi * c + yi * s;
341 const int n = V.
Nrows();
342 const int n2 = n / 2;
const int n4 = n * 4;
const int n21 = n2 + 1;
348 *x = *(--w); *y = 0.0;
349 int i = n2;
int k = 0;
354 *(++x) = vi * s + wi * c; *(++y) = - vi * c + wi * s;
360 while (i--) { *u++ = *a++; *u++ = -(*(--b)); }
368 const int n = V.
Nrows()-1;
369 const int n2 = n / 2;
const int n21 = n2 + 1;
374 Real vi = *v++; *x++ = vi; *y++ = 0.0;
375 Real sum1 = vi / 2.0;
Real sum2 = sum1; vi = *v++;
379 Real vi2 = *v++; sum1 += vi2 + vi; sum2 += vi2 - vi;
380 *x++ = vi2; vi2 = *v++; *y++ = vi - vi2; vi = vi2;
382 sum1 += vi; sum2 -= vi;
383 vi = *v; *x = vi; *y = 0.0; vi /= 2.0; sum1 += vi; sum2 += vi;
387 i = n2;
int k = 0; *u++ = sum1 / n2; *v-- = sum2 / n2;
390 Real s = sin(1.5707963267948966192 * (++k) / n2);
392 Real bz = (ai - bi) / 4 / s;
Real az = (ai + bi) / 2;
393 *u++ = az - bz; *v-- = az + bz;
403 V *= (V.
Nrows()-1)/2;
411 const int n = V.
Nrows()-1;
412 const int n2 = n / 2;
const int n21 = n2 + 1;
417 Real vi = *(++v); *x++ = 2 * vi; *y++ = 0.0;
419 while (i--) { *y++ = *(++v);
Real vi2 = *(++v); *x++ = vi2 - vi; vi = vi2; }
420 *x = -2 * vi; *y = 0.0;
424 i = n2;
int k = 0; *u++ = 0.0; *v-- = 0.0;
427 Real s = sin(1.5707963267948966192 * (++k) / n2);
429 Real az = (ai + bi) / 4 / s;
Real bz = (ai - bi) / 2;
430 *u++ = az - bz; *v-- = az + bz;
440 V *= (V.
Nrows()-1)/2;
449 if (m != V.
Nrows() || n != V.
Ncols() || m == 0 || n == 0)
453 for (i = 1; i <=
m; ++i)
456 X.
Row(i) = CVR.
t(); Y.
Row(i) = CVI.
t();
458 for (i = 1; i <= n; ++i)
void DST_inverse(const ColumnVector &V, ColumnVector &U)
void DST_II(const ColumnVector &U, ColumnVector &V)
void RealFFT(const ColumnVector &U, ColumnVector &X, ColumnVector &Y)
Miscellaneous exception (details in character string).
void DCT_II(const ColumnVector &U, ColumnVector &V)
void FFT(const ColumnVector &U, const ColumnVector &V, ColumnVector &X, ColumnVector &Y)
void FFT2(const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
void DCT_II_inverse(const ColumnVector &V, ColumnVector &U)
void FFT2I(const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
GetSubMatrix Column(int f) const
void FFTI(const ColumnVector &U, const ColumnVector &V, ColumnVector &X, ColumnVector &Y)
TransposedMatrix t() const
void DST(const ColumnVector &U, ColumnVector &V)
static void fftstep(ColumnVector &A, ColumnVector &B, ColumnVector &X, ColumnVector &Y, int after, int now, int before)
The usual rectangular matrix.
FloatVector FloatVector * a
Real trace(const BaseMatrix &B)
static void cossin(int n, int d, Real &c, Real &s)
static bool ar_1d_ft(int PTS, Real *X, Real *Y)
void RealFFTI(const ColumnVector &A, const ColumnVector &B, ColumnVector &U)
GetSubMatrix Row(int f) const
void DCT(const ColumnVector &U, ColumnVector &V)
void DCT_inverse(const ColumnVector &V, ColumnVector &U)
static bool CanFactor(int PTS)
void DST_II_inverse(const ColumnVector &V, ColumnVector &U)