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.h"
10 #include "_kiss_fft_guts.h"
11
12 struct kiss_fftr_state{
13 kiss_fft_cfg substate;
14 kiss_fft_cpx * tmpbuf;
15 kiss_fft_cpx * super_twiddles;
16 #ifdef USE_SIMD
17 void * pad;
18 #endif
19 };
20
kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)21 kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)
22 {
23 KISS_FFT_ALIGN_CHECK(mem)
24
25 int i;
26 kiss_fftr_cfg st = NULL;
27 size_t subsize = 0, memneeded;
28
29 if (nfft & 1) {
30 KISS_FFT_ERROR("Real FFT optimization must be even.");
31 return NULL;
32 }
33 nfft >>= 1;
34
35 kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize);
36 memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2);
37
38 if (lenmem == NULL) {
39 st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded);
40 } else {
41 if (*lenmem >= memneeded)
42 st = (kiss_fftr_cfg) mem;
43 *lenmem = memneeded;
44 }
45 if (!st)
46 return NULL;
47
48 st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */
49 st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize);
50 st->super_twiddles = st->tmpbuf + nfft;
51 kiss_fft_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
kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar * timedata,kiss_fft_cpx * freqdata)63 void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata)
64 {
65 /* input buffer timedata is stored row-wise */
66 int k,ncfft;
67 kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc;
68
69 if ( st->substate->inverse) {
70 KISS_FFT_ERROR("kiss fft usage error: improper alloc");
71 return;/* The caller did not call the correct function */
72 }
73
74 ncfft = st->substate->nfft;
75
76 /*perform the parallel fft of two real signals packed in real,imag*/
77 kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf );
78 /* The real part of the DC element of the frequency spectrum in st->tmpbuf
79 * contains the sum of the even-numbered elements of the input time sequence
80 * The imag part is the sum of the odd-numbered elements
81 *
82 * The sum of tdc.r and tdc.i is the sum of the input time sequence.
83 * yielding DC of input time sequence
84 * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1...
85 * yielding Nyquist bin of input time sequence
86 */
87
88 tdc.r = st->tmpbuf[0].r;
89 tdc.i = st->tmpbuf[0].i;
90 C_FIXDIV(tdc,2);
91 CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i);
92 CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i);
93 freqdata[0].r = tdc.r + tdc.i;
94 freqdata[ncfft].r = tdc.r - tdc.i;
95 #ifdef USE_SIMD
96 freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0);
97 #else
98 freqdata[ncfft].i = freqdata[0].i = 0;
99 #endif
100
101 for ( k=1;k <= ncfft/2 ; ++k ) {
102 fpk = st->tmpbuf[k];
103 fpnk.r = st->tmpbuf[ncfft-k].r;
104 fpnk.i = - st->tmpbuf[ncfft-k].i;
105 C_FIXDIV(fpk,2);
106 C_FIXDIV(fpnk,2);
107
108 C_ADD( f1k, fpk , fpnk );
109 C_SUB( f2k, fpk , fpnk );
110 C_MUL( tw , f2k , st->super_twiddles[k-1]);
111
112 freqdata[k].r = HALF_OF(f1k.r + tw.r);
113 freqdata[k].i = HALF_OF(f1k.i + tw.i);
114 freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r);
115 freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i);
116 }
117 }
118
kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx * freqdata,kiss_fft_scalar * timedata)119 void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
120 {
121 /* input buffer timedata is stored row-wise */
122 int k, ncfft;
123
124 if (st->substate->inverse == 0) {
125 KISS_FFT_ERROR("kiss fft usage error: improper alloc");
126 return;/* The caller did not call the correct function */
127 }
128
129 ncfft = st->substate->nfft;
130
131 st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r;
132 st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r;
133 C_FIXDIV(st->tmpbuf[0],2);
134
135 for (k = 1; k <= ncfft / 2; ++k) {
136 kiss_fft_cpx fk, fnkc, fek, fok, tmp;
137 fk = freqdata[k];
138 fnkc.r = freqdata[ncfft - k].r;
139 fnkc.i = -freqdata[ncfft - k].i;
140 C_FIXDIV( fk , 2 );
141 C_FIXDIV( fnkc , 2 );
142
143 C_ADD (fek, fk, fnkc);
144 C_SUB (tmp, fk, fnkc);
145 C_MUL (fok, tmp, st->super_twiddles[k-1]);
146 C_ADD (st->tmpbuf[k], fek, fok);
147 C_SUB (st->tmpbuf[ncfft - k], fek, fok);
148 #ifdef USE_SIMD
149 st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0);
150 #else
151 st->tmpbuf[ncfft - k].i *= -1;
152 #endif
153 }
154 kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata);
155 }
156