/* nag_fft_multiple_qtr_cosine(c06hdc) Example Program * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. */ #include #include #include #include #define MMAX 5 #define NMAX 20 #define X(I,J) x[(I)*n + (J)] int main(void) { double trig[2*NMAX], x[MMAX*NMAX]; Integer i, j, m, n; Vprintf("c06hdc Example Program Results\n"); Vscanf(" %*[^\n]"); /* Skip heading in data file */ while (scanf("%ld %ld", &m, &n) != EOF) if (m <= MMAX && n <= NMAX) { Vscanf(" %*[^\n]"); /* Skip text in data file */ Vscanf(" %*[^\n]"); for (i = 0; i < m; ++i) for (j = 0; j < n; ++j) Vscanf("%lf", &X(i,j)); Vprintf("\nOriginal data values\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < n; ++j) Vprintf(" %10.4f%s", X(i,j), (j%7==6 && j!=n-1 ? "\n" : "")); Vprintf("\n"); } /* Initialise trig array */ c06gzc(n, trig, NAGERR_DEFAULT); /* Compute transform */ c06hdc(Nag_ForwardTransform, m, n, x, trig, NAGERR_DEFAULT); Vprintf("\nDiscrete quarter-wave Fourier cosine transforms\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < n; ++j) Vprintf(" %10.4f%s", X(i,j), (j%7==6 && j!=n-1 ? "\n" : "")); Vprintf("\n"); } /* Compute inverse transform */ c06hdc(Nag_BackwardTransform, m, n, x, trig, NAGERR_DEFAULT); Vprintf("\nOriginal data as restored by inverse transform\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < n; ++j) Vprintf(" %10.4f%s", X(i,j), (j%7==6 && j!=n-1 ? "\n" : "")); Vprintf("\n"); } } else Vfprintf(stderr,"\nInvalid value of m or n.\n"); return EXIT_SUCCESS; }