seminaive.c

Go to the documentation of this file.

 
#include "s2kit/seminaive.h"
 
#include <math.h>
#include <stdio.h>
#include <string.h>
 
#include <fftw3.h>
 
#include "s2kit/cospml.h"
 
void InvDLTSemi(double* coeffs, const int bw, const int m, double* result, double* trans_cos_pml_table,
                double* sin_values, double* workspace, fftw_plan* plan) {
    int size = 2 * bw;
 
    // for paranoia, zero out arrays
    memset(workspace, 0, sizeof(double) * size);
    memset(result, 0, sizeof(double) * size);
 
    double* fcos = workspace;
    double* trans_tableptr = trans_cos_pml_table;
 
    double coeff = 0.5 / sqrt(bw);
    /* 
        main loop - compute each value of fcos
 
        Note that all zeroes have been stripped out of the
        `trans_cos_pml_table`, so indexing is somewhat complicated.
    */
    for (int i = 0; i < bw; ++i) {
        if (i == (bw - 1) && m % 2) {
            fcos[bw - 1] = 0.0;
            break;
        }
 
        double* assoc_offset;
        if (i > m)
            assoc_offset = coeffs + (i - m) + (m % 2);
        else
            assoc_offset = coeffs + (i % 2);
 
        int rowsize = Transpose_RowSize(i, m, bw);
 
        double value = 0.;
        for (int j = 0; j < rowsize; ++j)
            value += assoc_offset[2 * j] * trans_tableptr[j];
 
        fcos[i] = value * coeff; // scale coefficients prior to taking inverse DCT
 
        trans_tableptr += rowsize;
    }
 
    // special coeff for the first one
    fcos[0] /= (sqrt(size) * coeff);
 
    /*
        now we have the cosine series for the result,
        so now evaluate the cosine series at `2*bw` Chebyshev nodes
    */
 
    // take the inverse dct
    // Note that I am using the guru interface
    fftw_execute_r2r(*plan, fcos, result);
 
    if (!(m % 2))
        return;
 
    // if m is odd, then need to multiply by sin(x) at Chebyshev nodes
    for (int i = 0; i < size; ++i)
        result[i] *= sin_values[i];
}
 
void DLTSemi(double* data, const int bw, const int m, double* result, double* workspace, double* cos_pml_table,
             double* weights, fftw_plan* plan) {
    // TODO add check 0 <= `m` < `bw` (and to any other func)
    int size = 2 * bw;
 
    double* weighted_data = workspace;
    double* cos_data = weighted_data + size;
 
    // apply quadrature weights to the data and compute the cosine transform
    if (m % 2)
        for (int i = 0; i < size; ++i)
            weighted_data[i] = data[i] * weights[2 * bw + i];
    else
        for (int i = 0; i < size; ++i)
            weighted_data[i] = data[i] * weights[i];
 
    // smooth the weighted signal
    fftw_execute_r2r(*plan, weighted_data, cos_data);
 
    // normalize
    cos_data[0] *= M_SQRT1_2;
    double coeff = 1. / sqrt(2. * size);
    for (int i = 0; i < size; ++i)
        cos_data[i] *= coeff;
 
    /*
        Do the projections.
        Note that the `cos_pml_table` has had all the zeroes
        stripped out so the indexing is complicated somewhat
    */
    int toggle = 0;
    for (int i = m; i < bw; ++i) {
        double* pml_ptr = cos_pml_table + TableOffset(m, i);
 
        double value = 0.;
        for (int j = 0; j < (i / 2); ++j)
            value += cos_data[(2 * j) + toggle] * pml_ptr[j];
 
        if (!((i - m) % 2) || !(m % 2))
            value += cos_data[(2 * (i / 2)) + toggle] * pml_ptr[(i / 2)];
 
        result[i - m] = value;
 
        toggle = (toggle + 1) % 2;
    }
}