test_conv_semi_memo.c

Go to the documentation of this file.

 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#include <fftw3.h>
 
#include "s2kit/FST_semi_memo.h"
#include "s2kit/cospml.h"
 
int main(int argc, char** argv) {
    if (argc < 5) {
        fprintf(stdout, "Usage: test_conv_semi_memo signal_file filter_file output_file bw\n");
        exit(0);
    }
 
    int bw = atoi(argv[4]);
    int size = 2 * bw;
 
    int cutoff = bw; // seminaive all orders
    int legendreSize = Reduced_Naive_TableSize(bw, cutoff) + Reduced_SpharmonicTableSize(bw, cutoff);
 
    double* rsignal = (double*)malloc(sizeof(double) * size * size);
    double* rfilter = (double*)malloc(sizeof(double) * size * size);
 
    // read signal and filter
    fprintf(stdout, "Reading signal file...\n");
    FILE* fp = fopen(argv[1], "r");
    for (int i = 0; i < size * size; ++i)
        fscanf(fp, "%lf", rsignal + i);
    fclose(fp);
 
    fprintf(stdout, "Reading filter file...\n");
    fp = fopen(argv[2], "r");
    for (int i = 0; i < size * size; ++i)
        fscanf(fp, "%lf", rfilter + i);
    fclose(fp);
 
    // the imaginary parts are zeros,
    // since the data are strictly real-valued
    double* isignal = (double*)calloc((size_t)(size * size), sizeof(double));
    double* ifilter = (double*)calloc((size_t)(size * size), sizeof(double));
 
    double* rresult = (double*)malloc(sizeof(double) * size * size);
    double* iresult = (double*)malloc(sizeof(double) * size * size);
    double* workspace = (double*)malloc(sizeof(double) * (2 * legendreSize + 12 * bw * bw + 12 * bw));
 
    fprintf(stdout, "Calling ConvOn2SphereSemiMemo()\n");
    ConvOn2SphereSemiMemo(rsignal, isignal, rfilter, ifilter, rresult, iresult, bw, workspace);
 
    // convolving real functions results in real output,
    // so no need to write the imaginary array
    fprintf(stdout, "Writing output file...\n");
    fp = fopen(argv[3], "w");
    for (int i = 0; i < size * size; ++i)
        fprintf(fp, "%.16f\n", rresult[i]);
    fclose(fp);
 
    free(workspace);
    free(iresult);
    free(rresult);
    free(ifilter);
    free(rfilter);
    free(isignal);
    free(rsignal);
 
    return 0;
}