suseislet/pwd.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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
/*
 * pwd.c  --  Plane-wave definition (banded matrix representation)
 */

#include <stdbool.h> /* bool, true, false                          */
#include <stdlib.h>  /* malloc, free                               */

#include "su.h" /* ealloc1float, ealloc2float, err            */
#include "apfilt.h" /* apfilt_init, apfilt_close, passfilter      */
#include "pwd.h"


/* Full struct definition (private to this file) */
struct Pwd {
  int n, na;
  float **a; /* filter coefficients: a[na][n] */
  float *b;  /* temporary filter buffer: b[na] */
};

/* ================================================================== */
pwd pwd_init(int n1, int nw)
/*< initialize plane-wave filter of order nw for n1-sample traces >*/
{
  pwd w;

  w = (pwd)malloc(sizeof(*w)); /* was: sf_alloc(1,sizeof(*w)) */
  if (!w)
    err("pwd_init: malloc failed");

  w->n = n1;
  w->na = 2 * nw + 1;

  w->a = ealloc2float(n1, w->na); /* was: sf_floatalloc2(n1,w->na) */
  w->b = ealloc1float(w->na);     /* was: sf_floatalloc(w->na)     */

  apfilt_init(nw); /* still from apfilt.c */

  return w;
}

/* ================================================================== */
void pwd_close(pwd w)
/*< free allocated storage >*/
{
  free(w->a[0]);
  free(w->a);
  free(w->b);
  free(w);
  apfilt_close(); /* still from apfilt.c */
}

/* ================================================================== */
void pwd_define(bool adj, pwd w, const float *pp, float *diag, float **offd)
/*< fill banded matrix with plane-wave destruction contributions >*/
/*  diag and offd are modified (PWD terms are ADDED to existing values) */
{
  int i, j, k, m, n, nw;
  float am, aj;

  nw = (w->na - 1) / 2;
  n = w->n;

  /* Compute all-pass filter coefficients for each sample position */
  for (i = 0; i < n; i++) {
    passfilter(pp[i], w->b); /* still from apfilt.c */

    for (j = 0; j < w->na; j++) {
      /* adjoint reverses the filter */
      w->a[j][i] = adj ? w->b[w->na - 1 - j] : w->b[j];
    }
  }

  /* Add A^T A contributions to banded matrix */
  for (i = 0; i < n; i++) {
    for (j = 0; j < w->na; j++) {
      k = i + j - nw;
      if (k >= nw && k < n - nw) {
        aj = w->a[j][k];
        diag[i] += aj * aj;
      }
    }
    for (m = 0; m < 2 * nw; m++) {
      for (j = m + 1; j < w->na; j++) {
        k = i + j - nw;
        if (k >= nw && k < n - nw) {
          aj = w->a[j][k];
          am = w->a[j - m - 1][k];
          offd[m][i] += am * aj;
        }
      }
    }
  }
}

/* ================================================================== */
void pwd_set(bool adj, pwd w, float *inp, float *out, float *tmp)
/*< multiply by the PWD matrix (pure array operations, no RSF calls) >*/
{
  int i, j, k, n, nw;

  nw = (w->na - 1) / 2;
  n = w->n;

  if (adj) {
    for (i = 0; i < n; i++)
      tmp[i] = 0.0f;
    for (i = 0; i < n; i++) {
      for (j = 0; j < w->na; j++) {
        k = i + j - nw;
        if (k >= nw && k < n - nw)
          tmp[k] += w->a[j][k] * out[i];
      }
    }
    for (i = 0; i < n; i++)
      inp[i] = 0.0f;
    for (i = nw; i < n - nw; i++) {
      for (j = 0; j < w->na; j++) {
        k = i + j - nw;
        inp[k] += w->a[j][i] * tmp[i];
      }
    }
  } else {
    for (i = 0; i < n; i++)
      tmp[i] = 0.0f;
    for (i = nw; i < n - nw; i++) {
      for (j = 0; j < w->na; j++) {
        k = i + j - nw;
        tmp[i] += w->a[j][i] * inp[k];
      }
    }
    for (i = 0; i < n; i++) {
      out[i] = 0.0f;
      for (j = 0; j < w->na; j++) {
        k = i + j - nw;
        if (k >= nw && k < n - nw)
          out[i] += w->a[j][k] * tmp[k];
      }
    }
  }
}