LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
FourierOps.h
Go to the documentation of this file.
1 // -*- c++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008, 2009, 2010, 2011 LSST Corporation.
5  *
6  * This product includes software developed by the
7  * LSST Project (http://www.lsst.org/).
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the LSST License Statement and
20  * the GNU General Public License along with this program. If not,
21  * see <http://www.lsstcorp.org/LegalNotices/>.
22  */
23 #ifndef NDARRAY_FFT_FourierOps_h_INCLUDED
24 #define NDARRAY_FFT_FourierOps_h_INCLUDED
25 
32 #include <boost/noncopyable.hpp>
33 
34 #include "ndarray.h"
35 
36 namespace ndarray {
38 namespace detail {
39 
44 template <typename T, int N>
45 struct FourierOps {
46 
47  template <int C>
48  static void shift(
49  T const * offset,
50  std::complex<T> const & factor,
51  ArrayRef<std::complex<T>,N,C> const & array,
52  int const real_last_dim
53  ) {
54  typename ArrayRef<std::complex<T>,N,C>::Iterator iter = array.begin();
55  T u = -2.0 * M_PI * (*offset) / array.size();
56  int kMid = (array.size() + 1) / 2;
57  for (int k = 0; k < kMid; ++k, ++iter) {
58  FourierOps<T,N-1>::shift(offset+1, factor * std::polar(static_cast<T>(1), u * k),
59  *iter, real_last_dim);
60  }
61  if (array.size() % 2 == 0) {
62  FourierOps<T,N-1>::shift(offset+1, factor * std::cos(u * kMid), *iter, real_last_dim);
63  ++iter;
64  ++kMid;
65  }
66  for (int k_n = kMid - array.size(); k_n < 0; ++k_n, ++iter) {
67  FourierOps<T,N-1>::shift(offset+1, factor * std::polar(static_cast<T>(1), u * k_n),
68  *iter, real_last_dim);
69  }
70  }
71 
72  template <int C>
73  static void differentiate(int m, ArrayRef<std::complex<T>,N,C> const & array, int const real_last_dim) {
74  typename ArrayRef<std::complex<T>,N,C>::Iterator iter = array.begin();
75  int kMid = (array.size() + 1) / 2;
76  T u = 2.0 * M_PI / array.size();
77  for (int k = 0; k < kMid; ++k, ++iter) {
78  if (m == N) (*iter) *= std::complex<T>(static_cast<T>(0), u * T(k));
79  FourierOps<T,N-1>::differentiate(m, *iter, real_last_dim);
80  }
81  if (array.size() % 2 == 0) {
82  (*iter) = static_cast<T>(0);
83  ++iter;
84  ++kMid;
85  }
86  for (int k_n = kMid - array.size(); k_n < 0; ++k_n, ++iter) {
87  if (m == N) (*iter) *= std::complex<T>(static_cast<T>(0), u * T(k_n));
88  FourierOps<T,N-1>::differentiate(m, *iter, real_last_dim);
89  }
90  }
91 
92 };
93 
98 template <typename T>
99 struct FourierOps<T,1> {
100 
101  template <int C>
102  static void shift(
103  T const * offset,
104  std::complex<T> const & factor,
105  ArrayRef<std::complex<T>,1,C> const & array,
106  int const real_last_dim
107  ) {
108  typename ArrayRef<std::complex<T>,1,C>::Iterator iter = array.begin();
109  T u = -2.0 * M_PI * (*offset) / real_last_dim;
110  int kMid = (real_last_dim + 1) / 2;
111  for (int k = 0; k < kMid; ++k, ++iter) {
112  (*iter) *= factor * std::polar(1.0, u * T(k));
113  }
114  if (real_last_dim % 2 == 0) {
115  (*iter) *= factor * std::cos(u * kMid);
116  ++iter;
117  }
118  }
119 
120  template <int C>
121  static void differentiate(int m, ArrayRef<std::complex<T>,1,C> const & array, int const real_last_dim) {
122  typename ArrayRef<std::complex<T>,1,C>::Iterator iter = array.begin();
123  int kMid = (real_last_dim + 1) / 2;
124  if (m == 1) {
125  T u = 2.0 * M_PI / real_last_dim;
126  for (int k = 0; k < kMid; ++k, ++iter) {
127  (*iter) *= std::complex<T>(static_cast<T>(0), u * T(k));
128  }
129  }
130  if (real_last_dim % 2 == 0) {
131  array[kMid] = static_cast<T>(0);
132  }
133  }
134 
135 };
136 
137 } // namespace detail
139 
145 template <typename T, int N, int C>
146 void shift(
147  Vector<T,N> const & offset,
148  Array<std::complex<T>,N,C> const & array,
149  int const real_last_dim
150 ) {
152  offset.begin(),
153  static_cast< std::complex<T> >(1),
154  array.deep(),
155  real_last_dim
156  );
157 }
158 
164 template <typename T, int N, int C>
166  int n,
167  Array<std::complex<T>,N,C> const & array,
168  int const real_last_dim
169 ) {
170  detail::FourierOps<T,N>::differentiate(N-n, array.deep(), real_last_dim);
171 }
172 
173 } // namespace ndarray
174 
175 #endif // !NDARRAY_FFT_FourierOps_h_INCLUDED
int iter
void shift(Vector< T, N > const &offset, Array< std::complex< T >, N, C > const &array, int const real_last_dim)
Perform a Fourier-space translation transform.
Definition: FourierOps.h:146
void differentiate(int n, Array< std::complex< T >, N, C > const &array, int const real_last_dim)
Numerically differentiate the array in Fourier-space in the given dimension.
Definition: FourierOps.h:165
A fixed-size 1D array class.
Definition: Vector.h:93
iterator begin()
Return an iterator to the beginning of the Vector.
Definition: Vector.h:120
A multidimensional strided array.
Definition: Array.h:47
tuple m
Definition: lsstimport.py:48