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 | |
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 | |
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 | |
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 | |