/* nag_wav_dim1_init (c09aac) Example Program.
*
* Copyright 2023 Numerical Algorithms Group.
*
* Mark 29.0, 2023.
*/
/* Pre-processor includes */
#include <math.h>
#include <nag.h>
#include <stdio.h>
#include <string.h>
int main(void) {
/* Constants */
Integer licomm = 100;
/*Integer scalar and array declarations */
Integer exit_status = 0;
Integer i, n, nf, nnz, nwc, nwlmax, ny;
Integer *dwtlev = 0, *icomm = 0;
NagError fail;
Nag_Wavelet wavnamenum;
Nag_WaveletMode modenum;
/*Double scalar and array declarations */
double *c = 0, *x = 0, *y = 0;
/*Character scalar and array declarations */
char mode[24], wavnam[20];
INIT_FAIL(fail);
printf("nag_wav_dim1_init (c09aac) Example Program Results\n\n");
fflush(stdout);
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read n - length of input data sequence */
scanf("%" NAG_IFMT "%*[^\n] ", &n);
if (!(x = NAG_ALLOC(n, double)) || !(y = NAG_ALLOC(n, double)) ||
!(icomm = NAG_ALLOC(licomm, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read Wavelet name (wavnam) and end mode (mode) */
scanf("%19s%23s%*[^\n] ", wavnam, mode);
/*
* nag_enum_name_to_value (x04nac).
* Converts NAG enum member name to value
*/
wavnamenum = (Nag_Wavelet)nag_enum_name_to_value(wavnam);
modenum = (Nag_WaveletMode)nag_enum_name_to_value(mode);
if (n >= 2) {
printf(" Parameters read from file :: \n");
printf(" Wavelet :%15s\n", wavnam);
printf(" End mode :%15s\n", mode);
printf(" N :%15" NAG_IFMT "\n\n", n);
/* Read data array and write it out */
printf("%s\n", " Input Data X :");
for (i = 0; i < n; i++) {
scanf("%lf ", &x[i]);
printf("%8.3f%s", x[i], (i + 1) % 8 ? " " : "\n");
}
scanf("%*[^\n] ");
printf("\n");
/*
* nag_wav_dim1_init (c09aac)
* Wavelet filter query
*/
nag_wav_dim1_init(wavnamenum, Nag_MultiLevel, modenum, n, &nwlmax, &nf,
&nwc, icomm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_wav_dim1_init (c09aac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
if (!(c = NAG_ALLOC(nwc, double)) ||
!(dwtlev = NAG_ALLOC(nwlmax + 1, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Perform Discrete Wavelet transform */
/*
* nag_wav_dim1_multi_fwd (c09ccc)
* one-dimensional multi-level discrete wavelet transform (mldwt)
*/
nag_wav_dim1_multi_fwd(n, x, nwc, c, nwlmax, dwtlev, icomm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_wav_dim1_multi_fwd (c09ccc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf(" Length of wavelet filter : %10" NAG_IFMT "\n", nf);
printf(" Number of Levels : %10" NAG_IFMT "\n",
nwlmax);
printf(" Number of coefficients in each level : \n");
for (i = 0; i < nwlmax + 1; i++)
printf(" %8" NAG_IFMT "%s", dwtlev[i], (i + 1) % 8 ? " " : "\n");
printf("\n");
printf(" Total number of wavelet coefficients :%10" NAG_IFMT "\n", nwc);
nnz = 0;
for (i = 0; i < nwlmax + 1; i++)
nnz = nnz + dwtlev[i];
printf("\n");
printf(" Wavelet coefficients C : \n");
for (i = 0; i < nnz; i++)
printf("%8.3f%s", c[i], (i + 1) % 8 ? " " : "\n");
printf("\n");
/* Reconstruct original data */
ny = n;
/*
* nag_wav_dim1_multi_inv (c09cdc)
* one-dimensional inverse multi-level discrete wavelet transform
* (imldwt)
*/
nag_wav_dim1_multi_inv(nwlmax, nwc, c, n, y, icomm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_wav_dim1_multi_inv (c09cdc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf("\n");
printf(" Reconstruction Y : \n");
for (i = 0; i < ny; i++)
printf("%8.3f%s", y[i], (i + 1) % 8 ? " " : "\n");
}
END:
NAG_FREE(c);
NAG_FREE(dwtlev);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(icomm);
return exit_status;
}