Source

stockwelltransform / stockwell / st.c

Full commit
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <fftw3.h>

char *Wisfile = NULL;
char *Wistemplate = "%s/.fftwis";
#define WISLEN 8

void set_wisfile(void)
{
  char *home;

  if (Wisfile) return;
  home = getenv("HOME");
  Wisfile = (char *)malloc(strlen(home) + WISLEN + 1);
  sprintf(Wisfile, Wistemplate, home);
}

/* Convert frequencies in Hz into rows of the ST, given sampling rate and
   length. */

int st_freq(double f, int len, double srate)
{
  return floor(f * len / srate + .5);
}

static double gauss(int n, int m);

/* Stockwell transform of the real array data. The len argument is the
   number of time points, and it need not be a power of two. The lo and hi
   arguments specify the range of frequencies to return, in Hz. If they are
   both zero, they default to lo = 0 and hi = len / 2. The result is
   returned in the complex array result, which must be preallocated, with
   n rows and len columns, where n is hi - lo + 1. For the default values of
   lo and hi, n is len / 2 + 1.

   Clm: I don't think that lo and hi are actually specified in Hz but instead are in
   units related to the length of the signal to be analyzed


*/

void st(int len, int lo, int hi, double *data, double *result)
{
  int i, k, n, l2;
  double s, *p;
  FILE *wisdom;
  static int planlen = 0;
  static double *g;
  static fftw_plan p1, p2;
  static fftw_complex *h, *H, *G;

  /* Check for frequency defaults. */

  if (lo == 0 && hi == 0) {
    hi = len / 2;
  }

  /* Keep the arrays and plans around from last time, since this
     is a very common case. Reallocate them if they change. */

  if (len != planlen && planlen > 0) {
    fftw_destroy_plan(p1);
    fftw_destroy_plan(p2);
    fftw_free(h);
    fftw_free(H);
    fftw_free(G);
    free(g);
    planlen = 0;
  }

  if (planlen == 0) {
    planlen = len;
    h = fftw_malloc(sizeof(fftw_complex) * len);
    H = fftw_malloc(sizeof(fftw_complex) * len);
    G = fftw_malloc(sizeof(fftw_complex) * len);
    g = (double *)malloc(sizeof(double) * len);

    /* Get any accumulated wisdom. */

    set_wisfile();
    wisdom = fopen(Wisfile, "r");
    if (wisdom) {
      fftw_import_wisdom_from_file(wisdom);
      fclose(wisdom);
    }

    /* Set up the fftw plans. */

    p1 = fftw_plan_dft_1d(len, h, H, FFTW_FORWARD, FFTW_MEASURE);
    p2 = fftw_plan_dft_1d(len, G, h, FFTW_BACKWARD, FFTW_MEASURE);

    /* Save the wisdom. */

    wisdom = fopen(Wisfile, "w");
    if (wisdom) {
      fftw_export_wisdom_to_file(wisdom);
      fclose(wisdom);
    }
  }

  /* Convert the input to complex. Also compute the mean. */

  s = 0.;
  memset(h, 0, sizeof(fftw_complex) * len);
  for (i = 0; i < len; i++) {
    h[i][0] = data[i];
    s += data[i];
  }
  s /= len;

  /* FFT. */

  fftw_execute(p1); /* h -> H */

  /* Hilbert transform. The upper half-circle gets multiplied by
     two, and the lower half-circle gets set to zero.  The real axis
     is left alone. */

  l2 = (len + 1) / 2;
  for (i = 1; i < l2; i++) {
    H[i][0] *= 2.;
    H[i][1] *= 2.;
  }
  l2 = len / 2 + 1;
  for (i = l2; i < len; i++) {
    H[i][0] = 0.;
    H[i][1] = 0.;
  }

  /* Fill in rows of the result. */

  p = result;

  /* The row for lo == 0 contains the mean. */

  n = lo;
  if (n == 0) {
    for (i = 0; i < len; i++) {
      *p++ = s;
      *p++ = 0.;
    }
    n++;
  }

  /* Subsequent rows contain the inverse FFT of the spectrum
     multiplied with the FFT of scaled gaussians. */

  while (n <= hi) {

    /* Scale the FFT of the gaussian. Negative frequencies
       wrap around. */

    g[0] = gauss(n, 0);
    l2 = len / 2 + 1;
    for (i = 1; i < l2; i++) {
      g[i] = g[len - i] = gauss(n, i);
    }

    for (i = 0; i < len; i++) {
      s = g[i];
      k = n + i;
      if (k >= len) k -= len;
      G[i][0] = H[k][0] * s;
      G[i][1] = H[k][1] * s;
    }

    /* Inverse FFT the result to get the next row. */

    fftw_execute(p2); /* G -> h */
    for (i = 0; i < len; i++) {
      *p++ = h[i][0] / len;
      *p++ = h[i][1] / len;
    }

    /* Go to the next row. */

    n++;
  }
}

/* This is the Fourier Transform of a Gaussian. */

static double gauss(int n, int m)
{
  return exp(-2. * M_PI * M_PI * m * m / (n * n));
}

/* Inverse Stockwell transform. */

void ist(int len, int lo, int hi, double *data, double *result)
{
  int i, n, l2;
  double *p;
  FILE *wisdom;
  static int planlen = 0;
  static fftw_plan p2;
  static fftw_complex *h, *H;

  /* Check for frequency defaults. */

  if (lo == 0 && hi == 0) {
    hi = len / 2;
  }

  /* Keep the arrays and plans around from last time, since this
     is a very common case. Reallocate them if they change. */

  if (len != planlen && planlen > 0) {
    fftw_destroy_plan(p2);
    fftw_free(h);
    fftw_free(H);
    planlen = 0;
  }

  if (planlen == 0) {
    planlen = len;
    h = fftw_malloc(sizeof(fftw_complex) * len);
    H = fftw_malloc(sizeof(fftw_complex) * len);

    /* Get any accumulated wisdom. */

    set_wisfile();
    wisdom = fopen(Wisfile, "r");
    if (wisdom) {
      fftw_import_wisdom_from_file(wisdom);
      fclose(wisdom);
    }

    /* Set up the fftw plans. */

    p2 = fftw_plan_dft_1d(len, H, h, FFTW_BACKWARD, FFTW_MEASURE);

    /* Save the wisdom. */

    wisdom = fopen(Wisfile, "w");
    if (wisdom) {
      fftw_export_wisdom_to_file(wisdom);
      fclose(wisdom);
    }
  }

  /* Sum the complex array across time. */

  memset(H, 0, sizeof(fftw_complex) * len);
  p = data;
  for (n = lo; n <= hi; n++) {
    for (i = 0; i < len; i++) {
      H[n][0] += *p++;
      H[n][1] += *p++;
    }
  }

  /* Invert the Hilbert transform. */

  l2 = (len + 1) / 2;
  for (i = 1; i < l2; i++) {
    H[i][0] /= 2.;
    H[i][1] /= 2.;
  }
  l2 = len / 2 + 1;
  for (i = l2; i < len; i++) {
    H[i][0] = H[len - i][0];
    H[i][1] = -H[len - i][1];
  }

  /* Inverse FFT. */

  fftw_execute(p2); /* H -> h */
  p = result;
  for (i = 0; i < len; i++) {
    *p++ = h[i][0] / len;
  }
}

/* This does just the Hilbert transform. */

void hilbert(int len, double *data, double *result)
{
  int i, l2;
  double *p;
  FILE *wisdom;
  static int planlen = 0;
  static fftw_plan p1, p2;
  static fftw_complex *h, *H;

  /* Keep the arrays and plans around from last time, since this
     is a very common case. Reallocate them if they change. */

  if (len != planlen && planlen > 0) {
    fftw_destroy_plan(p1);
    fftw_destroy_plan(p2);
    fftw_free(h);
    fftw_free(H);
    planlen = 0;
  }

  if (planlen == 0) {
    planlen = len;
    h = fftw_malloc(sizeof(fftw_complex) * len);
    H = fftw_malloc(sizeof(fftw_complex) * len);

    /* Get any accumulated wisdom. */

    set_wisfile();
    wisdom = fopen(Wisfile, "r");
    if (wisdom) {
      fftw_import_wisdom_from_file(wisdom);
      fclose(wisdom);
    }

    /* Set up the fftw plans. */

    p1 = fftw_plan_dft_1d(len, h, H, FFTW_FORWARD, FFTW_MEASURE);
    p2 = fftw_plan_dft_1d(len, H, h, FFTW_BACKWARD, FFTW_MEASURE);

    /* Save the wisdom. */

    wisdom = fopen(Wisfile, "w");
    if (wisdom) {
      fftw_export_wisdom_to_file(wisdom);
      fclose(wisdom);
    }
  }

  /* Convert the input to complex. */

  memset(h, 0, sizeof(fftw_complex) * len);
  for (i = 0; i < len; i++) {
    h[i][0] = data[i];
  }

  /* FFT. */

  fftw_execute(p1); /* h -> H */

  /* Hilbert transform. The upper half-circle gets multiplied by
     two, and the lower half-circle gets set to zero.  The real axis
     is left alone. */

  l2 = (len + 1) / 2;
  for (i = 1; i < l2; i++) {
    H[i][0] *= 2.;
    H[i][1] *= 2.;
  }
  l2 = len / 2 + 1;
  for (i = l2; i < len; i++) {
    H[i][0] = 0.;
    H[i][1] = 0.;
  }

  /* Inverse FFT. */

  fftw_execute(p2); /* H -> h */

  /* Fill in the rows of the result. */

  p = result;
  for (i = 0; i < len; i++) {
    *p++ = h[i][0] / len;
    *p++ = h[i][1] / len;
  }
}