1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19
20 //! \addtogroup op_fft
21 //! @{
22
23
24
25 //
26 // op_fft_real
27
28
29
30 template<typename T1>
31 inline
32 void
apply(Mat<std::complex<typename T1::pod_type>> & out,const mtOp<std::complex<typename T1::pod_type>,T1,op_fft_real> & in)33 op_fft_real::apply( Mat< std::complex<typename T1::pod_type> >& out, const mtOp<std::complex<typename T1::pod_type>,T1,op_fft_real>& in )
34 {
35 arma_extra_debug_sigprint();
36
37 typedef typename T1::pod_type in_eT;
38 typedef typename std::complex<in_eT> out_eT;
39
40 const Proxy<T1> P(in.m);
41
42 const uword n_rows = P.get_n_rows();
43 const uword n_cols = P.get_n_cols();
44 const uword n_elem = P.get_n_elem();
45
46 const bool is_vec = ( (n_rows == 1) || (n_cols == 1) );
47
48 const uword N_orig = (is_vec) ? n_elem : n_rows;
49 const uword N_user = (in.aux_uword_b == 0) ? in.aux_uword_a : N_orig;
50
51 fft_engine<out_eT,false> worker(N_user);
52
53 // no need to worry about aliasing, as we're going from a real object to complex complex, which by definition cannot alias
54
55 if(is_vec)
56 {
57 (n_cols == 1) ? out.set_size(N_user, 1) : out.set_size(1, N_user);
58
59 if( (out.n_elem == 0) || (N_orig == 0) )
60 {
61 out.zeros();
62 return;
63 }
64
65 if( (N_user == 1) && (N_orig >= 1) )
66 {
67 out[0] = out_eT( P[0] );
68 return;
69 }
70
71 podarray<out_eT> data(N_user);
72
73 out_eT* data_mem = data.memptr();
74
75 if(N_user > N_orig) { arrayops::fill_zeros( &data_mem[N_orig], (N_user - N_orig) ); }
76
77 const uword N = (std::min)(N_user, N_orig);
78
79 if(Proxy<T1>::use_at == false)
80 {
81 typename Proxy<T1>::ea_type X = P.get_ea();
82
83 for(uword i=0; i < N; ++i) { data_mem[i] = out_eT( X[i], in_eT(0) ); }
84 }
85 else
86 {
87 if(n_cols == 1)
88 {
89 for(uword i=0; i < N; ++i) { data_mem[i] = out_eT( P.at(i,0), in_eT(0) ); }
90 }
91 else
92 {
93 for(uword i=0; i < N; ++i) { data_mem[i] = out_eT( P.at(0,i), in_eT(0) ); }
94 }
95 }
96
97 worker.run( out.memptr(), data_mem );
98 }
99 else
100 {
101 // process each column seperately
102
103 out.set_size(N_user, n_cols);
104
105 if( (out.n_elem == 0) || (N_orig == 0) )
106 {
107 out.zeros();
108 return;
109 }
110
111 if( (N_user == 1) && (N_orig >= 1) )
112 {
113 for(uword col=0; col < n_cols; ++col) { out.at(0,col) = out_eT( P.at(0,col) ); }
114
115 return;
116 }
117
118 podarray<out_eT> data(N_user);
119
120 out_eT* data_mem = data.memptr();
121
122 if(N_user > N_orig) { arrayops::fill_zeros( &data_mem[N_orig], (N_user - N_orig) ); }
123
124 const uword N = (std::min)(N_user, N_orig);
125
126 for(uword col=0; col < n_cols; ++col)
127 {
128 for(uword i=0; i < N; ++i) { data_mem[i] = P.at(i, col); }
129
130 worker.run( out.colptr(col), data_mem );
131 }
132 }
133 }
134
135
136
137 //
138 // op_fft_cx
139
140
141 template<typename T1>
142 inline
143 void
apply(Mat<typename T1::elem_type> & out,const Op<T1,op_fft_cx> & in)144 op_fft_cx::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_fft_cx>& in)
145 {
146 arma_extra_debug_sigprint();
147
148 typedef typename T1::elem_type eT;
149
150 const Proxy<T1> P(in.m);
151
152 if(P.is_alias(out) == false)
153 {
154 op_fft_cx::apply_noalias<T1,false>(out, P, in.aux_uword_a, in.aux_uword_b);
155 }
156 else
157 {
158 Mat<eT> tmp;
159
160 op_fft_cx::apply_noalias<T1,false>(tmp, P, in.aux_uword_a, in.aux_uword_b);
161
162 out.steal_mem(tmp);
163 }
164 }
165
166
167
168 template<typename T1, bool inverse>
169 inline
170 void
apply_noalias(Mat<typename T1::elem_type> & out,const Proxy<T1> & P,const uword a,const uword b)171 op_fft_cx::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword a, const uword b)
172 {
173 arma_extra_debug_sigprint();
174
175 typedef typename T1::elem_type eT;
176
177 const uword n_rows = P.get_n_rows();
178 const uword n_cols = P.get_n_cols();
179 const uword n_elem = P.get_n_elem();
180
181 const bool is_vec = ( (n_rows == 1) || (n_cols == 1) );
182
183 const uword N_orig = (is_vec) ? n_elem : n_rows;
184 const uword N_user = (b == 0) ? a : N_orig;
185
186 fft_engine<eT,inverse> worker(N_user);
187
188 if(is_vec)
189 {
190 (n_cols == 1) ? out.set_size(N_user, 1) : out.set_size(1, N_user);
191
192 if( (out.n_elem == 0) || (N_orig == 0) )
193 {
194 out.zeros();
195 return;
196 }
197
198 if( (N_user == 1) && (N_orig >= 1) )
199 {
200 out[0] = P[0];
201 return;
202 }
203
204 if( (N_user > N_orig) || (is_Mat<typename Proxy<T1>::stored_type>::value == false) )
205 {
206 podarray<eT> data(N_user);
207
208 eT* data_mem = data.memptr();
209
210 if(N_user > N_orig) { arrayops::fill_zeros( &data_mem[N_orig], (N_user - N_orig) ); }
211
212 op_fft_cx::copy_vec( data_mem, P, (std::min)(N_user, N_orig) );
213
214 worker.run( out.memptr(), data_mem );
215 }
216 else
217 {
218 const unwrap< typename Proxy<T1>::stored_type > tmp(P.Q);
219
220 worker.run( out.memptr(), tmp.M.memptr() );
221 }
222 }
223 else
224 {
225 // process each column seperately
226
227 out.set_size(N_user, n_cols);
228
229 if( (out.n_elem == 0) || (N_orig == 0) )
230 {
231 out.zeros();
232 return;
233 }
234
235 if( (N_user == 1) && (N_orig >= 1) )
236 {
237 for(uword col=0; col < n_cols; ++col) { out.at(0,col) = P.at(0,col); }
238
239 return;
240 }
241
242 if( (N_user > N_orig) || (is_Mat<typename Proxy<T1>::stored_type>::value == false) )
243 {
244 podarray<eT> data(N_user);
245
246 eT* data_mem = data.memptr();
247
248 if(N_user > N_orig) { arrayops::fill_zeros( &data_mem[N_orig], (N_user - N_orig) ); }
249
250 const uword N = (std::min)(N_user, N_orig);
251
252 for(uword col=0; col < n_cols; ++col)
253 {
254 for(uword i=0; i < N; ++i) { data_mem[i] = P.at(i, col); }
255
256 worker.run( out.colptr(col), data_mem );
257 }
258 }
259 else
260 {
261 const unwrap< typename Proxy<T1>::stored_type > tmp(P.Q);
262
263 for(uword col=0; col < n_cols; ++col)
264 {
265 worker.run( out.colptr(col), tmp.M.colptr(col) );
266 }
267 }
268 }
269
270
271 // correct the scaling for the inverse transform
272 if(inverse)
273 {
274 typedef typename get_pod_type<eT>::result T;
275
276 const T k = T(1) / T(N_user);
277
278 eT* out_mem = out.memptr();
279
280 const uword out_n_elem = out.n_elem;
281
282 for(uword i=0; i < out_n_elem; ++i) { out_mem[i] *= k; }
283 }
284 }
285
286
287
288 template<typename T1>
289 arma_hot
290 inline
291 void
copy_vec(typename Proxy<T1>::elem_type * dest,const Proxy<T1> & P,const uword N)292 op_fft_cx::copy_vec(typename Proxy<T1>::elem_type* dest, const Proxy<T1>& P, const uword N)
293 {
294 arma_extra_debug_sigprint();
295
296 if(is_Mat< typename Proxy<T1>::stored_type >::value)
297 {
298 op_fft_cx::copy_vec_unwrap(dest, P, N);
299 }
300 else
301 {
302 op_fft_cx::copy_vec_proxy(dest, P, N);
303 }
304 }
305
306
307
308 template<typename T1>
309 arma_hot
310 inline
311 void
copy_vec_unwrap(typename Proxy<T1>::elem_type * dest,const Proxy<T1> & P,const uword N)312 op_fft_cx::copy_vec_unwrap(typename Proxy<T1>::elem_type* dest, const Proxy<T1>& P, const uword N)
313 {
314 arma_extra_debug_sigprint();
315
316 const unwrap< typename Proxy<T1>::stored_type > tmp(P.Q);
317
318 arrayops::copy(dest, tmp.M.memptr(), N);
319 }
320
321
322
323 template<typename T1>
324 arma_hot
325 inline
326 void
copy_vec_proxy(typename Proxy<T1>::elem_type * dest,const Proxy<T1> & P,const uword N)327 op_fft_cx::copy_vec_proxy(typename Proxy<T1>::elem_type* dest, const Proxy<T1>& P, const uword N)
328 {
329 arma_extra_debug_sigprint();
330
331 if(Proxy<T1>::use_at == false)
332 {
333 typename Proxy<T1>::ea_type X = P.get_ea();
334
335 for(uword i=0; i < N; ++i) { dest[i] = X[i]; }
336 }
337 else
338 {
339 if(P.get_n_cols() == 1)
340 {
341 for(uword i=0; i < N; ++i) { dest[i] = P.at(i,0); }
342 }
343 else
344 {
345 for(uword i=0; i < N; ++i) { dest[i] = P.at(0,i); }
346 }
347 }
348 }
349
350
351
352 //
353 // op_ifft_cx
354
355
356 template<typename T1>
357 inline
358 void
apply(Mat<typename T1::elem_type> & out,const Op<T1,op_ifft_cx> & in)359 op_ifft_cx::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_ifft_cx>& in)
360 {
361 arma_extra_debug_sigprint();
362
363 typedef typename T1::elem_type eT;
364
365 const Proxy<T1> P(in.m);
366
367 if(P.is_alias(out) == false)
368 {
369 op_fft_cx::apply_noalias<T1,true>(out, P, in.aux_uword_a, in.aux_uword_b);
370 }
371 else
372 {
373 Mat<eT> tmp;
374
375 op_fft_cx::apply_noalias<T1,true>(tmp, P, in.aux_uword_a, in.aux_uword_b);
376
377 out.steal_mem(tmp);
378 }
379 }
380
381
382
383 //! @}
384