suseislet/apfilt.c at Main · vandanue/suseislet · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
/*
 * apfilt.c  --  All-pass plane-wave destruction filter coefficients
 *
 * NOTE: Uses double precision internally for numerical accuracy.
 * The filter coefficients implement Wilson's all-pass approximation
 * to a fractional-sample shift (slope p), computed using binomial
 * expansion. These are called by pwd.c (passfilter, aderfilter).
 */

#include "su.h"     /* err()        */
#include "apfilt.h"
#include <stdlib.h> /* malloc, free */

static int n;     /* filter order = 2*nw */
static double *b; /* binomial coefficients [n+1] */

void apfilt_init(int nw)
/*< initialize filter of order nw; call once before passfilter/aderfilter >*/
{
  int j, k;
  double bk;

  n = nw * 2;

  /* was: b = (double*) sf_alloc(n+1, sizeof(double)) */
  b = (double *)malloc((n + 1) * sizeof(double));
  if (!b)
    err("apfilt_init: malloc failed");

  /* Compute binomial base coefficients b[k] */
  for (k = 0; k <= n; k++) {
    bk = 1.0;
    for (j = 0; j < n; j++) {
      if (j < n - k) {
        bk *= (k + j + 1.0) / (2 * (2 * j + 1) * (j + 1));
      } else {
        bk *= 1.0 / (2 * (2 * j + 1));
      }
    }
    b[k] = bk;
  }
}

void apfilt_close(void)
/*< free allocated storage >*/
{
  free(b);
}

void passfilter(float p, float *a)
/*< compute all-pass filter coefficients for slope p >*/
/*  output: a[0..n] (n+1 coefficients) */
{
  int j, k;
  double ak;

  for (k = 0; k <= n; k++) {
    ak = b[k];
    for (j = 0; j < n; j++) {
      if (j < n - k) {
        ak *= (n - j - p);
      } else {
        ak *= (p + j + 1);
      }
    }
    a[k] = (float)ak; /* downcast to float for use in pwd.c / allp2.c */
  }
}

void aderfilter(float p, float *a)
/*< compute derivative of all-pass coefficients with respect to slope p >*/
/*  output: a[0..n] (n+1 coefficients) */
{
  int i, j, k;
  double ak, ai;

  for (k = 0; k <= n; k++) {
    ak = 0.0;
    for (i = 0; i < n; i++) {
      ai = -1.0;
      for (j = 0; j < n; j++) {
        if (j != i) {
          if (j < n - k) {
            ai *= (n - j - p);
          } else {
            ai *= (p + j + 1);
          }
        } else if (j < n - k) {
          ai *= (-1);
        }
      }
      ak += ai;
    }
    a[k] = (float)(ak * b[k]);
  }
}