fftwaffar.cc 8.63 KB
Newer Older
thomas.forbriger's avatar
thomas.forbriger committed
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
/*! \file fftwaffar.cc
 * \brief engine to transfrom several signals at once (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
 * $Id$
 * \author Thomas Forbriger
 * \date 13/05/2011
 * 
 * engine to transfrom several signals at once (implementation)
 * 
 * Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach) 
 *
 * ----
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version. 
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 * ----
 * 
 * REVISIONS and CHANGES 
 *  - 13/05/2011   V1.0   Thomas Forbriger
 * 
 * ============================================================================
 */
#define TF_FFTWAFFAR_CC_VERSION \
  "TF_FFTWAFFAR_CC   V1.0   "
#define TF_FFTWAFFAR_CC_CVSID \
  "$Id$"

40
41
#include <fourier/fftwaffar.h>
#include <fourier/error.h>
42
#include <aff/Carray.h>
thomas.forbriger's avatar
thomas.forbriger committed
43
44
45
46
47
48
49

namespace fourier {

  namespace fft {

    DRFFTWAFFArrayEngine::DRFFTWAFFArrayEngine(const int& nseis,
                                               const int& nsamp)
50
      : Mseriesarray(TAseries(nsamp,nseis)),
51
      Mspectrumarray(TAspectrum(DRFFTWAFFArrayEngine::ncoeff(nsamp),nseis)),
thomas.forbriger's avatar
thomas.forbriger committed
52
53
      Mplanr2c(0), Mplanc2r(0) 
    {
54
      this->checkconsistency();
thomas.forbriger's avatar
thomas.forbriger committed
55
56
57
58
59
60
    }

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

    DRFFTWAFFArrayEngine::DRFFTWAFFArrayEngine(const TAseries& series,
                                               const TAspectrum& spec)
61
      : Mseriesarray(series), Mspectrumarray(spec),
thomas.forbriger's avatar
thomas.forbriger committed
62
63
64
65
66
67
68
69
70
71
72
73
74
75
      Mplanr2c(0), Mplanc2r(0) 
    { this->checkconsistency(); }

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

    //! delete plan.
    DRFFTWAFFArrayEngine::~DRFFTWAFFArrayEngine()
    {
      this->delete_plans();
    } // DRFFTWAFFArrayEngine::~DRFFTWAFFArrayEngine()

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

    //! delete plans.
76
    void DRFFTWAFFArrayEngine::delete_plans() 
thomas.forbriger's avatar
thomas.forbriger committed
77
78
79
80
81
82
83
84
85
    {
      if (Mplanr2c != 0) { fftw_destroy_plan(Mplanr2c); }
      if (Mplanc2r != 0) { fftw_destroy_plan(Mplanc2r); }
      Mplanr2c=0;
      Mplanc2r=0;
    } // DRFFTWAFFArrayEngine::delete_plans()

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

86
87
88
89
90
91
92
93
94
95
96
97
    //! execute c2r plan
    void DRFFTWAFFArrayEngine::c2r() 
    {
      if (Mplanc2r == 0)
      {
        this->createplanc2r();
      }
      fftw_execute(Mplanc2r);
    } // DRFFTWAFFArrayEngine::c2r()

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

thomas.forbriger's avatar
thomas.forbriger committed
98
99
100
101
102
    //! execute r2c plan
    void DRFFTWAFFArrayEngine::r2c() 
    {
      if (Mplanr2c == 0)
      {
103
        this->createplanr2c();
thomas.forbriger's avatar
thomas.forbriger committed
104
      }
105
      fftw_execute(Mplanr2c);
thomas.forbriger's avatar
thomas.forbriger committed
106
107
108
109
    } // DRFFTWAFFArrayEngine::r2c()

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

110
    void DRFFTWAFFArrayEngine::checkconsistency()
thomas.forbriger's avatar
thomas.forbriger committed
111
    {
112
      FOURIER_assert(Mseriesarray.size(0)==Mspectrumarray.size(0),
thomas.forbriger's avatar
thomas.forbriger committed
113
                     "ERROR: inconsistent number of signals");
114
115
116
      FOURIER_assert(Mseriesarray.size(1)==
                     DRFFTWAFFArrayEngine::ncoeff(Mspectrumarray.size(1)),
                     "ERROR: sample size inconsistent");
117
      FOURIER_assert(((Mseriesarray.size(2)==1) && (Mspectrumarray.size(2)==1)),
thomas.forbriger's avatar
thomas.forbriger committed
118
                     "ERROR: two-dimensional arrays only");
119
      FOURIER_assert(((Mseriesarray.size(3)==1) && (Mspectrumarray.size(3)==1)),
thomas.forbriger's avatar
thomas.forbriger committed
120
                     "ERROR: two-dimensional arrays only");
121
    } // void DRFFTWAFFArrayEngine::checkconsistency()
thomas.forbriger's avatar
thomas.forbriger committed
122
123
124

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

125
    void DRFFTWAFFArrayEngine::createplanr2c()
thomas.forbriger's avatar
thomas.forbriger committed
126
127
128
    {
      if (Mplanr2c==0)
      {
129
130
        // overall parameters
        // ------------------
thomas.forbriger's avatar
thomas.forbriger committed
131
132
        // tranformation only along one dimension
        const int rank=1;
133
134
        // size of each series
        int n=Mseriesarray.size(0);
thomas.forbriger's avatar
thomas.forbriger committed
135
        // number of signals to be transformed
136
        const int howmany=Mseriesarray.size(1);
thomas.forbriger's avatar
thomas.forbriger committed
137
138
        // quick design
        const unsigned flags=FFTW_ESTIMATE;
139
140
141

        // input array
        // -----------
142
        aff::CArray<TAseries::Tvalue> Cseriesarray(Mseriesarray);
143
144
145
        // one-dimensional transform: use default
        int* inembed=0;
        // distance to next sample
146
        const int istride=Cseriesarray.stride(0);
147
        // distance to next signal
148
        const int idist=Cseriesarray.stride(1);
149
        // casted pointer
150
        double* in=Cseriesarray.castedpointer<double>();
151
152
153

        // output array
        // ------------
154
        aff::CArray<TAspectrum::Tvalue> Cspectrumarray(Mspectrumarray);
155
156
157
        // one-dimensional transform: use default
        int* onembed=0;
        // distance to next sample
158
        const int ostride=Cspectrumarray.stride(0);
159
        // distance to next signal
160
        const int odist=Cspectrumarray.stride(1);
161
        // casted pointer
162
        fftw_complex* out=Cspectrumarray.castedpointer<fftw_complex>();
163
164
165
166
167
168

        // create plan
        Mplanr2c=fftw_plan_many_dft_r2c(rank, &n, howmany,
                                        in, inembed, istride, idist,
                                        out, onembed, ostride, odist,
                                        flags);
thomas.forbriger's avatar
thomas.forbriger committed
169
170
        FOURIER_assert(Mplanr2c!=0, "ERROR: creating r2c plan");
      }
171
    } // void DRFFTWAFFArrayEngine::createplanr2c()
thomas.forbriger's avatar
thomas.forbriger committed
172
173
174

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

175
    void DRFFTWAFFArrayEngine::createplanc2r()
thomas.forbriger's avatar
thomas.forbriger committed
176
177
178
    {
      if (Mplanc2r==0)
      {
179
180
181
182
183
184
185
186
187
188
189
190
191
        // overall parameters
        // ------------------
        // tranformation only along one dimension
        const int rank=1;
        // size of each series
        int n=Mseriesarray.size(0);
        // number of signals to be transformed
        const int howmany=Mseriesarray.size(1);
        // quick design
        const unsigned flags=FFTW_ESTIMATE;

        // input array
        // -----------
192
        aff::CArray<TAspectrum::Tvalue> Cspectrumarray(Mspectrumarray);
193
194
195
        // one-dimensional transform: use default
        int* inembed=0;
        // distance to next sample
196
        const int istride=Cspectrumarray.stride(0);
197
        // distance to next signal
198
        const int idist=Cspectrumarray.stride(1);
199
        // casted pointer
200
        fftw_complex* in=Cspectrumarray.castedpointer<fftw_complex>();
201
202
203

        // output array
        // ------------
204
        aff::CArray<TAseries::Tvalue> Cseriesarray(Mseriesarray);
205
206
207
        // one-dimensional transform: use default
        int* onembed=0;
        // distance to next sample
208
        const int ostride=Cseriesarray.stride(0);
209
        // distance to next signal
210
        const int odist=Cseriesarray.stride(1);
211
        // casted pointer
212
        double* out=Cseriesarray.castedpointer<double>();
213
214
215
216
217
218

        // create plan
        Mplanr2c=fftw_plan_many_dft_c2r(rank, &n, howmany,
                                        in, inembed, istride, idist,
                                        out, onembed, ostride, odist,
                                        flags);
thomas.forbriger's avatar
thomas.forbriger committed
219
220
        FOURIER_assert(Mplanr2c!=0, "ERROR: creating c2r plan");
      }
221
    } // void DRFFTWAFFArrayEngine::createplanc2r()
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
    
    /*----------------------------------------------------------------------*/

    //! Return appropriate scaling factor for sampling interval dt.
    DRFFTWAFFArrayEngine::Tsample
      DRFFTWAFFArrayEngine::scale_series(const Tsample& dt) const
    {
      return(1./(Mseriesarray.size(0)*dt));
    } // Tsample DRFFTWAFFArrayEngine::scale_series(const Tsample& dt) const

    /*----------------------------------------------------------------------*/
    
    //! Return appropriate scaling factor for sampling interval dt.
    DRFFTWAFFArrayEngine::Tsample
      DRFFTWAFFArrayEngine::scale_spectrum(const Tsample& dt) const
    {
      return(dt);
    } // Tsample DRFFTWAFFArrayEngine::scale_spectrum(const Tsample& dt) const

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

    unsigned int DRFFTWAFFArrayEngine::nseries() const
    {
      return Mseriesarray.size(1);
    } // unsigned int DRFFTWAFFArrayEngine::nseries() const
thomas.forbriger's avatar
thomas.forbriger committed
247
248
249
250
251
252

  } // namespace fft

} // namespace fourier

/* ----- END OF fftwaffar.cc ----- */