← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/fftpack/src/convolve.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/fftpack/src/convolve.c
 
1
/*
 
2
  Generic functions for computing 1D convolutions of periodic sequences.
 
3
 
 
4
  Supported FFT libraries:
 
5
    DJBFFT  - optional, used for power-of-two length arrays
 
6
    FFTW    - optional
 
7
    FFTPACK - used if any of the above libraries is not available
 
8
 
 
9
  Author: Pearu Peterson, September 2002
 
10
 */
 
11
 
 
12
#include "fftpack.h"
 
13
 
 
14
/**************** DJBFFT *****************************/
 
15
#ifdef WITH_DJBFFT
 
16
GEN_CACHE(ddjbfft,(int n)
 
17
          ,double* ptr;
 
18
          ,(caches_ddjbfft[i].n==n)
 
19
          ,caches_ddjbfft[id].ptr = (double*)malloc(sizeof(double)*n);
 
20
          ,free(caches_ddjbfft[id].ptr);
 
21
          ,20)
 
22
#endif
 
23
 
 
24
/**************** FFTW *****************************/
 
25
#ifdef WITH_FFTW
 
26
GEN_CACHE(drfftw,(int n)
 
27
          ,rfftw_plan plan1;
 
28
           rfftw_plan plan2;
 
29
          ,(caches_drfftw[i].n==n)
 
30
          ,caches_drfftw[id].plan1 = rfftw_create_plan(n,
 
31
                FFTW_REAL_TO_COMPLEX,
 
32
                FFTW_IN_PLACE|FFTW_ESTIMATE);
 
33
           caches_drfftw[id].plan2 = rfftw_create_plan(n,
 
34
                FFTW_COMPLEX_TO_REAL,
 
35
                FFTW_IN_PLACE|FFTW_ESTIMATE);
 
36
          ,rfftw_destroy_plan(caches_drfftw[id].plan1);
 
37
           rfftw_destroy_plan(caches_drfftw[id].plan2);
 
38
          ,20)
 
39
#else
 
40
/**************** FFTPACK ZFFT **********************/
 
41
extern void F_FUNC(dfftf,DFFTF)(int*,double*,double*);
 
42
extern void F_FUNC(dfftb,DFFTB)(int*,double*,double*);
 
43
extern void F_FUNC(dffti,DFFTI)(int*,double*);
 
44
GEN_CACHE(dfftpack,(int n)
 
45
          ,double* wsave;
 
46
          ,(caches_dfftpack[i].n==n)
 
47
          ,caches_dfftpack[id].wsave = (double*)malloc(sizeof(double)*(2*n+15));
 
48
           F_FUNC(dffti,DFFTI)(&n,caches_dfftpack[id].wsave);
 
49
          ,free(caches_dfftpack[id].wsave);
 
50
          ,20)
 
51
#endif
 
52
extern void destroy_convolve_cache(void) {
 
53
#ifdef WITH_DJBFFT
 
54
  destroy_ddjbfft_caches();
 
55
#endif
 
56
#ifdef WITH_FFTW
 
57
  destroy_drfftw_caches();
 
58
#else
 
59
  destroy_dfftpack_caches();
 
60
#endif
 
61
}
 
62
 
 
63
/**************** convolve **********************/
 
64
extern
 
65
void convolve(int n,double* inout,double* omega,int swap_real_imag) {
 
66
  int i;
 
67
#ifdef WITH_DJBFFT
 
68
  double* ptr = NULL;
 
69
#endif
 
70
#ifdef WITH_FFTW
 
71
  rfftw_plan plan1 = NULL;
 
72
  rfftw_plan plan2 = NULL;
 
73
#else
 
74
  double* wsave = NULL;
 
75
#endif
 
76
#ifdef WITH_DJBFFT
 
77
  switch (n) {
 
78
  case 2:;case 4:;case 8:;case 16:;case 32:;case 64:;case 128:;case 256:;
 
79
  case 512:;case 1024:;case 2048:;case 4096:;case 8192:
 
80
    i = get_cache_id_ddjbfft(n);
 
81
    ptr = caches_ddjbfft[i].ptr;
 
82
    COPYSTD2DJB(inout,ptr,n);
 
83
    switch (n) {
 
84
#define TMPCASE(N) case N: fftr8_##N(ptr); break
 
85
      TMPCASE(2);TMPCASE(4);TMPCASE(8);TMPCASE(16);TMPCASE(32);
 
86
      TMPCASE(64);TMPCASE(128);TMPCASE(256);TMPCASE(512);
 
87
      TMPCASE(1024);TMPCASE(2048);TMPCASE(4096);TMPCASE(8192);
 
88
#undef TMPCASE
 
89
    }
 
90
    if (swap_real_imag) {
 
91
      int n1 = n-1;
 
92
      double c;
 
93
      ptr[0] *= omega[0];
 
94
      ptr[1] *= omega[1];
 
95
      for(i=2;i<n1;i+=2) {
 
96
        c = ptr[i] * omega[i];
 
97
        ptr[i] = ptr[i+1] * omega[i+1];
 
98
        ptr[i+1] = c;
 
99
      }
 
100
    }
 
101
    else
 
102
      for(i=0;i<n;++i)
 
103
        ptr[i] *= omega[i];
 
104
    switch (n) {
 
105
#define TMPCASE(N)case N:fftr8_un##N(ptr);break
 
106
      TMPCASE(2);TMPCASE(4);TMPCASE(8);TMPCASE(16);TMPCASE(32);
 
107
      TMPCASE(64);TMPCASE(128);TMPCASE(256);TMPCASE(512);
 
108
      TMPCASE(1024);TMPCASE(2048);TMPCASE(4096);TMPCASE(8192);
 
109
#undef TMPCASE
 
110
    }
 
111
    COPYINVDJB2STD2(ptr,inout,n);
 
112
    return;
 
113
  }
 
114
#endif
 
115
  {
 
116
#ifdef WITH_FFTW
 
117
    int l = (n-1)/2+1;
 
118
    i = get_cache_id_drfftw(n);
 
119
    plan1 = caches_drfftw[i].plan1;
 
120
    plan2 = caches_drfftw[i].plan2;
 
121
    rfftw_one(plan1, (fftw_real *)inout, NULL);
 
122
    if (swap_real_imag) {
 
123
      double c;
 
124
      inout[0] *= omega[0];
 
125
      if (!(n%2))
 
126
        inout[n/2] *= omega[n/2];
 
127
      for(i=1;i<l;++i) {
 
128
        c = inout[i] * omega[i];
 
129
        inout[i] = omega[n-i] * inout[n-i];
 
130
        inout[n-i] = c;
 
131
      }
 
132
    }
 
133
    else
 
134
      for(i=0;i<n;++i)
 
135
        inout[i] *= omega[i];
 
136
    rfftw_one(plan2, (fftw_real *)inout, NULL);
 
137
#else
 
138
    i = get_cache_id_dfftpack(n);
 
139
    wsave = caches_dfftpack[i].wsave;
 
140
    F_FUNC(dfftf,DFFTF)(&n,inout,wsave);
 
141
    if (swap_real_imag) {
 
142
      double c;
 
143
      int n1 = n-1;
 
144
      inout[0] *= omega[0];
 
145
      if (!(n%2))
 
146
        inout[n-1] *= omega[n-1];
 
147
      for(i=1;i<n1;i+=2) {
 
148
        c = inout[i] * omega[i];
 
149
        inout[i] = inout[i+1] * omega[i+1];
 
150
        inout[i+1] = c;
 
151
      }
 
152
    }
 
153
    else
 
154
      for(i=0;i<n;++i)
 
155
        inout[i] *= omega[i];
 
156
    F_FUNC(dfftb,DFFTB)(&n,inout,wsave);
 
157
#endif
 
158
  }
 
159
}
 
160
 
 
161
/**************** convolve **********************/
 
162
extern
 
163
void convolve_z(int n,double* inout,double* omega_real,double* omega_imag) {
 
164
  int i;
 
165
#ifdef WITH_DJBFFT
 
166
  double* ptr = NULL;
 
167
#endif
 
168
#ifdef WITH_FFTW
 
169
  rfftw_plan plan1 = NULL;
 
170
  rfftw_plan plan2 = NULL;
 
171
#else
 
172
  double* wsave = NULL;
 
173
#endif
 
174
#ifdef WITH_DJBFFT
 
175
  switch (n) {
 
176
  case 2:;case 4:;case 8:;case 16:;case 32:;case 64:;case 128:;case 256:;
 
177
  case 512:;case 1024:;case 2048:;case 4096:;case 8192:
 
178
    i = get_cache_id_ddjbfft(n);
 
179
    ptr = caches_ddjbfft[i].ptr;
 
180
    COPYSTD2DJB(inout,ptr,n);
 
181
    switch (n) {
 
182
#define TMPCASE(N) case N: fftr8_##N(ptr); break
 
183
      TMPCASE(2);TMPCASE(4);TMPCASE(8);TMPCASE(16);TMPCASE(32);
 
184
      TMPCASE(64);TMPCASE(128);TMPCASE(256);TMPCASE(512);
 
185
      TMPCASE(1024);TMPCASE(2048);TMPCASE(4096);TMPCASE(8192);
 
186
#undef TMPCASE
 
187
    }
 
188
    {
 
189
      int n1 = n-1;
 
190
      double c;
 
191
      ptr[0] *= (omega_real[0]+omega_imag[0]);
 
192
      ptr[1] *= (omega_real[1]+omega_imag[1]);
 
193
      for(i=2;i<n1;i+=2) {
 
194
        c = ptr[i] * omega_imag[i];
 
195
        ptr[i] *= omega_real[i];
 
196
        ptr[i] += ptr[i+1] * omega_imag[i+1];
 
197
        ptr[i+1] *= omega_real[i+1];
 
198
        ptr[i+1] += c;
 
199
      }
 
200
    }
 
201
    switch (n) {
 
202
#define TMPCASE(N)case N:fftr8_un##N(ptr);break
 
203
      TMPCASE(2);TMPCASE(4);TMPCASE(8);TMPCASE(16);TMPCASE(32);
 
204
      TMPCASE(64);TMPCASE(128);TMPCASE(256);TMPCASE(512);
 
205
      TMPCASE(1024);TMPCASE(2048);TMPCASE(4096);TMPCASE(8192);
 
206
#undef TMPCASE
 
207
    }
 
208
    COPYINVDJB2STD2(ptr,inout,n);
 
209
    return;
 
210
  }
 
211
#endif
 
212
  {
 
213
#ifdef WITH_FFTW
 
214
    int l = (n-1)/2+1;
 
215
    i = get_cache_id_drfftw(n);
 
216
    plan1 = caches_drfftw[i].plan1;
 
217
    plan2 = caches_drfftw[i].plan2;
 
218
    rfftw_one(plan1, (fftw_real *)inout, NULL);
 
219
    {
 
220
      double c;
 
221
      inout[0] *= (omega_real[0]+omega_imag[0]);
 
222
      if (!(n%2))
 
223
        inout[n/2] *= (omega_real[n/2]+omega_imag[n/2]);
 
224
      for(i=1;i<l;++i) {
 
225
        c = inout[i] * omega_imag[i];
 
226
        inout[i] *= omega_real[i];
 
227
        inout[i] += omega_imag[n-i] * inout[n-i];
 
228
        inout[n-i] *= omega_real[n-i];
 
229
        inout[n-i] += c;
 
230
      }
 
231
    }
 
232
    rfftw_one(plan2, (fftw_real *)inout, NULL);
 
233
#else
 
234
    i = get_cache_id_dfftpack(n);
 
235
    wsave = caches_dfftpack[i].wsave;
 
236
    F_FUNC(dfftf,DFFTF)(&n,inout,wsave);
 
237
    {
 
238
      double c;
 
239
      int n1 = n-1;
 
240
      inout[0] *= (omega_real[0]+omega_imag[0]);
 
241
      if (!(n%2))
 
242
        inout[n-1] *= (omega_real[n-1]+omega_imag[n-1]);
 
243
      for(i=1;i<n1;i+=2) {
 
244
        c = inout[i] * omega_imag[i];
 
245
        inout[i] *= omega_real[i];
 
246
        inout[i] += inout[i+1] * omega_imag[i+1];
 
247
        inout[i+1] *= omega_real[i+1];
 
248
        inout[i+1] += c;
 
249
      }
 
250
    }
 
251
    F_FUNC(dfftb,DFFTB)(&n,inout,wsave);
 
252
#endif
 
253
  }
 
254
}
 
255
 
 
256
extern
 
257
void init_convolution_kernel(int n,double* omega, int d,
 
258
                             double (*kernel_func)(int),
 
259
                             int zero_nyquist) {
 
260
  /*
 
261
    omega[k] = pow(sqrt(-1),d) * kernel_func(k)
 
262
    omega[0] = kernel_func(0)
 
263
    conjugate(omega[-k]) == omega[k]
 
264
   */
 
265
#ifdef WITH_DJBFFT
 
266
  switch (n) {
 
267
  case 2:;case 4:;case 8:;case 16:;case 32:;case 64:;case 128:;case 256:;
 
268
  case 512:;case 1024:;case 2048:;case 4096:;case 8192:
 
269
    {
 
270
      int k,n2=n/2, *f = (int*)malloc(sizeof(int)*(n));
 
271
      fftfreq_rtable(f,n);
 
272
      for (k=1;k<n;++k)
 
273
        if (f[k]>n2) f[k] -= n;
 
274
      omega[0] = (*kernel_func)(0)/n;
 
275
      switch (d%4) {
 
276
      case 0:
 
277
        for (k=2;k<n-1;k+=2) {
 
278
          omega[k] = (*kernel_func)(f[k])/n2;
 
279
          omega[k+1] = -omega[k];
 
280
        }
 
281
        omega[1] = (zero_nyquist?0.0:(*kernel_func)(n2)/n);
 
282
        break;
 
283
      case 1:;case -3:
 
284
        for (k=2;k<n-1;k+=2)
 
285
          omega[k] = omega[k+1] = -(*kernel_func)(f[k])/n2;
 
286
        omega[1] = (zero_nyquist?0.0:(*kernel_func)(n2)/n);
 
287
        break;
 
288
      case 2:;case -2:
 
289
        for (k=2;k<n-1;k+=2) {
 
290
          omega[k] = -(*kernel_func)(f[k])/n2;
 
291
          omega[k+1] = -omega[k];
 
292
        }
 
293
        omega[1] = (zero_nyquist?0.0:-(*kernel_func)(n2)/n);
 
294
        break;
 
295
      case 3:;case -1:
 
296
        for (k=2;k<n-1;k+=2)
 
297
          omega[k] = omega[k+1] = (*kernel_func)(f[k])/n2;
 
298
        omega[1] = (zero_nyquist?0.0:-(*kernel_func)(n2)/n);
 
299
        break;
 
300
      }
 
301
      free(f);
 
302
    }
 
303
    return;
 
304
  }
 
305
#endif
 
306
#ifdef WITH_FFTW
 
307
  {
 
308
    int k,l=(n-1)/2+1;
 
309
    omega[0] = (*kernel_func)(0)/n;;
 
310
    switch (d%4) {
 
311
    case 0:
 
312
      for (k=1;k<l;++k)
 
313
        omega[k] = omega[n-k] = (*kernel_func)(k)/n;
 
314
      if (!(n%2)) 
 
315
        omega[n/2] = (zero_nyquist?0.0:(*kernel_func)(n/2)/n);
 
316
      break;
 
317
    case 1:;case -3:
 
318
      for (k=1;k<l;++k) {
 
319
        omega[k] = (*kernel_func)(k)/n;
 
320
        omega[n-k] = -omega[k];
 
321
      }
 
322
      if (!(n%2))
 
323
        omega[n/2] = (zero_nyquist?0.0:(*kernel_func)(n/2)/n);
 
324
      break;
 
325
    case 2:;case -2:
 
326
      for (k=1;k<l;++k)
 
327
        omega[k] = omega[n-k] = -(*kernel_func)(k)/n;
 
328
      if (!(n%2))
 
329
        omega[n/2] = (zero_nyquist?0.0:-(*kernel_func)(n/2)/n);
 
330
      break;
 
331
    case 3:;case -1:
 
332
      for (k=1;k<l;++k) {
 
333
        omega[k] = -(*kernel_func)(k)/n;
 
334
        omega[n-k] = -omega[k];
 
335
      }
 
336
      if (!(n%2))
 
337
        omega[n/2] = (zero_nyquist?0.0:-(*kernel_func)(n/2)/n);
 
338
      break;
 
339
    }
 
340
  }
 
341
#else
 
342
  {
 
343
    int j,k,l=(n%2?n:n-1);
 
344
    omega[0] = (*kernel_func)(0)/n;
 
345
    switch (d%4) {
 
346
    case 0:
 
347
      for (k=j=1;j<l;j+=2,++k)
 
348
        omega[j] = omega[j+1] = (*kernel_func)(k)/n;
 
349
      if (!(n%2))
 
350
        omega[n-1] = (zero_nyquist?0.0:(*kernel_func)(k)/n);
 
351
      break;
 
352
    case 1:;case -3:
 
353
      for (k=j=1;j<l;j+=2,++k) {
 
354
        omega[j] = (*kernel_func)(k)/n;
 
355
        omega[j+1] = -omega[j];
 
356
      }
 
357
      if (!(n%2))
 
358
        omega[n-1] = (zero_nyquist?0.0:(*kernel_func)(k)/n);
 
359
      break;
 
360
    case 2:;case -2:
 
361
      for (k=j=1;j<l;j+=2,++k)
 
362
        omega[j] = omega[j+1] = -(*kernel_func)(k)/n;
 
363
      if (!(n%2))
 
364
        omega[n-1] = (zero_nyquist?0.0:-(*kernel_func)(k)/n);
 
365
      break;
 
366
    case 3:;case -1:
 
367
      for (k=j=1;j<l;j+=2,++k) {
 
368
        omega[j] = -(*kernel_func)(k)/n;
 
369
        omega[j+1] = -omega[j];
 
370
      }
 
371
      if (!(n%2))
 
372
        omega[n-1] = (zero_nyquist?0.0:-(*kernel_func)(k)/n);
 
373
      break;
 
374
    }
 
375
  }
 
376
#endif
 
377
}