20 using namespace NEWMAT;
30 Real f = 6.2831853071795864769/n;
31 for (
int j=1; j<=n; j++)
34 for (
int k=1; k<=n; k++)
36 Real theta = - (j-1) * (k-1) * f;
37 Real c = cos(theta);
Real s = sin(theta);
38 sumx += c *
a(k) - s * b(k); sumy += s *
a(k) + c * b(k);
40 x(j) = sumx; y(j) = sumy;
47 Real f = 6.2831853071795864769 / (4*n);
54 for (
int j=1; j<=n; j++) sum += cos(k1 * (2*j-1) * f) *
a(j);
61 for (
int j=1; j<=n; j++) sum += sin(k * (2*j-1) * f) *
a(j);
68 int n1 = a.
Nrows();
int n = n1 - 1;
70 Real f = 6.2831853071795864769 / (2*n);
77 for (
int j=2; j<=n; j++) sum += cos((j-1) * (k-1) * f) *
a(j);
78 c(k) = sum + (
a(1) + sign *
a(n1)) / 2.0;
85 for (
int j=2; j<=n; j++) sum += sin((j-1) * (k-1) * f) *
a(j);
95 if (m != V.
Nrows() || n != V.
Ncols() || m == 0 || n == 0)
98 const Real pi2 = atan(1.0) * 8.0;
99 for (
int i = 0; i <
m; ++i)
for (
int j = 0; j < n; ++j)
101 Real sumr = 0.0, sumi = 0.0;
102 for (
int k = 0; k <
m; ++k)
for (
int l = 0; l < n; ++l)
107 sumr += cs * U(k+1,l+1) - sn * V(k+1,l+1);
108 sumi += cs * V(k+1,l+1) + sn * U(k+1,l+1);
110 X(i+1,j+1) = sumr; Y(i+1,j+1) = sumi;
121 for (
int i=0; i<n; i++)
123 FFT(A, B, X, Y);
FFTI(X, Y, X, Y);
124 X = X - A; Y = Y - B;
130 Tracer et(
"Test RealFFT");
141 X2 = X1 - X.
Rows(1,n2); Y2 = Y1 - Y.
Rows(1,n2);
151 Tracer et(
"cf FFT and slow FT");
155 for (
int i=0; i<n; i++)
159 X = X - X1; Y = Y - Y1;
165 Tracer et(
"cf slow and fast DCT_II and DST_II");
182 for (
int i=0; i<n; i++) A1.
element(i) = mwc.
Next();
185 P = 0.0; Q = 0.0; P.
Rows(1,n) = A1;
188 for (
int k=0; k<n; k++)
190 Real arg = k * 6.2831853071795864769 / (4 * n);
191 B1(k+1) = P(k+1) * cos(arg) + Q(k+1) * sin(arg);
209 for (
int i=0; i<n; i++) A1.
element(i) = mwc.
Next();
212 P = 0.0; Q = 0.0; P.
Rows(1,n) = A1;
215 for (
int k=1; k<=n; k++)
217 Real arg = k * 6.2831853071795864769 / (4 * n);
218 B1(k) = P(k+1) * sin(arg) - Q(k+1) * cos(arg);
237 for (
int i=1; i<n; i++) A1.
element(i) = mwc.
Next();
261 for (
int i=0; i<=n; i++) A1.
element(i) = mwc.
Next();
265 P(1) /= 2.0; P(n+1) /= 2.0;
280 Tracer et(
"cf slow and fast DCT and DST");
284 for (
int i=0; i<=n; i++) A.
element(i) = mwc.
Next();
294 Tracer et(
"cf FFT2 and slow FT2");
297 Matrix A(m,n), B(m,n), X, Y, X1, Y1;
298 for (
int i=0; i<
m; i++)
for (
int j=0; j<n; j++)
302 X = X - X1; Y = Y - Y1;
304 FFT2I(X1, Y1, X1, Y1);
314 Tracer et(
"Fifteenth test of Matrix package");
319 for (i = 1; i <=12; i++)
323 + .2 * cos(6.2831853071795864769 * 4.0 * i1 / 12)
324 + .3 * sin(6.2831853071795864769 * 3.0 * i1 / 12);
326 + .5 * sin(6.2831853071795864769 * 2.0 * i1 / 12)
327 + .4 * cos(6.2831853071795864769 * 1.0 * i1 / 12);
330 A(1) -= 8.4; A(3) -= 3.0; A(5) -= 1.2; A(9) -= 1.2; A(11) += 3.0;
331 B(1) -= 10.8; B(2) -= 2.4; B(4) += 1.8; B(10) -= 1.8; B(12) -= 2.4;
374 for (i = 1; i <=12; i++)
378 + .2 * cos(6.2831853071795864769 * 4.0 * i1 / 12)
379 + .3 * sin(6.2831853071795864769 * 3.0 * i1 / 12);
381 + .5 * sin(6.2831853071795864769 * 2.0 * i1 / 12)
382 + .4 * cos(6.2831853071795864769 * 1.0 * i1 / 12);
385 A(1) -= 8.4; A(3) -= 3.0; A(5) -= 1.2; A(9) -= 1.2; A(11) += 3.0;
386 B(1) -= 10.8; B(2) -= 2.4; B(4) += 1.8; B(10) -= 1.8; B(12) -= 2.4;
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 SlowFT2(const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
static void SlowFT(const ColumnVector &a, const ColumnVector &b, ColumnVector &x, ColumnVector &y)
virtual void ReSize(int m, int n)
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)
void Clean(Matrix &A, Real c)
void FFTI(const ColumnVector &U, const ColumnVector &V, ColumnVector &X, ColumnVector &Y)
void DST(const ColumnVector &U, ColumnVector &V)
The usual rectangular matrix.
static void test9(int m, int n)
FloatVector FloatVector * a
Real trace(const BaseMatrix &B)
void RealFFTI(const ColumnVector &A, const ColumnVector &B, ColumnVector &U)
static void SlowDTT_II(const ColumnVector &a, ColumnVector &c, ColumnVector &s)
Real sum(const BaseMatrix &B)
void DCT(const ColumnVector &U, ColumnVector &V)
void DCT_inverse(const ColumnVector &V, ColumnVector &U)
void Print(const Matrix &X)
static void SlowDTT(const ColumnVector &a, ColumnVector &c, ColumnVector &s)
void DST_II_inverse(const ColumnVector &V, ColumnVector &U)
GetSubMatrix Rows(int f, int l) const