← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/fftpack/src/zfft.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/fftpack/src/zfft.c
 
1
/*
 
2
  Interface to various FFT libraries.
 
3
  Double complex FFT and IFFT.
 
4
  Author: Pearu Peterson, August 2002
 
5
 */
 
6
 
 
7
#include "fftpack.h"
 
8
 
 
9
/**************** FFTWORK *****************************/
 
10
 
 
11
#ifdef WITH_FFTWORK
 
12
GEN_CACHE(zfftwork,(int n)
 
13
          ,coef_dbl* coef;
 
14
          ,caches_zfftwork[i].n==n
 
15
          ,caches_zfftwork[id].coef = (coef_dbl*)malloc(sizeof(coef_dbl)*(n));
 
16
           fft_coef_dbl(caches_zfftwork[id].coef,n);
 
17
          ,free(caches_zfftwork[id].coef);
 
18
          ,10)
 
19
#endif
 
20
 
 
21
/**************** DJBFFT *****************************/
 
22
#ifdef WITH_DJBFFT
 
23
GEN_CACHE(zdjbfft,(int n)
 
24
          ,unsigned int* f;
 
25
           double* ptr;
 
26
          ,caches_zdjbfft[i].n==n
 
27
          ,caches_zdjbfft[id].f = (unsigned int*)malloc(sizeof(unsigned int)*(n));
 
28
           caches_zdjbfft[id].ptr = (double*)malloc(sizeof(double)*(2*n));
 
29
           fftfreq_ctable(caches_zdjbfft[id].f,n);
 
30
           for(i=0;i<n;++i) 
 
31
             caches_zdjbfft[id].f[i] = (n-caches_zdjbfft[id].f[i])%n;
 
32
          ,free(caches_zdjbfft[id].f);
 
33
           free(caches_zdjbfft[id].ptr);
 
34
          ,10)
 
35
#endif
 
36
 
 
37
/**************** FFTW *****************************/
 
38
#ifdef WITH_FFTW
 
39
GEN_CACHE(zfftw,(int n,int d)
 
40
          ,int direction;
 
41
           fftw_plan plan;
 
42
          ,((caches_zfftw[i].n==n) && 
 
43
            (caches_zfftw[i].direction==d))
 
44
          ,caches_zfftw[id].direction = d;
 
45
           caches_zfftw[id].plan = fftw_create_plan(n,
 
46
                (d>0?FFTW_FORWARD:FFTW_BACKWARD),
 
47
                FFTW_IN_PLACE|FFTW_ESTIMATE);
 
48
          ,fftw_destroy_plan(caches_zfftw[id].plan);
 
49
          ,10)
 
50
#else
 
51
/**************** FFTPACK ZFFT **********************/
 
52
extern void F_FUNC(zfftf,ZFFTF)(int*,double*,double*);
 
53
extern void F_FUNC(zfftb,ZFFTB)(int*,double*,double*);
 
54
extern void F_FUNC(zffti,ZFFTI)(int*,double*);
 
55
GEN_CACHE(zfftpack,(int n)
 
56
          ,double* wsave;
 
57
          ,(caches_zfftpack[i].n==n)
 
58
          ,caches_zfftpack[id].wsave = (double*)malloc(sizeof(double)*(4*n+15));
 
59
           F_FUNC(zffti,ZFFTI)(&n,caches_zfftpack[id].wsave);
 
60
          ,free(caches_zfftpack[id].wsave);
 
61
          ,10)
 
62
#endif
 
63
 
 
64
extern void destroy_zfft_cache(void) {
 
65
#ifdef WITH_FFTWORK
 
66
  destroy_zfftwork_caches();
 
67
#endif
 
68
#ifdef WITH_DJBFFT
 
69
  destroy_zdjbfft_caches();
 
70
#endif
 
71
#ifdef WITH_FFTW
 
72
  destroy_zfftw_caches();
 
73
#else
 
74
  destroy_zfftpack_caches();
 
75
#endif
 
76
}
 
77
 
 
78
/**************** ZFFT function **********************/
 
79
extern void zfft(complex_double *inout,
 
80
                 int n,int direction,int howmany,int normalize) {
 
81
  int i;
 
82
  complex_double *ptr = inout;
 
83
#ifdef WITH_FFTW
 
84
  fftw_plan plan = NULL;
 
85
#else
 
86
  double* wsave = NULL;
 
87
#endif
 
88
#ifdef WITH_FFTWORK
 
89
  coef_dbl* coef = NULL;
 
90
#endif
 
91
#ifdef WITH_DJBFFT
 
92
  int j;
 
93
  complex_double *ptrc = NULL;
 
94
  unsigned int *f = NULL;
 
95
#endif
 
96
#ifdef WITH_FFTWORK
 
97
  if (ispow2le2e30(n)) {
 
98
    i = get_cache_id_zfftwork(n);
 
99
    coef = caches_zfftwork[i].coef;
 
100
  } else
 
101
#endif
 
102
#ifdef WITH_DJBFFT
 
103
  switch (n) {
 
104
  case 2:;case 4:;case 8:;case 16:;case 32:;case 64:;case 128:;case 256:;
 
105
  case 512:;case 1024:;case 2048:;case 4096:;case 8192:
 
106
    i = get_cache_id_zdjbfft(n);
 
107
    f = caches_zdjbfft[i].f;
 
108
    ptrc = (complex_double*)caches_zdjbfft[i].ptr;
 
109
  }
 
110
  if (f==0)
 
111
#endif
 
112
#ifdef WITH_FFTW
 
113
    plan = caches_zfftw[get_cache_id_zfftw(n,direction)].plan;
 
114
#else
 
115
    wsave = caches_zfftpack[get_cache_id_zfftpack(n)].wsave;
 
116
#endif
 
117
 
 
118
  switch (direction) {
 
119
 
 
120
  case 1:
 
121
    for (i=0;i<howmany;++i,ptr+=n) {
 
122
#ifdef WITH_FFTWORK
 
123
      if (coef!=NULL) {
 
124
        fft_for_cplx_flt((cplx_dbl*)ptr,coef,n);
 
125
      } else
 
126
#endif
 
127
#ifdef WITH_DJBFFT
 
128
      if (f!=NULL) {
 
129
        memcpy(ptrc,ptr,2*n*sizeof(double));
 
130
        switch (n) {
 
131
#define TMPCASE(N) case N:  fftc8_##N(ptrc); break
 
132
          TMPCASE(2);TMPCASE(4);TMPCASE(8);TMPCASE(16);TMPCASE(32);
 
133
          TMPCASE(64);TMPCASE(128);TMPCASE(256);TMPCASE(512);
 
134
          TMPCASE(1024);TMPCASE(2048);TMPCASE(4096);TMPCASE(8192);
 
135
#undef TMPCASE
 
136
        }
 
137
        for (j=0;j<n;++j) *(ptr+f[j]) = *(ptrc+j);
 
138
      } else
 
139
#endif
 
140
#ifdef WITH_FFTW
 
141
        fftw_one(plan,(fftw_complex*)ptr,NULL);
 
142
#else
 
143
        F_FUNC(zfftf,ZFFTF)(&n,(double*)(ptr),wsave);
 
144
#endif
 
145
    }
 
146
    break;
 
147
 
 
148
  case -1:
 
149
    for (i=0;i<howmany;++i,ptr+=n) {
 
150
#ifdef WITH_FFTWORK
 
151
      if (coef!=NULL) {
 
152
        fft_bak_cplx_flt((cplx_dbl*)ptr,coef,n);
 
153
      } else
 
154
#endif
 
155
#ifdef WITH_DJBFFT
 
156
      if (f!=NULL) {
 
157
        for (j=0;j<n;++j) *(ptrc+j) = *(ptr+f[j]);
 
158
        switch (n) {
 
159
#define TMPCASE(N) case N:  fftc8_un##N(ptrc); break
 
160
          TMPCASE(2);TMPCASE(4);TMPCASE(8);TMPCASE(16);TMPCASE(32);
 
161
          TMPCASE(64);TMPCASE(128);TMPCASE(256);TMPCASE(512);
 
162
          TMPCASE(1024);TMPCASE(2048);TMPCASE(4096);TMPCASE(8192);
 
163
#undef TMPCASE
 
164
        }
 
165
        memcpy(ptr,ptrc,2*n*sizeof(double));
 
166
      } else
 
167
#endif
 
168
#ifdef WITH_FFTW
 
169
        fftw_one(plan,(fftw_complex*)ptr,NULL);
 
170
#else
 
171
        F_FUNC(zfftb,ZFFTB)(&n,(double*)(ptr),wsave);
 
172
#endif
 
173
    }
 
174
    break;
 
175
 
 
176
  default:
 
177
    fprintf(stderr,"zfft: invalid direction=%d\n",direction);
 
178
  }
 
179
 
 
180
  if (normalize) {
 
181
    ptr = inout;
 
182
#ifdef WITH_DJBFFT
 
183
    if (f!=NULL) {
 
184
      for (i=0;i<howmany;++i,ptr+=n) {
 
185
        switch (n) {
 
186
#define TMPCASE(N) case N:  fftc8_scale##N(ptr); break
 
187
          TMPCASE(2);TMPCASE(4);TMPCASE(8);TMPCASE(16);TMPCASE(32);
 
188
          TMPCASE(64);TMPCASE(128);TMPCASE(256);TMPCASE(512);
 
189
          TMPCASE(1024);TMPCASE(2048);TMPCASE(4096);TMPCASE(8192);
 
190
#undef TMPCASE
 
191
        }
 
192
      }
 
193
    } else
 
194
#endif
 
195
      for (i=n*howmany-1;i>=0;--i) {
 
196
        *((double*)(ptr)) /= n;
 
197
        *((double*)(ptr++)+1) /= n;
 
198
      }
 
199
  }
 
200
}