lib file was not embedded in the binary version
[scilab.git] / scilab / modules / signal_processing / sci_gateway / c / sci_fft.c
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2006 - INRIA - Allan CORNET
4  * Copyright (C) 2009 - Digiteo - Vincent LIARD
5  *
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.  The terms
9  * are also available at
10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11  *
12  */
13
14 #include "api_scilab.h"
15
16
17 #include "gw_signal.h"
18 #include "MALLOC.h"
19 #include "stack-c.h"
20 #include "Scierror.h"
21 #include "localization.h"
22 #include <math.h>
23
24 extern void C2F(fft842)(int *inverse, int *signal_length,
25                         double *signal_real, double *signal_imaginary,
26                         int *error);
27 extern void C2F(dfft2)(double *signal_real, double *signal_imaginary,
28                        int *nseg, int *n, int *nspn,
29                        int *inverse, int *error,
30                        int *buffer_data, int *buffer_size);
31
32 int POW2_15 = 32768;
33
34 int my_log2(int n);
35 int my_pow2(int n);
36 int is_pow2(int n);
37 int maxfactor(int n);
38 int fft_1dim(double *signal_real, double *signal_imaginary, int signal_length,
39              int inverse, int *buffer_data, int buffer_size);
40 int fft_2dim(double *signal_real, double *signal_imaginary,
41              int signal_rows, int signal_cols,
42              int inverse, int *buffer_data, int buffer_size);
43 int fft_ndim(double *signal_real, double *signal_imaginary, int signal_length,
44              int dimensions_length, int dimension_stride,
45              int inverse, int *buffer_data, int buffer_size);
46
47 int sci_fft(char *fname,  int* _piKey)
48 {
49   double
50     *signal_real, *signal_imaginary, /* input signal */
51     *frequencies_real, *frequencies_imaginary, /* freq domain output */
52     *argument;
53   int
54     inverse = -1, /* -1 for forward transform, 1 for inverse */
55     /* [4213] should likely hold an array */
56     dimensions_length = 0, /* length of each dimension */
57     /* [4213] should likely hold an array */
58     /* assert: dimension_stride(1) = product(dimensions_length(j),j=1..i-1) */
59     dimension_stride = 0, /* stride for one step along each dimension */
60     dimension_count = 0; /* number of dimensions */
61   int signal_rows = 0, signal_cols = 0, rows = 0, cols = 0;
62
63   int *buffer_data = NULL;
64   int signal_length = 0, error = 0, buffer_size = 0;
65   int one = 1;
66   int *p;
67
68   CheckRhs(1,4); /* but RHS=3 invalid too: checked afterwars */
69   /** @todo check RHS = 3 verified */
70   CheckLhs(1,1);
71
72   /* collecting and checking arguments */
73   switch (Rhs) {   /* no breaks to collect all arguments */
74   case 4:
75     // GetRhsVarMatrixDouble(4, &rows, &cols, &argument);
76     getVarAddressFromPosition(_piKey, 4, &p);
77     getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
78     dimension_stride = (int)argument[0];
79     // GetRhsVarMatrixDouble(3, &rows, &cols, &argument);
80     getVarAddressFromPosition(_piKey, 3, &p);
81     getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
82     dimensions_length = (int)argument[0];
83     dimension_count = 3; /* any value > 2 (used as a flag) */
84   case 2: /* FALLTHROUGH */
85     // GetRhsVarMatrixDouble(2, &rows, &cols, &argument);
86     getVarAddressFromPosition(_piKey, 2, &p);
87     getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
88     inverse = argument[0];
89     if (rows != 1 || cols != 1) {
90       Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), fname, 2);
91       return 1;
92     }
93     if (inverse != 1 && inverse != -1) {
94       Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"), fname, 2, "-1 1");
95       return 1;
96     }
97   case 1: /* FALLTHROUGH */
98     // GetVarDimension(1, &signal_rows, &signal_cols);
99     getVarAddressFromPosition(_piKey, 1, &p);
100     getVarDimension(_piKey, p, &signal_rows, &signal_cols);
101     dimension_count = Max(dimension_count,
102                           (signal_rows == 1 || signal_cols == 1) ? 1 : 2);
103     signal_length = signal_rows * signal_cols;
104     if (signal_length <= 1) {
105       Scierror(999, _("%s: Wrong size for input argument #%d: A vector expected.\n"), fname, 1);
106       return 1;
107     }
108     // iAllocComplexMatrixOfDouble(Rhs + 1, signal_rows, signal_cols, &frequencies_real, &frequencies_imaginary);
109     allocComplexMatrixOfDouble(_piKey, Rhs + 1, signal_rows, signal_cols, &frequencies_real, &frequencies_imaginary);
110     createComplexMatrixOfDouble(_piKey, Rhs + 1, signal_rows, signal_cols, frequencies_real, frequencies_imaginary);
111     // if (iIsComplex(1)) {
112     if (isVarComplex(_piKey, p)) {
113       // GetRhsVarMatrixComplex(1, &signal_rows, &signal_cols, &signal_real, &signal_imaginary);
114       getComplexMatrixOfDouble(_piKey, p, &signal_rows, &signal_cols, &signal_real, &signal_imaginary);
115       C2F(dcopy)(&signal_length, signal_real,
116                                 &one, frequencies_real, &one);
117       C2F(dcopy)(&signal_length, signal_imaginary,
118                                 &one, frequencies_imaginary, &one);
119     }
120     else {
121       double zero = 0L;
122       // GetRhsVarMatrixDouble(1, &signal_rows, &signal_cols, &signal_real);
123       getMatrixOfDouble(_piKey, p, &signal_rows, &signal_cols, &signal_real);
124       C2F(dcopy)(&signal_length, signal_real,
125                                 &one, frequencies_real, &one);
126       C2F(dset)(&signal_length, &zero, frequencies_imaginary, &one);
127     }
128     break;
129   default:
130     Scierror(999, _("%s: Wrong number of input arguments: %d to %d expected.\n"), fname, 1, 4);
131     return 1;
132   }
133
134   /* upper bound on the workspace required by dfft2 */
135   buffer_size = 8 * maxfactor(signal_length) + 24;
136   buffer_data = (int *)CALLOC(buffer_size, sizeof(int));
137   if (buffer_data == NULL) {
138     Scierror(999, _("%s : Memory allocation error.\n"), fname);
139     return 1;
140   }
141
142   /* how to compute depends on the dimension */
143   switch (dimension_count) {
144   case 1:
145     error = fft_1dim(frequencies_real, frequencies_imaginary,
146                      signal_length, inverse, buffer_data, buffer_size);
147     break;
148   case 2:
149     error = fft_2dim(frequencies_real, frequencies_imaginary,
150                      signal_rows, signal_cols, inverse,
151                      buffer_data, buffer_size);
152     if (error == 1) {
153       Scierror(999, _("%s : Memory allocation error.\n"), fname);
154     }
155     break;
156   default:
157     error = fft_ndim(frequencies_real, frequencies_imaginary,
158                      signal_length, dimensions_length, dimension_stride,
159                      inverse, buffer_data, buffer_size);
160     break;
161   }
162
163   /* preparing the return value */
164   LhsVar(1) = Rhs + 1;
165   PutLhsVar();
166   FREE(buffer_data);
167   buffer_data = NULL;
168   return 0;
169 }
170
171 int fft_1dim(double *signal_real, double *signal_imaginary,
172              int signal_length, int inverse,
173              int *buffer_data, int buffer_size) {
174   int error = 0, one = 1;
175   if (is_pow2(signal_length) && signal_length < POW2_15) {
176     C2F(fft842)(&inverse, &signal_length,
177                 signal_real, signal_imaginary, &error);
178   }
179   else {
180     C2F(dfft2)(signal_real, signal_imaginary,
181                &one, &signal_length, &one,
182                &inverse, &error,
183                buffer_data, &buffer_size);
184   }
185   return error;
186 }
187
188 int fft_2dim(double *signal_real, double *signal_imaginary,
189                     int signal_rows, int signal_cols,
190                     int inverse,
191                     int *buffer_data, int buffer_size) {
192   int error = 0, i, one = 1;
193   if (is_pow2(signal_rows) && signal_rows < POW2_15) {
194     for (i = 0 ; i < signal_cols ; i++) {
195       C2F(fft842)(&inverse, &signal_rows,
196                   &(signal_real[signal_rows * i]),
197                   &(signal_imaginary[signal_rows * i]), &error);
198     }
199   }
200   else {
201     C2F(dfft2)(signal_real, signal_imaginary,
202                &signal_cols, &signal_rows, &one,
203                &inverse, &error, buffer_data, &buffer_size);
204   }
205   if (is_pow2(signal_cols) && signal_cols < POW2_15) {
206     double *temp_real, *temp_imaginary;
207     temp_real = (double *)MALLOC(signal_cols * sizeof(double));
208     temp_imaginary = (double *)MALLOC(signal_cols * sizeof(double));
209     if (temp_real == NULL || temp_imaginary == NULL) {
210       return 1;
211     }
212     for (i = 0 ; i < signal_rows ; i++) {
213       C2F(dcopy)(&signal_cols, &(signal_real[i]),
214                  &signal_rows, temp_real, &one);
215       C2F(dcopy)(&signal_cols, &(signal_imaginary[i]),
216                  &signal_rows, temp_imaginary, &one);
217       C2F(fft842)(&inverse, &signal_cols,
218                   temp_real, temp_imaginary, &error);
219       C2F(dcopy)(&signal_cols, temp_real, &one,
220                  &(signal_real[i]), &signal_rows);
221       C2F(dcopy)(&signal_cols, temp_imaginary, &one,
222                  &(signal_imaginary[i]), &signal_rows);
223     }
224     FREE(temp_imaginary);
225     temp_imaginary = NULL;
226     free(temp_real);
227     temp_real = NULL;
228   }
229   else { /* erroneous implementation suspected */
230     C2F(dfft2)(signal_real, signal_imaginary,
231                &one, &signal_cols, &signal_rows,
232                &inverse, &error, buffer_data, &buffer_size);
233   }
234   return error;
235 }
236
237 int fft_ndim(double *signal_real, double *signal_imaginary,
238              int signal_length, int dimensions_length, int dimension_stride,
239              int inverse, int *buffer_data, int buffer_size) {
240   /* translated litterally from Fortran... but... wtf ?! */
241   int error = 0;
242   int nseg = signal_length / dimensions_length / dimension_stride;
243   C2F(dfft2)(signal_real, signal_imaginary, &nseg,
244              &dimensions_length, &dimension_stride, &inverse, &error,
245              buffer_data, &buffer_size);
246   return error;
247 }
248
249 /* maximum factor of the prime decomposition of n
250    - extracted and translated from dfftbi.f
251    - pityfully needed to allocate the appropriate mem
252    - translated litterally to make sure the result is the same
253    TODO: remove once allocation moved to algorithm
254  */
255 int maxfactor(int n) {
256   int nfac[15];
257   int m = 0, j = 0, jj = 0, k = 0, kt = 0, max = 0;
258
259   for (k = abs(n) ; k % 16 == 0 ; k /= 16) {
260     m++;
261     nfac[m-1] = 4;
262   }
263   for (j = 3, jj = 9 ; jj <= k ; j += 2, jj = j*j) {
264     if (k % jj != 0) continue;
265     m++;
266     nfac[m-1] = j;
267     k /= jj;
268   }
269   if (k <= 4) {
270     kt = m;
271     nfac[m] = k;
272     if (k != 1)
273       m++;
274   }
275   else {
276     if (k % 4 == 0) {
277       m++;
278       nfac[m-1] = 2;
279       k /= 4;
280     }
281     kt = m;
282     for (j = 2 ; j <= k ; j = ((j + 1) / 2) * 2 + 1) {
283       if (k % j != 0)
284         continue;
285       m++;
286       nfac[m-1] = j;
287       k /= j;
288     }
289   }
290   if (kt != 0) {
291     for (j = kt ; j > 0 ; j--) {
292       m++;
293       nfac[m-1] = nfac[j-1];
294     }
295   }
296   /* get nfac maximum */
297   for (j = 0, max = nfac[0]; j < m ; j++) {
298     max = (nfac[j] > max) ? nfac[j] : max;
299   }
300   return max;
301 }
302
303 /* because log2 is only C99... */
304 int my_log2(int n) {
305   return (int)(log(n)/log(2));
306 }
307
308 int my_pow2(int n) {
309   return 1 << n;
310 }
311
312 int is_pow2(int n) {
313   return my_pow2(my_log2(n)) == n;
314 }