fftwaffar.cc 7.34 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
40
/*! \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$"

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

namespace fourier {

  namespace fft {

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

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

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

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

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

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

    //! delete plans.
77
    void DRFFTWAFFArrayEngine::delete_plans() 
thomas.forbriger's avatar
thomas.forbriger committed
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
    {
      if (Mplanr2c != 0) { fftw_destroy_plan(Mplanr2c); }
      if (Mplanc2r != 0) { fftw_destroy_plan(Mplanc2r); }
      Mplanr2c=0;
      Mplanc2r=0;
    } // DRFFTWAFFArrayEngine::delete_plans()

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

    //! execute r2c plan
    void DRFFTWAFFArrayEngine::r2c() 
    {
      if (Mplanr2c == 0)
      {

      }

    } // DRFFTWAFFArrayEngine::r2c()

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

99
    void DRFFTWAFFArrayEngine::checkconsistency()
thomas.forbriger's avatar
thomas.forbriger committed
100
    {
101
      FOURIER_assert(Mseriesarray.size(0)==Mspectrumarray.size(0),
thomas.forbriger's avatar
thomas.forbriger committed
102
                     "ERROR: inconsistent number of signals");
103
104
105
      FOURIER_assert(Mseriesarray.size(1)==
                     DRFFTWAFFArrayEngine::ncoeff(Mspectrumarray.size(1)),
                     "ERROR: sample size inconsistent");
106
      FOURIER_assert(((Mseriesarray.size(2)==1) && (Mspectrumarray.size(2)==1)),
thomas.forbriger's avatar
thomas.forbriger committed
107
                     "ERROR: two-dimensional arrays only");
108
      FOURIER_assert(((Mseriesarray.size(3)==1) && (Mspectrumarray.size(3)==1)),
thomas.forbriger's avatar
thomas.forbriger committed
109
                     "ERROR: two-dimensional arrays only");
110
    } // void DRFFTWAFFArrayEngine::checkconsistency()
thomas.forbriger's avatar
thomas.forbriger committed
111
112
113

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

114
    void DRFFTWAFFArrayEngine::createplanr2c()
thomas.forbriger's avatar
thomas.forbriger committed
115
116
117
    {
      if (Mplanr2c==0)
      {
118
119
        // overall parameters
        // ------------------
thomas.forbriger's avatar
thomas.forbriger committed
120
121
        // tranformation only along one dimension
        const int rank=1;
122
123
        // size of each series
        int n=Mseriesarray.size(0);
thomas.forbriger's avatar
thomas.forbriger committed
124
        // number of signals to be transformed
125
        const int howmany=Mseriesarray.size(1);
thomas.forbriger's avatar
thomas.forbriger committed
126
127
        // quick design
        const unsigned flags=FFTW_ESTIMATE;
128
129
130

        // input array
        // -----------
131
        aff::CArray<TAseries::Tvalue> Cseriesarray(Mseriesarray);
132
133
134
        // one-dimensional transform: use default
        int* inembed=0;
        // distance to next sample
135
        const int istride=Cseriesarray.stride(0);
136
        // distance to next signal
137
        const int idist=Cseriesarray.stride(1);
138
        // casted pointer
139
        double* in=Cseriesarray.castedpointer<double>();
140
141
142

        // output array
        // ------------
143
        aff::CArray<TAspectrum::Tvalue> Cspectrumarray(Mspectrumarray);
144
145
146
        // one-dimensional transform: use default
        int* onembed=0;
        // distance to next sample
147
        const int ostride=Cspectrumarray.stride(0);
148
        // distance to next signal
149
        const int odist=Cspectrumarray.stride(1);
150
        // casted pointer
151
        fftw_complex* out=Cspectrumarray.castedpointer<fftw_complex>();
152
153
154
155
156
157

        // 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
158
159
        FOURIER_assert(Mplanr2c!=0, "ERROR: creating r2c plan");
      }
160
    } // void DRFFTWAFFArrayEngine::createplanr2c()
thomas.forbriger's avatar
thomas.forbriger committed
161
162
163

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

164
    void DRFFTWAFFArrayEngine::createplanc2r()
thomas.forbriger's avatar
thomas.forbriger committed
165
166
167
    {
      if (Mplanc2r==0)
      {
168
169
170
171
172
173
174
175
176
177
178
179
180
        // 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
        // -----------
181
        aff::CArray<TAspectrum::Tvalue> Cspectrumarray(Mspectrumarray);
182
183
184
        // one-dimensional transform: use default
        int* inembed=0;
        // distance to next sample
185
        const int istride=Cspectrumarray.stride(0);
186
        // distance to next signal
187
        const int idist=Cspectrumarray.stride(1);
188
        // casted pointer
189
        fftw_complex* in=Cspectrumarray.castedpointer<fftw_complex>();
190
191
192

        // output array
        // ------------
193
        aff::CArray<TAseries::Tvalue> Cseriesarray(Mseriesarray);
194
195
196
        // one-dimensional transform: use default
        int* onembed=0;
        // distance to next sample
197
        const int ostride=Cseriesarray.stride(0);
198
        // distance to next signal
199
        const int odist=Cseriesarray.stride(1);
200
        // casted pointer
201
        double* out=Cseriesarray.castedpointer<double>();
202
203
204
205
206
207

        // 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
208
209
        FOURIER_assert(Mplanr2c!=0, "ERROR: creating c2r plan");
      }
210
    } // void DRFFTWAFFArrayEngine::createplanc2r()
thomas.forbriger's avatar
thomas.forbriger committed
211
212
213
214
215
216

  } // namespace fft

} // namespace fourier

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