1 /*
2 * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
3 * This file is part of KISS FFT - https://github.com/mborgerding/kissfft
4 *
5 * SPDX-License-Identifier: BSD-3-Clause
6 * See COPYING file for more information.
7 */
8
9 #include "kiss_fftr_s16.h"
10 #include "_kiss_fft_guts_s16.h"
11
12 struct kiss_fftr_s16_state
13 {
14 kiss_fft_s16_cfg substate;
15 kiss_fft_s16_cpx *tmpbuf;
16 kiss_fft_s16_cpx *super_twiddles;
17 #ifdef USE_SIMD
18 void *pad;
19 #endif
20 };
21
22 kiss_fftr_s16_cfg
kiss_fftr_s16_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)23 kiss_fftr_s16_alloc (int nfft, int inverse_fft, void *mem, size_t * lenmem)
24 {
25 int i;
26 kiss_fftr_s16_cfg st = NULL;
27 size_t subsize = 0, memneeded;
28
29 g_return_val_if_fail ((nfft & 1) == 0, NULL);
30 nfft >>= 1;
31
32 kiss_fft_s16_alloc (nfft, inverse_fft, NULL, &subsize);
33 memneeded =
34 ALIGN_STRUCT (sizeof (struct kiss_fftr_s16_state)) +
35 ALIGN_STRUCT (subsize) + sizeof (kiss_fft_s16_cpx) * (nfft * 3 / 2);
36
37 if (lenmem == NULL) {
38 st = (kiss_fftr_s16_cfg) KISS_FFT_S16_MALLOC (memneeded);
39 } else {
40 if (*lenmem >= memneeded)
41 st = (kiss_fftr_s16_cfg) mem;
42 *lenmem = memneeded;
43 }
44 if (!st)
45 return NULL;
46
47 st->substate = (kiss_fft_s16_cfg) (((char *) st) + ALIGN_STRUCT (sizeof (struct kiss_fftr_s16_state))); /*just beyond kiss_fftr_s16_state struct */
48 st->tmpbuf =
49 (kiss_fft_s16_cpx *) (((char *) st->substate) + ALIGN_STRUCT (subsize));
50 st->super_twiddles = st->tmpbuf + nfft;
51 kiss_fft_s16_alloc (nfft, inverse_fft, st->substate, &subsize);
52
53 for (i = 0; i < nfft / 2; ++i) {
54 double phase =
55 -3.14159265358979323846264338327 * ((double) (i + 1) / nfft + .5);
56 if (inverse_fft)
57 phase *= -1;
58 kf_cexp (st->super_twiddles + i, phase);
59 }
60 return st;
61 }
62
63 void
kiss_fftr_s16(kiss_fftr_s16_cfg st,const kiss_fft_s16_scalar * timedata,kiss_fft_s16_cpx * freqdata)64 kiss_fftr_s16 (kiss_fftr_s16_cfg st, const kiss_fft_s16_scalar * timedata,
65 kiss_fft_s16_cpx * freqdata)
66 {
67 /* input buffer timedata is stored row-wise */
68 int k, ncfft;
69 kiss_fft_s16_cpx fpnk, fpk, f1k, f2k, tw, tdc;
70
71 g_return_if_fail (!st->substate->inverse);
72
73 ncfft = st->substate->nfft;
74
75 /*perform the parallel fft of two real signals packed in real,imag */
76 kiss_fft_s16 (st->substate, (const kiss_fft_s16_cpx *) timedata, st->tmpbuf);
77 /* The real part of the DC element of the frequency spectrum in st->tmpbuf
78 * contains the sum of the even-numbered elements of the input time sequence
79 * The imag part is the sum of the odd-numbered elements
80 *
81 * The sum of tdc.r and tdc.i is the sum of the input time sequence.
82 * yielding DC of input time sequence
83 * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1...
84 * yielding Nyquist bin of input time sequence
85 */
86
87 tdc.r = st->tmpbuf[0].r;
88 tdc.i = st->tmpbuf[0].i;
89 C_FIXDIV (tdc, 2);
90 CHECK_OVERFLOW_OP (tdc.r, +, tdc.i);
91 CHECK_OVERFLOW_OP (tdc.r, -, tdc.i);
92 freqdata[0].r = tdc.r + tdc.i;
93 freqdata[ncfft].r = tdc.r - tdc.i;
94 #ifdef USE_SIMD
95 freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps (0);
96 #else
97 freqdata[ncfft].i = freqdata[0].i = 0;
98 #endif
99
100 for (k = 1; k <= ncfft / 2; ++k) {
101 fpk = st->tmpbuf[k];
102 fpnk.r = st->tmpbuf[ncfft - k].r;
103 fpnk.i = -st->tmpbuf[ncfft - k].i;
104 C_FIXDIV (fpk, 2);
105 C_FIXDIV (fpnk, 2);
106
107 C_ADD (f1k, fpk, fpnk);
108 C_SUB (f2k, fpk, fpnk);
109 C_MUL (tw, f2k, st->super_twiddles[k - 1]);
110
111 freqdata[k].r = HALF_OF (f1k.r + tw.r);
112 freqdata[k].i = HALF_OF (f1k.i + tw.i);
113 freqdata[ncfft - k].r = HALF_OF (f1k.r - tw.r);
114 freqdata[ncfft - k].i = HALF_OF (tw.i - f1k.i);
115 }
116 }
117
118 void
kiss_fftri_s16(kiss_fftr_s16_cfg st,const kiss_fft_s16_cpx * freqdata,kiss_fft_s16_scalar * timedata)119 kiss_fftri_s16 (kiss_fftr_s16_cfg st, const kiss_fft_s16_cpx * freqdata,
120 kiss_fft_s16_scalar * timedata)
121 {
122 /* input buffer timedata is stored row-wise */
123 int k, ncfft;
124
125 g_return_if_fail (st->substate->inverse);
126
127 ncfft = st->substate->nfft;
128
129 st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r;
130 st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r;
131 C_FIXDIV (st->tmpbuf[0], 2);
132
133 for (k = 1; k <= ncfft / 2; ++k) {
134 kiss_fft_s16_cpx fk, fnkc, fek, fok, tmp;
135 fk = freqdata[k];
136 fnkc.r = freqdata[ncfft - k].r;
137 fnkc.i = -freqdata[ncfft - k].i;
138 C_FIXDIV (fk, 2);
139 C_FIXDIV (fnkc, 2);
140
141 C_ADD (fek, fk, fnkc);
142 C_SUB (tmp, fk, fnkc);
143 C_MUL (fok, tmp, st->super_twiddles[k - 1]);
144 C_ADD (st->tmpbuf[k], fek, fok);
145 C_SUB (st->tmpbuf[ncfft - k], fek, fok);
146 #ifdef USE_SIMD
147 st->tmpbuf[ncfft - k].i *= _mm_set1_ps (-1.0);
148 #else
149 st->tmpbuf[ncfft - k].i *= -1;
150 #endif
151 }
152 kiss_fft_s16 (st->substate, st->tmpbuf, (kiss_fft_s16_cpx *) timedata);
153 }
154