Replace Min, Max and Abs by std::min, std::max and std::abs
[scilab.git] / scilab / modules / signal_processing / sci_gateway / cpp / sci_fft.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Antoine ELIAS
4 *
5 * This file must be used under the terms of the CeCILL.
6 * This source file is licensed as described in the file COPYING, which
7 * you should have received as part of this distribution.  The terms
8 * are also available at
9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12 /*--------------------------------------------------------------------------*/
13
14 #include "signal_gw.hxx"
15 #include "double.hxx"
16 #include "function.hxx"
17
18 extern "C"
19 {
20 #include "sci_malloc.h"
21 #include "localization.h"
22 #include "Scierror.h"
23 #include "elem_common.h"
24
25     extern void C2F(fft842)(int *inverse, int *signal_length, double *signal_real, double *signal_imaginary, int *error);
26     extern void C2F(dfft2)(double *signal_real, double *signal_imaginary, int *nseg, int *n, int *nspn, int *inverse, int *error, int *buffer_data, int *buffer_size);
27 }
28
29 #define POW2_15     32768
30
31 bool isPowerOf2(int _iVal);
32
33 int maxfactor(int n);
34 int fft_1dim(double *signal_real, double *signal_imaginary, int signal_length, int inverse, int *buffer_data, int buffer_size);
35 int fft_2dim(double *signal_real, double *signal_imaginary, int signal_rows, int signal_cols, int inverse, int *buffer_data, int buffer_size);
36 int fft_ndim(double *signal_real, double *signal_imaginary, int signal_length, int dimensions_length, int dimension_stride, int inverse, int *buffer_data, int buffer_size);
37
38 /*--------------------------------------------------------------------------*/
39 types::Function::ReturnValue sci_fft(types::typed_list &in, int _iRetCount, types::typed_list &out)
40 {
41     int iDimLength          = 0;
42     int iDimCount           = 0;
43     int iInc                = 0;
44     int iWay                = -1;
45     int iSize               = 0;
46     int iOne                = 1;
47     int iErr                = 0;
48     double dblZero          = 0;
49
50     //workspace
51     int iWS                 = 0;
52     int* piWS               = NULL;
53     types::Double* pIn1     = NULL;
54
55     //check input parameters
56     if (in.size() != 1 && in.size() != 2 && in.size() != 4)
57     {
58         Scierror(77, _("%s: Wrong number of input argument(s): %d or %d expected.\n"), "fft", 1, 4);
59         return types::Function::Error;
60     }
61
62     switch (in.size())
63     {
64         case 4 :
65             //check fourth input parameter : inc
66             if (in[3]->isDouble() == false || in[3]->getAs<types::Double>()->isScalar() == false)
67             {
68                 Scierror(999, _("%s: Wrong type for input argument #%d: Scalar expected.\n"), "fft", 4);
69                 return types::Function::Error;
70             }
71
72             iInc = (int)in[3]->getAs<types::Double>()->get(0);
73
74             //check third input parameter : dim
75             if (in[2]->isDouble() == false || in[2]->getAs<types::Double>()->isScalar() == false)
76             {
77                 Scierror(999, _("%s: Wrong type for input argument #%d: Scalar expected.\n"), "fft", 3);
78                 return types::Function::Error;
79             }
80
81             iDimLength = (int)in[2]->getAs<types::Double>()->get(0);
82             iDimCount = 3; //any value > 2 (used as a flag)
83
84         case 2 :
85             //check third input parameter : way
86             if (in[1]->isDouble() == false || in[1]->getAs<types::Double>()->isScalar() == false)
87             {
88                 Scierror(999, _("%s: Wrong type for input argument #%d: Scalar expected.\n"), "fft", 2);
89                 return types::Function::Error;
90             }
91
92             iWay = (int)in[1]->getAs<types::Double>()->get(0);
93             if (iWay != -1 && iWay != 1)
94             {
95                 Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"), "fft", 2, "-1 1");
96                 return types::Function::Error;
97             }
98
99         case 1:
100             if (in[0]->isDouble() == false)
101             {
102                 Scierror(999, _("%s: Wrong type for input argument #%d: Scalar expected.\n"), "fft", 1);
103                 return types::Function::Error;
104             }
105
106             pIn1 = in[0]->getAs<types::Double>();
107
108             iDimCount = std::max(iDimCount, ((pIn1->getRows() == 1 || pIn1->getCols() == 1) ? 1 : 2));
109             iSize = pIn1->getSize();
110             break;
111         default :
112         {
113             Scierror(77, _("%s: Wrong number of input argument(s): %d eor %d xpected.\n"), "fft", 1, 4);
114             return types::Function::Error;
115         }
116     }
117
118     types::Double* pOut = pIn1->clone()->getAs<types::Double>();
119     pOut->setComplex(true);
120
121     //alloc workspace required by dfft2
122     iWS = 8 * maxfactor(iSize) + 24;
123     piWS = (int*)MALLOC(iWS * sizeof(int));
124     if (piWS == NULL)
125     {
126         Scierror(999, _("%s : Memory allocation error.\n"), "fft");
127         return types::Function::Error;
128     }
129
130     switch (iDimCount)
131     {
132         case 1 :
133             iErr = fft_1dim(pOut->getReal(), pOut->getImg(), iSize, iWay, piWS, iWS);
134             break;
135         case 2 :
136             iErr = fft_2dim(pOut->getReal(), pOut->getImg(), pOut->getRows(), pOut->getCols(), iWay, piWS, iWS);
137             if (iErr == 1)
138             {
139                 Scierror(999, _("%s : Memory allocation error.\n"), "fft");
140                 return types::Function::Error;
141             }
142             break;
143         default :
144             iErr = fft_ndim(pOut->getReal(), pOut->getImg(), iSize, iDimLength, iInc, iWay, piWS, iWS);
145             break;
146     }
147
148     FREE(piWS);
149
150     out.push_back(pOut);
151     return types::Function::OK;
152 }
153
154 int fft_1dim(double *signal_real, double *signal_imaginary, int signal_length, int inverse, int *buffer_data, int buffer_size)
155 {
156     int error = 0;
157     int one = 1;
158
159     if (isPowerOf2(signal_length) && signal_length < POW2_15)
160     {
161         C2F(fft842)(&inverse, &signal_length, signal_real, signal_imaginary, &error);
162     }
163     else
164     {
165         C2F(dfft2)(signal_real, signal_imaginary, &one, &signal_length, &one, &inverse, &error, buffer_data, &buffer_size);
166     }
167
168     return error;
169 }
170
171 int fft_2dim(double *signal_real, double *signal_imaginary, int signal_rows, int signal_cols, int inverse, int *buffer_data, int buffer_size)
172 {
173     int error = 0;
174     int one = 1;
175
176     if (isPowerOf2(signal_rows) && signal_rows < POW2_15)
177     {
178         for (int i = 0 ; i < signal_cols ; i++)
179         {
180             C2F(fft842)(&inverse, &signal_rows, &(signal_real[signal_rows * i]), &(signal_imaginary[signal_rows * i]), &error);
181         }
182     }
183     else
184     {
185         C2F(dfft2)(signal_real, signal_imaginary, &signal_cols, &signal_rows, &one, &inverse, &error, buffer_data, &buffer_size);
186     }
187
188     if (isPowerOf2(signal_cols) && signal_cols < POW2_15)
189     {
190         double* temp_real = NULL;
191         double* temp_imaginary = NULL;
192
193         temp_real = (double *)MALLOC(signal_cols * sizeof(double));
194         temp_imaginary = (double *)MALLOC(signal_cols * sizeof(double));
195
196         if (temp_real == NULL || temp_imaginary == NULL)
197         {
198             return 1;
199         }
200
201         for (int i = 0 ; i < signal_rows ; i++)
202         {
203             C2F(dcopy)(&signal_cols, &(signal_real[i]), &signal_rows, temp_real, &one);
204             C2F(dcopy)(&signal_cols, &(signal_imaginary[i]), &signal_rows, temp_imaginary, &one);
205             C2F(fft842)(&inverse, &signal_cols, temp_real, temp_imaginary, &error);
206             C2F(dcopy)(&signal_cols, temp_real, &one, &(signal_real[i]), &signal_rows);
207             C2F(dcopy)(&signal_cols, temp_imaginary, &one, &(signal_imaginary[i]), &signal_rows);
208         }
209
210         FREE(temp_imaginary);
211         temp_imaginary = NULL;
212         free(temp_real);
213         temp_real = NULL;
214     }
215     else
216     {
217         /* erroneous implementation suspected */
218         C2F(dfft2)(signal_real, signal_imaginary, &one, &signal_cols, &signal_rows, &inverse, &error, buffer_data, &buffer_size);
219     }
220
221     return error;
222 }
223
224 int fft_ndim(double *signal_real, double *signal_imaginary, int signal_length, int dimensions_length, int dimension_stride, int inverse, int *buffer_data, int buffer_size)
225 {
226     /* translated litterally from Fortran... but... wtf ?! */
227     int error = 0;
228     int nseg = signal_length / dimensions_length / dimension_stride;
229     C2F(dfft2)(signal_real, signal_imaginary, &nseg, &dimensions_length, &dimension_stride, &inverse, &error, buffer_data, &buffer_size);
230     return error;
231 }
232
233
234 int maxfactor(int n)
235 {
236     int nfac[15];
237     int m = 0, j = 0, jj = 0, k = 0, kt = 0, max = 0;
238
239     for (k = abs(n) ; k % 16 == 0 ; k /= 16)
240     {
241         m++;
242         nfac[m - 1] = 4;
243     }
244
245     for (j = 3, jj = 9 ; jj <= k ; j += 2, jj = j * j)
246     {
247         if (k % jj != 0)
248         {
249             continue;
250         }
251         m++;
252         nfac[m - 1] = j;
253         k /= jj;
254     }
255
256     if (k <= 4)
257     {
258         kt = m;
259         nfac[m] = k;
260         if (k != 1)
261         {
262             m++;
263         }
264     }
265     else
266     {
267         if (k % 4 == 0)
268         {
269             m++;
270             nfac[m - 1] = 2;
271             k /= 4;
272         }
273
274         kt = m;
275         for (j = 2 ; j <= k ; j = ((j + 1) / 2) * 2 + 1)
276         {
277             if (k % j != 0)
278             {
279                 continue;
280             }
281
282             m++;
283             nfac[m - 1] = j;
284             k /= j;
285         }
286     }
287
288     if (kt != 0)
289     {
290         for (j = kt ; j > 0 ; j--)
291         {
292             m++;
293             nfac[m - 1] = nfac[j - 1];
294         }
295     }
296
297     /* get nfac maximum */
298     for (j = 0, max = nfac[0]; j < m ; j++)
299     {
300         max = (nfac[j] > max) ? nfac[j] : max;
301     }
302
303     return max;
304 }
305
306 bool isPowerOf2(int _iVal)
307 {
308     return (_iVal & (_iVal - 1)) == 0;
309 }