comparison armadillo-3.900.4/include/armadillo_bits/op_fft_meat.hpp @ 49:1ec0e2823891

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