utf: module signal_processing
[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 or %d expected.\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(iDimLength == 0 ? iSize : iDimLength) + 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     switch (iDimCount)
130     {
131         case 1 :
132             iErr = fft_1dim(pOut->getReal(), pOut->getImg(), iSize, iWay, piWS, iWS);
133             break;
134         case 2 :
135             iErr = fft_2dim(pOut->getReal(), pOut->getImg(), pOut->getRows(), pOut->getCols(), iWay, piWS, iWS);
136             if (iErr == 1)
137             {
138                 Scierror(999, _("%s : Memory allocation error.\n"), "fft");
139                 return types::Function::Error;
140             }
141             break;
142         default :
143             iErr = fft_ndim(pOut->getReal(), pOut->getImg(), iSize, iDimLength, iInc, iWay, piWS, iWS);
144             break;
145     }
146     double *df = pOut->getImg();
147     bool cplx = false;
148     for (int i = 0; i < iSize; i++)
149     {
150         if (df[i] != 0)
151         {
152             cplx = true;
153             break;
154         }
155     }
156     if (cplx == false)
157     {
158         pOut->setComplex(false);
159     }
160
161     FREE(piWS);
162
163     out.push_back(pOut);
164     return types::Function::OK;
165 }
166
167 int fft_1dim(double *signal_real, double *signal_imaginary, int signal_length, int inverse, int *buffer_data, int buffer_size)
168 {
169     int error = 0;
170     int one = 1;
171
172     if (isPowerOf2(signal_length) && signal_length < POW2_15)
173     {
174         C2F(fft842)(&inverse, &signal_length, signal_real, signal_imaginary, &error);
175     }
176     else
177     {
178         C2F(dfft2)(signal_real, signal_imaginary, &one, &signal_length, &one, &inverse, &error, buffer_data, &buffer_size);
179     }
180
181     return error;
182 }
183
184 int fft_2dim(double *signal_real, double *signal_imaginary, int signal_rows, int signal_cols, int inverse, int *buffer_data, int buffer_size)
185 {
186     int error = 0;
187     int one = 1;
188
189     if (isPowerOf2(signal_rows) && signal_rows < POW2_15)
190     {
191         for (int i = 0 ; i < signal_cols ; i++)
192         {
193             C2F(fft842)(&inverse, &signal_rows, &(signal_real[signal_rows * i]), &(signal_imaginary[signal_rows * i]), &error);
194         }
195     }
196     else
197     {
198         C2F(dfft2)(signal_real, signal_imaginary, &signal_cols, &signal_rows, &one, &inverse, &error, buffer_data, &buffer_size);
199     }
200
201     if (isPowerOf2(signal_cols) && signal_cols < POW2_15)
202     {
203         double* temp_real = NULL;
204         double* temp_imaginary = NULL;
205
206         temp_real = (double *)MALLOC(signal_cols * sizeof(double));
207         temp_imaginary = (double *)MALLOC(signal_cols * sizeof(double));
208
209         if (temp_real == NULL || temp_imaginary == NULL)
210         {
211             return 1;
212         }
213
214         for (int i = 0 ; i < signal_rows ; i++)
215         {
216             C2F(dcopy)(&signal_cols, &(signal_real[i]), &signal_rows, temp_real, &one);
217             C2F(dcopy)(&signal_cols, &(signal_imaginary[i]), &signal_rows, temp_imaginary, &one);
218             C2F(fft842)(&inverse, &signal_cols, temp_real, temp_imaginary, &error);
219             C2F(dcopy)(&signal_cols, temp_real, &one, &(signal_real[i]), &signal_rows);
220             C2F(dcopy)(&signal_cols, temp_imaginary, &one, &(signal_imaginary[i]), &signal_rows);
221         }
222
223         FREE(temp_imaginary);
224         temp_imaginary = NULL;
225         free(temp_real);
226         temp_real = NULL;
227     }
228     else
229     {
230         /* erroneous implementation suspected */
231         C2F(dfft2)(signal_real, signal_imaginary, &one, &signal_cols, &signal_rows, &inverse, &error, buffer_data, &buffer_size);
232     }
233
234     return error;
235 }
236
237 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)
238 {
239     /* translated litterally from Fortran... but... wtf ?! */
240     int error = 0;
241     int nseg = signal_length / dimensions_length / dimension_stride;
242     C2F(dfft2)(signal_real, signal_imaginary, &nseg, &dimensions_length, &dimension_stride, &inverse, &error, buffer_data, &buffer_size);
243
244     return error;
245 }
246
247
248 int maxfactor(int n)
249 {
250     int nfac[15];
251     int m = 0, j = 0, jj = 0, k = 0, kt = 0, max = 0;
252
253     for (k = abs(n) ; k % 16 == 0 ; k /= 16)
254     {
255         m++;
256         nfac[m - 1] = 4;
257     }
258
259     for (j = 3, jj = 9 ; jj <= k ; j += 2, jj = j * j)
260     {
261         if (k % jj != 0)
262         {
263             continue;
264         }
265         m++;
266         nfac[m - 1] = j;
267         k /= jj;
268     }
269
270     if (k <= 4)
271     {
272         kt = m;
273         nfac[m] = k;
274         if (k != 1)
275         {
276             m++;
277         }
278     }
279     else
280     {
281         if (k % 4 == 0)
282         {
283             m++;
284             nfac[m - 1] = 2;
285             k /= 4;
286         }
287
288         kt = m;
289         for (j = 2 ; j <= k ; j = ((j + 1) / 2) * 2 + 1)
290         {
291             if (k % j != 0)
292             {
293                 continue;
294             }
295
296             m++;
297             nfac[m - 1] = j;
298             k /= j;
299         }
300     }
301
302     if (kt != 0)
303     {
304         for (j = kt ; j > 0 ; j--)
305         {
306             m++;
307             nfac[m - 1] = nfac[j - 1];
308         }
309     }
310
311     /* get nfac maximum */
312     for (j = 0, max = nfac[0]; j < m ; j++)
313     {
314         max = (nfac[j] > max) ? nfac[j] : max;
315     }
316
317     return max;
318 }
319
320 bool isPowerOf2(int _iVal)
321 {
322     return (_iVal & (_iVal - 1)) == 0;
323 }