← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/signal/firfilter.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/signal/firfilter.c
 
1
#include <string.h>
 
2
#include <stdlib.h>
 
3
#include "sigtools.h"
 
4
 
 
5
#include "Python.h"                 /* only needed for defining unsigned types or not */
 
6
#include "Numeric/arrayobject.h"    
 
7
 
 
8
static int elsizes[] = {sizeof(char),
 
9
                        sizeof(unsigned char),
 
10
                        sizeof(signed char),
 
11
                        sizeof(short),
 
12
#ifdef PyArray_UNSIGNED_TYPES
 
13
                        sizeof(unsigned short),
 
14
#endif
 
15
                        sizeof(int),
 
16
#ifdef PyArray_UNSIGNED_TYPES
 
17
                        sizeof(unsigned int),
 
18
#endif
 
19
                        sizeof(long),
 
20
                        sizeof(float),
 
21
                        sizeof(double),
 
22
                        2*sizeof(float),
 
23
                        2*sizeof(double),
 
24
                        0};
 
25
 
 
26
typedef void (OneMultAddFunction) (char *, char *, char *);
 
27
 
 
28
#define MAKE_ONEMULTADD(fname, type) \
 
29
static void fname ## _onemultadd(char *sum, char *term1, char *term2) { \
 
30
  (*((type *) sum)) += (*((type *) term1)) * \
 
31
  (*((type *) term2)); return; }
 
32
 
 
33
#ifdef PyArray_UNSIGNED_TYPES
 
34
MAKE_ONEMULTADD(USHORT, unsigned short)
 
35
MAKE_ONEMULTADD(UINT, unsigned int)
 
36
#endif
 
37
MAKE_ONEMULTADD(UCHAR, unsigned char)
 
38
MAKE_ONEMULTADD(SCHAR, signed char)
 
39
MAKE_ONEMULTADD(SHORT, short)
 
40
MAKE_ONEMULTADD(INT, int)
 
41
MAKE_ONEMULTADD(LONG, long)
 
42
MAKE_ONEMULTADD(FLOAT, float)
 
43
MAKE_ONEMULTADD(DOUBLE, double)
 
44
 
 
45
#ifdef __GNUC__
 
46
MAKE_ONEMULTADD(CFLOAT, __complex__ float)
 
47
MAKE_ONEMULTADD(CDOUBLE, __complex__ double)
 
48
#else
 
49
#define MAKE_C_ONEMULTADD(fname, type) \
 
50
static void fname ## _onemultadd(char *sum, char *term1, char *term2) { \
 
51
  ((type *) sum)[0] += ((type *) term1)[0] * ((type *) term2)[0] \
 
52
    - ((type *) term1)[1] * ((type *) term2)[1]; \
 
53
  ((type *) sum)[1] += ((type *) term1)[0] * ((type *) term2)[1] \
 
54
    + ((type *) term1)[1] * ((type *) term2)[0]; \
 
55
  return; }
 
56
MAKE_C_ONEMULTADD(CFLOAT, float)
 
57
MAKE_C_ONEMULTADD(CDOUBLE, double)
 
58
#endif /* __GNUC__ */
 
59
 
 
60
static OneMultAddFunction *OneMultAdd[]={NULL,
 
61
                                         UCHAR_onemultadd,
 
62
                                         SCHAR_onemultadd,
 
63
                                         SHORT_onemultadd,
 
64
#ifdef PyArray_UNSIGNED_TYPES
 
65
                                         USHORT_onemultadd,
 
66
#endif
 
67
                                         INT_onemultadd,
 
68
#ifdef PyArray_UNSIGNED_TYPES
 
69
                                         UINT_onemultadd,
 
70
#endif
 
71
                                         LONG_onemultadd,
 
72
                                         FLOAT_onemultadd,
 
73
                                         DOUBLE_onemultadd,
 
74
                                         CFLOAT_onemultadd,
 
75
                                         CDOUBLE_onemultadd,
 
76
                                         NULL};
 
77
 
 
78
 
 
79
/* This could definitely be more optimized... */
 
80
 
 
81
int pylab_convolve_2d (char  *in,        /* Input data Ns[0] x Ns[1] */
 
82
                       int   *instr,     /* Input strides */
 
83
                       char  *out,       /* Output data */
 
84
                       int   *outstr,    /* Ouput strides */
 
85
                       char  *hvals,     /* coefficients in filter */
 
86
                       int   *hstr,      /* coefficients strides */ 
 
87
                       int   *Nwin,      /* Size of kernel Nwin[0] x Nwin[1] */
 
88
                       int   *Ns,        /* Size of image Ns[0] x Ns[1] */
 
89
                       int   flag,       /* convolution parameters */
 
90
                       char  *fillvalue) /* fill value */
 
91
{
 
92
  int bounds_pad_flag = 0;
 
93
  int m, n, j, k, ind0, ind1;
 
94
  int Os[2];
 
95
  char *sum=NULL, *value=NULL;
 
96
  int new_m, new_n, ind0_memory=0;
 
97
  int boundary, outsize, convolve, type_num, type_size;
 
98
  OneMultAddFunction *mult_and_add;
 
99
 
 
100
  boundary = flag & BOUNDARY_MASK;  /* flag can be fill, reflecting, circular */
 
101
  outsize = flag & OUTSIZE_MASK;
 
102
  convolve = flag & FLIP_MASK;
 
103
  type_num = (flag & TYPE_MASK) >> TYPE_SHIFT;
 
104
  /*type_size*/
 
105
 
 
106
  mult_and_add = OneMultAdd[type_num];
 
107
  if (mult_and_add == NULL) return -5;  /* Not available for this type */
 
108
 
 
109
  if (type_num < 0 || type_num > MAXTYPES) return -4;  /* Invalid type */
 
110
  type_size = elsizes[type_num];
 
111
 
 
112
  if ((sum = calloc(type_size,2))==NULL) return -3; /* No memory */
 
113
  value = sum + type_size;
 
114
 
 
115
  if (outsize == FULL) {Os[0] = Ns[0]+Nwin[0]-1; Os[1] = Ns[0]+Nwin[1]-1;}
 
116
  else if (outsize == SAME) {Os[0] = Ns[0]; Os[1] = Ns[1];}
 
117
  else if (outsize == VALID) {Os[0] = Ns[0]-Nwin[0]+1; Os[1] = Ns[1]-Nwin[1]+1;}
 
118
  else return -1;  /* Invalid output flag */  
 
119
  
 
120
  if ((boundary != PAD) && (boundary != REFLECT) && (boundary != CIRCULAR)) 
 
121
    return -2;   /* Invalid boundary flag */
 
122
 
 
123
  for (m=0; m < Os[0]; m++) {
 
124
    /* Reposition index into input image based on requested output size */
 
125
    if (outsize == FULL) new_m = convolve ? m : (m-Nwin[0]+1);
 
126
    else if (outsize == SAME) new_m = convolve ? (m+((Nwin[0]-1)>>1)) : (m-((Nwin[0]-1) >> 1));
 
127
    else new_m = convolve ? (m+Nwin[0]-1) : m; /* VALID */
 
128
 
 
129
    for (n=0; n < Os[1]; n++) {  /* loop over columns */
 
130
      memset(sum, 0, type_size); /* sum = 0.0; */
 
131
 
 
132
      if (outsize == FULL) new_n = convolve ? n : (n-Nwin[1]+1);
 
133
      else if (outsize == SAME) new_n = convolve ? (n+((Nwin[1]-1)>>1)) : (n-((Nwin[1]-1) >> 1));
 
134
      else new_n = convolve ? (n+Nwin[1]-1) : n;
 
135
 
 
136
      /* Sum over kernel, if index into image is out of bounds
 
137
         handle it according to boundary flag */
 
138
      for (j=0; j < Nwin[0]; j++) {
 
139
        ind0 = convolve ? (new_m-j): (new_m+j);
 
140
        bounds_pad_flag = 0;
 
141
 
 
142
        if (ind0 < 0) {
 
143
          if (boundary == REFLECT) ind0 = -1-ind0;
 
144
          else if (boundary == CIRCULAR) ind0 = Ns[0] + ind0;
 
145
          else bounds_pad_flag = 1;
 
146
        }
 
147
        else if (ind0 >= Ns[0]) {
 
148
          if (boundary == REFLECT) ind0 = Ns[0]+Ns[0]-1-ind0;
 
149
          else if (boundary == CIRCULAR) ind0 = ind0 - Ns[0];
 
150
          else bounds_pad_flag = 1;
 
151
        }
 
152
        
 
153
        if (!bounds_pad_flag) ind0_memory = ind0*instr[0];
 
154
 
 
155
        for (k=0; k < Nwin[1]; k++) {
 
156
          if (bounds_pad_flag) memcpy(value,fillvalue,type_size);
 
157
          else {
 
158
            ind1 = convolve ? (new_n-k) : (new_n+k);
 
159
            if (ind1 < 0) {
 
160
              if (boundary == REFLECT) ind1 = -1-ind1;
 
161
              else if (boundary == CIRCULAR) ind1 = Ns[1] + ind1;
 
162
              else bounds_pad_flag = 1;
 
163
            }
 
164
            else if (ind1 >= Ns[1]) {
 
165
              if (boundary == REFLECT) ind1 = Ns[1]+Ns[1]-1-ind1;
 
166
              else if (boundary == CIRCULAR) ind1 = ind1 - Ns[1];
 
167
              else bounds_pad_flag = 1;
 
168
            }
 
169
           
 
170
            if (bounds_pad_flag) memcpy(value, fillvalue, type_size);
 
171
            else memcpy(value, in+ind0_memory+ind1*instr[1], type_size);
 
172
            bounds_pad_flag = 0;
 
173
          }
 
174
          mult_and_add(sum, hvals+j*hstr[0]+k*hstr[1], value);
 
175
        }
 
176
        memcpy(out+m*outstr[0]+n*outstr[1], sum, type_size);
 
177
      }
 
178
    }
 
179
  }
 
180
  free(sum);
 
181
  return 0;
 
182
}
 
183
 
 
184
 
 
185