fftwaff.cc 10.3 KB
Newer Older
thomas.forbriger's avatar
thomas.forbriger committed
1
2
3
4
5
/*! \file fftwaff.cc
 * \brief use fftw together with aff containers (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
thomas.forbriger's avatar
thomas.forbriger committed
6
 * $Id$
thomas.forbriger's avatar
thomas.forbriger committed
7
8
9
10
11
12
13
14
15
 * \author Thomas Forbriger
 * \date 11/07/2006
 * 
 * use fftw together with aff containers (implementation)
 * 
 * Copyright (c) 2006 by Thomas Forbriger (BFO Schiltach) 
 * 
 * REVISIONS and CHANGES 
 *  - 11/07/2006   V1.0   Thomas Forbriger
thomas.forbriger's avatar
thomas.forbriger committed
16
 *  - 12/09/2007   V1.1   first running version
thomas.forbriger's avatar
thomas.forbriger committed
17
18
 *  - 07/10/2010   V1.2   
 *                        - migrate to FFTW3
thomas.forbriger's avatar
thomas.forbriger committed
19
20
 *  - 13/05/2011   V1.3 
 *                        - properly delete plans in delete_arrays() function
thomas.forbriger's avatar
thomas.forbriger committed
21
 *
thomas.forbriger's avatar
thomas.forbriger committed
22
23
24
 * ============================================================================
 */
#define TF_FFTWAFF_CC_VERSION \
thomas.forbriger's avatar
thomas.forbriger committed
25
  "TF_FFTWAFF_CC   V1.3"
thomas.forbriger's avatar
thomas.forbriger committed
26
#define TF_FFTWAFF_CC_CVSID \
thomas.forbriger's avatar
thomas.forbriger committed
27
  "$Id$"
thomas.forbriger's avatar
thomas.forbriger committed
28
29

#include <fourier/fftwaff.h>
thomas.forbriger's avatar
thomas.forbriger committed
30
#include <fourier/error.h>
thomas.forbriger's avatar
thomas.forbriger committed
31
#include <fourier/error.h>
thomas.forbriger's avatar
thomas.forbriger committed
32
#include <aff/seriesoperators.h>
thomas.forbriger's avatar
thomas.forbriger committed
33
#include <aff/dump.h>
thomas.forbriger's avatar
thomas.forbriger committed
34
#include <cmath>
thomas.forbriger's avatar
thomas.forbriger committed
35
36
37
38
39
40
41

namespace fourier {

  /*! All Fourier transform stuff is collected here.
   */
  namespace fft {

thomas.forbriger's avatar
thomas.forbriger committed
42
    //! create plan.
thomas.forbriger's avatar
thomas.forbriger committed
43
44
45
46
    void DRFFTWAFF::create_plan_forward() const
    {
      if (Mplan_forward==0)
      {
thomas.forbriger's avatar
thomas.forbriger committed
47
        //Mplan_forward=rfftw_create_plan(Msize, FFTW_FORWARD, 0);
thomas.forbriger's avatar
thomas.forbriger committed
48
#ifdef FFTWFALLBACK
thomas.forbriger's avatar
thomas.forbriger committed
49
50
        Mplan_forward=rfftw_create_plan(Msize, FFTW_REAL_TO_COMPLEX,
                                        FFTW_ESTIMATE);
thomas.forbriger's avatar
thomas.forbriger committed
51
52
53
54
55
56
57
#else
        this->create_arrays();
        Mplan_forward=fftw_plan_dft_r2c_1d(Msize, Mseriesarray,
                                           Mspectrumarray,
                                           FFTW_ESTIMATE);
#endif
        FOURIER_assert(Mplan_forward!=0, 
thomas.forbriger's avatar
thomas.forbriger committed
58
59
60
61
62
63
64
                    "Error (DRFFTWAFF::create_plan_forward): "
                    "could not create plan!")
      }
    } // void DRFFTWAFF::create_plan_forward() const

    /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
65
    //! create plan.
thomas.forbriger's avatar
thomas.forbriger committed
66
67
68
69
    void DRFFTWAFF::create_plan_backward() const
    {
      if (Mplan_backward==0)
      {
thomas.forbriger's avatar
thomas.forbriger committed
70
#ifdef FFTWFALLBACK
thomas.forbriger's avatar
thomas.forbriger committed
71
        Mplan_backward=rfftw_create_plan(Msize, FFTW_BACKWARD, 0);
thomas.forbriger's avatar
thomas.forbriger committed
72
73
74
75
76
77
78
#else
        this->create_arrays();
        Mplan_backward=fftw_plan_dft_c2r_1d(Msize, Mspectrumarray,
                                           Mseriesarray,
                                           FFTW_ESTIMATE);
#endif
        FOURIER_assert(Mplan_backward!=0, 
thomas.forbriger's avatar
thomas.forbriger committed
79
80
81
82
83
84
85
                    "Error (DRFFTWAFF::create_plan_backward): "
                    "could not create plan!")
      }
    } // void DRFFTWAFF::create_plan_backward() const

    /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
86
    //! delete plan.
thomas.forbriger's avatar
thomas.forbriger committed
87
    DRFFTWAFF::~DRFFTWAFF()
thomas.forbriger's avatar
thomas.forbriger committed
88
89
    {
      this->delete_plans();
thomas.forbriger's avatar
thomas.forbriger committed
90
91
92
#ifndef FFTWFALLBACK
      fftw_cleanup();
#endif
thomas.forbriger's avatar
thomas.forbriger committed
93
94
95
96
97
98
99
100
101
102
103
    } // DRFFTWAFF::~DRFFTWAFF()

    /*----------------------------------------------------------------------*/

    //! prepare FFT settings for size n.
    void DRFFTWAFF::set_size(const int& n) const
    {
      if (n != Msize) 
      {
        Msize=n;
        this->delete_plans();
104
#ifndef FFTWFALLBACK
thomas.forbriger's avatar
thomas.forbriger committed
105
106
        this->delete_arrays();
#endif
thomas.forbriger's avatar
thomas.forbriger committed
107
108
109
110
111
112
113
      }
    } // DRFFTWAFF::~DRFFTWAFF()

    /*----------------------------------------------------------------------*/

    //! delete plans.
    void DRFFTWAFF::delete_plans() const
thomas.forbriger's avatar
thomas.forbriger committed
114
    {
thomas.forbriger's avatar
thomas.forbriger committed
115
#ifdef FFTWFALLBACK
thomas.forbriger's avatar
thomas.forbriger committed
116
117
      if (Mplan_forward != 0) { rfftw_destroy_plan(Mplan_forward); }
      if (Mplan_backward != 0) { rfftw_destroy_plan(Mplan_backward); }
thomas.forbriger's avatar
thomas.forbriger committed
118
119
120
121
#else
      if (Mplan_forward != 0) { fftw_destroy_plan(Mplan_forward); }
      if (Mplan_backward != 0) { fftw_destroy_plan(Mplan_backward); }
#endif
thomas.forbriger's avatar
thomas.forbriger committed
122
123
      Mplan_forward=0;
      Mplan_backward=0;
thomas.forbriger's avatar
thomas.forbriger committed
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
    } // DRFFTWAFF::delete_plans()

    /*----------------------------------------------------------------------*/

#ifndef FFTWFALLBACK
    //! create plans.
    void DRFFTWAFF::create_arrays() const
    {
      this->delete_arrays();
      Mseriesarray = (double *) fftw_malloc(sizeof(double)*Msize);
      FOURIER_assert(Mseriesarray!=0, 
                  "Error (DRFFTWAFF::create_plan_forward): "
                  "could not create series array!")
      Mseries=Tseries(Tseries::Trepresentation(Mseriesarray, Msize));
      Mspectrumarray = (fftw_complex *)
thomas.forbriger's avatar
thomas.forbriger committed
139
        fftw_malloc(sizeof(fftw_complex)*this->ssize());
thomas.forbriger's avatar
thomas.forbriger committed
140
141
142
143
144
145
146
147
148
149
150
151
152
      FOURIER_assert(Mspectrumarray!=0, 
                  "Error (DRFFTWAFF::create_plan_forward): "
                  "could not create spectrum array!")
      Mspectrum=Tspectrum(Tspectrum::Trepresentation(
                                  reinterpret_cast<Tcoeff*>(Mspectrumarray),
                                                     this->ssize()));
    } // DRFFTWAFF::create_arrays()

    /*----------------------------------------------------------------------*/

    //! delete arrays.
    void DRFFTWAFF::delete_arrays() const
    {
thomas.forbriger's avatar
thomas.forbriger committed
153
      this->delete_plans();
thomas.forbriger's avatar
thomas.forbriger committed
154
155
156
157
      if (Mseriesarray != 0) { fftw_free(Mseriesarray); }
      if (Mspectrumarray != 0) { fftw_free(Mspectrumarray); }
      Mseries=Tseries(0);
      Mspectrum=Tspectrum(0);
thomas.forbriger's avatar
thomas.forbriger committed
158
159
      // Mplan_forward=0;
      // Mplan_backward=0;
thomas.forbriger's avatar
thomas.forbriger committed
160
161
    } // DRFFTWAFF::delete_arrays()
#endif
thomas.forbriger's avatar
thomas.forbriger committed
162
163
164

    /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
165
166
167
168
    /*! \brief Transform time series to Fourier coefficients.
     *
     * No scaling is applied.
     */
thomas.forbriger's avatar
thomas.forbriger committed
169
170
    DRFFTWAFF::Tspectrum DRFFTWAFF::operator()(const Tseries::Tcoc& s,
                                               const bool& debug) const
thomas.forbriger's avatar
thomas.forbriger committed
171
    {
thomas.forbriger's avatar
thomas.forbriger committed
172
      this->set_size(s.size());
thomas.forbriger's avatar
thomas.forbriger committed
173
#ifdef FFTWFALLBACK
174
175
      FOURIER_debug(debug, "DRFFTWAFF::operator()",
                 "use fallback code (FFTW2)");
thomas.forbriger's avatar
thomas.forbriger committed
176
177
178
179
180
      Tspectrum retval(this->Msize/2+1);
      aff::Series<fftw_real> out(this->Msize);
      aff::Series<fftw_real> in(this->Msize);
      fftw_real* pout=out.pointer();
      fftw_real* pin=in.pointer();
thomas.forbriger's avatar
thomas.forbriger committed
181
      FOURIER_debug(debug, "DRFFTWAFF::operator()",
thomas.forbriger's avatar
thomas.forbriger committed
182
                 "processing arrays are created; copy in series");
thomas.forbriger's avatar
thomas.forbriger committed
183
      in.copyin(s);
thomas.forbriger's avatar
thomas.forbriger committed
184
      
thomas.forbriger's avatar
thomas.forbriger committed
185
      FOURIER_debug(debug, "DRFFTWAFF::operator()",
thomas.forbriger's avatar
thomas.forbriger committed
186
                 "create plan forward");
thomas.forbriger's avatar
thomas.forbriger committed
187
      this->create_plan_forward();
thomas.forbriger's avatar
thomas.forbriger committed
188
      FOURIER_debug(debug, "DRFFTWAFF::operator()",
thomas.forbriger's avatar
thomas.forbriger committed
189
                 "call rfftw_one");
thomas.forbriger's avatar
thomas.forbriger committed
190
      rfftw_one(Mplan_forward, pin, pout);
thomas.forbriger's avatar
thomas.forbriger committed
191
      FOURIER_debug(debug, "DRFFTWAFF::operator()",
thomas.forbriger's avatar
thomas.forbriger committed
192
                 "copy results to output");
thomas.forbriger's avatar
thomas.forbriger committed
193
194
195
196
197
198
199
200
201
      retval(0)=out(0);
      for (int i=1; i<((Msize+1)/2); ++i)
      {
        retval(i)=Tcoeff(out(i),out(Msize-i));
      }
      if ((Msize % 2) == 0)
      {
        retval(Msize/2)=out(Msize/2);
      }
thomas.forbriger's avatar
thomas.forbriger committed
202
      /*
thomas.forbriger's avatar
thomas.forbriger committed
203
204
205
206
      if (debug)
      {
        DUMP(in);
        DUMP(out);
thomas.forbriger's avatar
thomas.forbriger committed
207
208
209
      }
      */
#else
210
211
      FOURIER_debug(debug, "DRFFTWAFF::operator()",
                 "use recent code (FFTW3)");
thomas.forbriger's avatar
thomas.forbriger committed
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
      Tspectrum retval(this->ssize());
      this->create_plan_forward();
      Mseries.copyin(s);
      /*
      for (int i=0; i<this->size(); ++i)
      { Mseries(Mseries.first()+i)=s(s.first()+i); }
      */
      fftw_execute(Mplan_forward);
      retval.copyin(Mspectrum);
      /*
      for (int i=0; <this->ssize(); ++i)
      { retval(retval.first()+i)=Mspectrum(Mspectrum.first()+i); }
      */
#endif
      /*
      if (debug)
      {
        DUMP(s);
thomas.forbriger's avatar
thomas.forbriger committed
230
231
232
        DUMP(retval);
      }
      */
thomas.forbriger's avatar
thomas.forbriger committed
233
      FOURIER_debug(debug, "DRFFTWAFF::operator()",
thomas.forbriger's avatar
thomas.forbriger committed
234
                 "finished; return");
thomas.forbriger's avatar
thomas.forbriger committed
235
236
237
238
239
      return(retval);
    } // Tspectrum DRFFTWAFF::operator()(const Tseries::Tcoc& s) const

    /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
240
241
242
243
    /*! \brief Transform Fourier coefficients to time series.
     *
     * No scaling is applied.
     */
thomas.forbriger's avatar
thomas.forbriger committed
244
245
    DRFFTWAFF::Tseries DRFFTWAFF::operator()(const Tspectrum::Tcoc& s,
                                             const bool& debug) const
thomas.forbriger's avatar
thomas.forbriger committed
246
    {
thomas.forbriger's avatar
thomas.forbriger committed
247
248
249
250
251
252
253
      // adjust FFT size
      int seriessize=s.size()*2-1;
      // is Nyquits coefficients real?
      if (std::abs(s(s.size()).imag()) < 1.e-8*std::abs(s(s.size()).real()))
      { --seriessize; }
      this->set_size(seriessize);

thomas.forbriger's avatar
thomas.forbriger committed
254
#ifdef FFTWFALLBACK
thomas.forbriger's avatar
thomas.forbriger committed
255
      Tseries retval(Msize);
thomas.forbriger's avatar
thomas.forbriger committed
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
      aff::Series<fftw_real> out(Msize);
      aff::Series<fftw_real> in(Msize);
      fftw_real* pout=out.pointer();
      fftw_real* pin=in.pointer();
      in(0)=s(0).real();
      for (int i=1; i<((Msize+1)/2); ++i)
      {
        in(i)=s(i).real();
        in(Msize-i)=s(i).imag();
      }
      if ((Msize % 2) == 0)
      {
        in(Msize/2)=s(Msize/2).real();
      }
      this->create_plan_backward();
      rfftw_one(Mplan_backward, pin, pout);
      retval.copyin(out);
thomas.forbriger's avatar
thomas.forbriger committed
273
274
275
276
277
278
279
#else
      Tseries retval(this->size());
      this->create_plan_backward();
      Mspectrum.copyin(s);
      fftw_execute(Mplan_backward);
      retval.copyin(Mseries);
#endif
thomas.forbriger's avatar
thomas.forbriger committed
280
281
282
283
284
      return(retval);
    } // Tseries DRFFTWAFF::operator()(const Tspectrum::Tcoc& s) const

    /*----------------------------------------------------------------------*/
    
thomas.forbriger's avatar
thomas.forbriger committed
285
    //! Return appropriate scaling factor for sampling interval dt.
thomas.forbriger's avatar
thomas.forbriger committed
286
287
288
289
290
291
292
    DRFFTWAFF::Tsample DRFFTWAFF::scale_series(const Tsample& dt) const
    {
      return(1./(Msize*dt));
    } // Tsample DRFFTWAFF::scale_series(const Tsample& dt) const

    /*----------------------------------------------------------------------*/
    
thomas.forbriger's avatar
thomas.forbriger committed
293
    //! Return appropriate scaling factor for sampling interval dt.
thomas.forbriger's avatar
thomas.forbriger committed
294
295
296
297
    DRFFTWAFF::Tsample DRFFTWAFF::scale_spectrum(const Tsample& dt) const
    {
      return(dt);
    } // Tsample DRFFTWAFF::scale_spectrum(const Tsample& dt) const
thomas.forbriger's avatar
thomas.forbriger committed
298

thomas.forbriger's avatar
thomas.forbriger committed
299
300
    /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
301
302
303
304
305
    /*! \brief Transform time series to Fourier coefficients and scale.
     *
     * Appropriate scaling is applied for sampling interval dt.
     */
    DRFFTWAFF::Tspectrum DRFFTWAFF::operator()(const Tseries::Tcoc& s,
thomas.forbriger's avatar
thomas.forbriger committed
306
307
                            const double& dt,
                            const bool& debug) const
thomas.forbriger's avatar
thomas.forbriger committed
308
    {
thomas.forbriger's avatar
thomas.forbriger committed
309
      Tspectrum retval=this->operator()(s, debug);
thomas.forbriger's avatar
thomas.forbriger committed
310
311
312
313
314
315
      retval *= this->scale_spectrum(dt);
      return(retval);
    }

    /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
316
317
318
319
320
    /*! \brief Transform Fourier coefficients to time series and scale.
     *
     * Appropriate scaling is applied for sampling interval dt.
     */
    DRFFTWAFF::Tseries DRFFTWAFF::operator()(const Tspectrum::Tcoc& s,
thomas.forbriger's avatar
thomas.forbriger committed
321
322
                            const double& dt,
                            const bool& debug) const
thomas.forbriger's avatar
thomas.forbriger committed
323
324
    {
      return(this->scale_series(dt)*
thomas.forbriger's avatar
thomas.forbriger committed
325
             this->operator()(s, debug));
thomas.forbriger's avatar
thomas.forbriger committed
326
327
328
329
330
331
332
    }

  } // namespace ftt

} // namespace fourier

/* ----- END OF fftwaff.cc ----- */