1 /* NUMfft_d.cpp
2 *
3 * Copyright (C) 1997-2011 David Weenink, Paul Boersma 2016-2018,2020
4 *
5 * This code is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This code is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this work. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 /* djmw 20020813 GPL header
20 djmw 20040511 Added n>1 test for compatibility with old behaviour.
21 djmw 20110308 struct renaming
22 */
23
24 #include "NUM2.h"
25 #include "melder.h"
26
27 #define FFT_DATA_TYPE double
28 #include "NUMfft_core.h"
29
NUMforwardRealFastFourierTransform(VEC data)30 void NUMforwardRealFastFourierTransform (VEC data) {
31 autoNUMfft_Table table;
32 NUMfft_Table_init (& table, data.size);
33 NUMfft_forward (& table, data);
34 if (data.size > 1) {
35 /*
36 To be compatible with old behaviour.
37 */
38 double tmp = data [data.size];
39 for (integer i = data.size; i > 2; i --)
40 data [i] = data [i - 1];
41 data [2] = tmp;
42 }
43 }
44
NUMreverseRealFastFourierTransform(VEC data)45 void NUMreverseRealFastFourierTransform (VEC data) {
46 autoNUMfft_Table table;
47 if (data.size > 1) {
48 /*
49 To be compatible with old behaviour.
50 */
51 double tmp = data [2];
52 for (integer i = 2; i < data.size; i ++)
53 data [i] = data [i + 1];
54 data [data.size] = tmp;
55 }
56 NUMfft_Table_init (& table, data.size);
57 NUMfft_backward (& table, data);
58 }
59
NUMfft_forward(NUMfft_Table me,VEC data)60 void NUMfft_forward (NUMfft_Table me, VEC data) {
61 if (my n == 1)
62 return;
63 Melder_assert (my n == data.size);
64 drftf1 (my n, data.asArgumentToFunctionThatExpectsZeroBasedArray(),
65 my trigcache.asArgumentToFunctionThatExpectsZeroBasedArray(),
66 my trigcache.asArgumentToFunctionThatExpectsZeroBasedArray() + my n,
67 my splitcache.asArgumentToFunctionThatExpectsZeroBasedArray()
68 );
69 }
70
NUMfft_backward(NUMfft_Table me,VEC data)71 void NUMfft_backward (NUMfft_Table me, VEC data) {
72 if (my n == 1)
73 return;
74 Melder_assert (my n == data.size);
75 drftb1 (my n, data.asArgumentToFunctionThatExpectsZeroBasedArray(),
76 my trigcache.asArgumentToFunctionThatExpectsZeroBasedArray(),
77 my trigcache.asArgumentToFunctionThatExpectsZeroBasedArray() + my n,
78 my splitcache.asArgumentToFunctionThatExpectsZeroBasedArray()
79 );
80 }
81
NUMfft_Table_init(NUMfft_Table me,integer n)82 void NUMfft_Table_init (NUMfft_Table me, integer n) {
83 my n = n;
84 my trigcache = zero_VEC (3 * n);
85 my splitcache = zero_INTVEC (32);
86 NUMrffti (n, my trigcache.asArgumentToFunctionThatExpectsZeroBasedArray(),
87 my splitcache.asArgumentToFunctionThatExpectsZeroBasedArray()
88 );
89 }
90
NUMrealft(VEC data,integer isign)91 void NUMrealft (VEC data, integer isign) {
92 if (isign == 1)
93 NUMforwardRealFastFourierTransform (data);
94 else
95 NUMreverseRealFastFourierTransform (data);
96 }
97
98 /* End of file NUMfft.cpp */
99