/* nag_fft_multiple_sine(c06hac) 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)*row_len + (J)] int main(void) { double trig[2*NMAX], x[MMAX*NMAX]; Integer i, j, m, n, row_len; Vprintf("c06hac Example Program Results\n"); Vscanf(" %*[^\n]"); /* Skip heading in data file */ while (scanf("%ld %ld", &m, &n) != EOF) if (m <= MMAX && n <= NMAX) { row_len = n - 1; Vscanf(" %*[^\n]"); /* Skip text in data file */ Vscanf(" %*[^\n]"); for (i = 0; i < m; ++i) for (j = 0; j < row_len; ++j) Vscanf("%lf", &X(i,j)); Vprintf("\nOriginal data values\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < row_len; ++j) Vprintf(" %10.4f%s", X(i,j), (j%7==6 && j!=row_len-1 ? "\n" : "")); Vprintf("\n"); } c06gzc(n, trig, NAGERR_DEFAULT); /* Initialise trig array */ c06hac(m, n, x, trig, NAGERR_DEFAULT); /* Compute transform */ Vprintf("\nDiscrete Fourier sine transforms\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < row_len; ++j) Vprintf(" %10.4f%s", X(i,j), (j%7==6 && j!=row_len-1 ? "\n" : "")); Vprintf("\n"); } c06hac(m, n, x, trig, NAGERR_DEFAULT); /* Compute inverse transform */ Vprintf("\nOriginal data as restored by inverse transform\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < row_len; ++j) Vprintf(" %10.4f%s", X(i,j), (j%7==6 && j!=row_len-1 ? "\n" : "")); Vprintf("\n"); } } else Vfprintf(stderr,"\nInvalid value of m or n.\n"); return EXIT_SUCCESS; }