fourier.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
/* (C) Copr. 1986-92 Numerical Recipes Software #.,. */
 
/* Algorithms for Fourier Transform Spectral Methods */
 
#define _USE_MATH_DEFINES
 
#include <math.h>
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
 
// to print diagnostic statements
// #define DEBUG_SPCTRM
// otherwise
#undef DEBUG_SPCTRM
 
#undef myfabs
#if MAC
#if __GNUC__ < 4
#define myfabs std::fabs
#else
#define myfabs ::fabs
#endif
#else
#define myfabs fabs
#endif
 
// Apple LLVM version 5.1 (clang-503.0.38) generates "unsequenced modification warning".
// #define SQUARE(a) ((x_=(a)) == 0.0 ? 0.0 : x_*x_)
static inline double SQUARE(double a) {
    return a * a;
}
 
extern "C" {
void hoc_execerror(const char*, const char*);
extern void nrn_exit(int);
}  // extern "C"
 
#include "nrngsl_real_radix2.cpp"
#include "nrngsl_hc_radix2.cpp"
 
void nrngsl_realft(double* data, unsigned long n, int direction) {
    if (direction == 1) {
        nrngsl_fft_real_radix2_transform(data, 1, n);
    } else {
        nrngsl_fft_halfcomplex_radix2_inverse(data, 1, n);
    }
}
 
/*
  convlv()  -- convolution of a read data set with a response function
  the respns function must be stored in wrap-around order
*/
 
/* frequency domain format copy from gsl fft to NRC fft */
void nrn_gsl2nrc(double* x, double* y, unsigned long n) {
    unsigned long i, n2;
    n2 = n / 2;
    y[0] = x[0];
    if (n < 2) {
        return;
    }
    y[1] = x[n2];
    for (i = 1; i < n2; ++i) {
        y[2 * i] = x[i];
        y[2 * i + 1] = -x[n - i];
    }
}
void nrn_nrc2gsl(double* x, double* y, unsigned long n) {
    double xn, xn2;
    unsigned long i, n2;
    n2 = n / 2;
    xn = n;
    xn2 = .5 * xn;
    y[0] = x[0] * xn2;
    if (n < 2) {
        return;
    }
    y[n2] = x[1] * xn2;
    for (i = 1; i < n2; ++i) {
        y[i] = x[2 * i] * xn2;
        y[n - i] = -x[2 * i + 1] * xn2;
    }
}
 
void nrn_convlv(double* data,
                unsigned long n,
                double* respns,
                unsigned long m,
                int isign,
                double* ans) {
    // data and respns are modified.
    unsigned long i;
    double scl, x_;
 
    int n2 = n / 2;
    for (i = 1; i <= (m - 1) / 2; i++) {
        respns[n - i] = respns[m - i];
    }
    for (i = (m + 1) / 2; i < n - (m - 1) / 2; i++) {
        respns[i] = 0.0;
    }
    nrngsl_realft(data, n, 1);
    nrngsl_realft(respns, n, 1);
    ans[0] = data[0] * respns[0];
    for (i = 1; i < n2; ++i) {
        if (isign == 1) {
            ans[i] = data[i] * respns[i] - data[n - i] * respns[n - i];
            ans[n - i] = data[i] * respns[n - i] + data[n - i] * respns[i];
        } else if (isign == -1) {
            if ((scl = SQUARE(ans[i - 1]) + SQUARE(ans[i])) == 0.0) {
                hoc_execerror("Deconvolving at response zero in nrn_convlv", 0);
            }
            scl *= 2.0;
            ans[i] = (data[i] * respns[i] + data[n - i] * respns[n - i]) / scl;
            ans[i] = (data[i] * respns[n - i] - data[n - i] * respns[i]) / scl;
        } else {
            hoc_execerror("No meaning for isign in nrn_convlv", 0);
        }
    }
    ans[n2] = data[n2] * respns[n2];
    nrngsl_realft(ans, n, -1);
}
 
 
/*
  correl()  -- correlation of two real data sets
  see http://www.cv.nrao.edu/course/astr534/FourierTransforms.html SF16
*/
void nrn_correl(double* x, double* y, unsigned long n, double* z) {
    nrngsl_realft(x, n, 1);
    nrngsl_realft(y, n, 1);
    z[0] = x[0] * y[0];
    int n2 = n / 2;
    for (int i = 1; i < n2; ++i) {
        z[i] = x[i] * y[i] + x[n - i] * y[n - i];
        // with following sign, produces same result as old NRC code
        z[n - i] = y[i] * x[n - i] - y[n - i] * x[i];
    }
    z[n2] = x[n2] * y[n2];
    nrngsl_realft(z, n, -1);
}
 
 
/*
  spctrm()  -- power spectrum using fourier transforms
  modified to take array rather than data stream.
*/
 
#define WINDOW(j, a, b) (1.0 - myfabs((((j) -1) - (a)) * (b))) /* Bartlett */
 
// void nrn_spctrm(double* data, double* p, int m, int k)
void nrn_spctrm(double* data, double* psd, int setsize, int numsegpairs) {
    int j, k, cx, n;
    double a, ainv, wfac, x_;
    double* fftv;
 
    // 0 is index of first meaningful element of power spectral density
    for (j = 0; j <= setsize - 1; j++)
        psd[j] = 0.0;  // zero out any prior result
 
    // calc factor used to correct for window's effect on psd
    a = setsize;
    ainv = 1.0 / a;
 
    n = 2 * setsize;
    wfac = 0.0;
    for (j = 1; j <= n; j++)
        wfac += SQUARE(WINDOW(j, a, 1 / a));
 
#ifdef DEBUG_SPCTRM
    printf("setsize %d, numsegpairs %d\n", setsize, numsegpairs);
    printf("=============\n");
    printf("window values\n");
    for (j = 1; j <= n; j++)
        printf("%d  %f\n", j, WINDOW(j, a, 1 / a));
    printf("initial wfac %f\n", wfac);
#endif
 
    // storage for a data segment
    fftv = (double*) malloc((size_t) (n * sizeof(double)));
 
    cx = 0;  // data index
    for (k = 1; k <= numsegpairs * 2; k++) {
#ifdef DEBUG_SPCTRM
        printf("==============\n");
        printf("segment %d\n", k);
#endif
 
        // fill fftv with a segment of data
        for (j = 0; j < n; j++)
            fftv[j] = data[cx++];
#ifdef DEBUG_SPCTRM
        for (j = 0; j < n; j++)
            printf("datum %d  %f\n", j, fftv[j]);
        printf("--------------\n");
#endif
 
        cx -= setsize;  // next segment overlaps by this much
 
        // apply window
        for (j = 1; j <= n; j++)
            fftv[j - 1] *= WINDOW(j, a, 1 / a);
 
        // apply fft
        nrngsl_realft(fftv, n, 1);
#ifdef DEBUG_SPCTRM
        printf("\nfftv:\n");
        for (j = 0; j < n; j++)
            printf("%f  ", fftv[j]);
        printf("\n");
#endif
 
        // add to psd
        psd[0] += SQUARE(fftv[0]);
        for (j = 1; j < setsize; j++) {
            psd[j] += (SQUARE(fftv[j]) + SQUARE(fftv[n - j]));
        }
    }
 
    wfac = 1 / (wfac * n * numsegpairs);
    for (j = 0; j < setsize; j++)
        psd[j] *= wfac;
    psd[0] *= 0.5;
 
#ifdef DEBUG_SPCTRM
    printf("1/wfac for all but DC = %f\n", 1. / wfac);
    printf("1/wfac for DC = %f\n", 2. / wfac);
    printf("final results--\n");
    for (j = 0; j < setsize; j++)
        printf("%d  %f\n", j, psd[j]);
    printf("\ndone\n");
#endif
 
    free(fftv);
}
#undef WINDOW