Mercurial > hg > dmlib

comparison src/dmfft.c @ 1691:e68f7b0cb536

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleanups.
author Matti Hamalainen <ccr@tnsp.org>
date Tue, 05 Jun 2018 09:37:37 +0300
parents f077a709d491
children 69a5af2eb1ea
comparison
equal deleted inserted replaced
1690:12d4da502fa4 1691:e68f7b0cb536
11 */ 11 */
12 #include "dmfft.h" 12 #include "dmfft.h"
13 #include <math.h> 13 #include <math.h>
14 14
15 15
16 #define DM_PW(a, b, c) (sqrt((a) * (a) + (b) * (b)) * (c))
17
18
16 int dmInitializeFFT(DMFFTContext * ctx, const int fftlen) 19 int dmInitializeFFT(DMFFTContext * ctx, const int fftlen)
17 { 20 {
18 int i;
19
20 if (ctx == NULL) 21 if (ctx == NULL)
21 return DMERR_NULLPTR; 22 return DMERR_NULLPTR;
22 23
23 if (fftlen <= 16) 24 if (fftlen <= 16)
24 return DMERR_INVALID_ARGS; 25 return DMERR_INVALID_ARGS;
36 dmFree(ctx->sinTable); 37 dmFree(ctx->sinTable);
37 return DMERR_MALLOC; 38 return DMERR_MALLOC;
38 } 39 }
39 40
40 // Calculate data 41 // Calculate data
41 for (i = 0; i < ctx->npoints; i++) 42 for (int i = 0; i < ctx->npoints; i++)
42 { 43 {
43 int mask, tmp; 44 int mask, tmp;
44 for (tmp = 0, mask = ctx->npoints / 2; mask > 0; mask >>= 1) 45 for (tmp = 0, mask = ctx->npoints / 2; mask > 0; mask >>= 1)
45 { 46 {
46 tmp = (tmp >> 1) + ((i & mask) ? ctx->npoints : 0); 47 tmp = (tmp >> 1) + ((i & mask) ? ctx->npoints : 0);
47 } 48 }
48 49
49 ctx->breversed[i] = tmp; 50 ctx->breversed[i] = tmp;
50 } 51 }
51 52
52 for (i = 0; i < ctx->npoints; i++) 53 for (int i = 0; i < ctx->npoints; i++)
53 { 54 {
54 ctx->sinTable[ctx->breversed[i]] = -sin((2.0f * DM_PI * i) / fftlen); 55 ctx->sinTable[ctx->breversed[i]] = -sin((2.0f * DM_PI * i) / fftlen);
55 ctx->sinTable[ctx->breversed[i] + 1] = 56 ctx->sinTable[ctx->breversed[i] + 1] =
56 -cos((2.0f * DM_PI * i) / fftlen); 57 -cos((2.0f * DM_PI * i) / fftlen);
57 } 58 }
109 110
110 rounds >>= 1; 111 rounds >>= 1;
111 } 112 }
112 113
113 // Massage output to get the output for a real input sequence. 114 // Massage output to get the output for a real input sequence.
114
115 while (br1 <= br2) 115 while (br1 <= br2)
116 { 116 {
117 const DMFFTType vsin = ctx->sinTable[*br1]; 117 const DMFFTType vsin = ctx->sinTable[*br1];
118 const DMFFTType vcos = ctx->sinTable[*br1 + 1]; 118 const DMFFTType vcos = ctx->sinTable[*br1 + 1];
119 DMFFTType *A = buffer + *br1, 119 DMFFTType *A = buffer + *br1,
150 150
151 151
152 int dmConvertFFTtoFreqDomain(DMFFTContext *ctx, DMFFTType *buffer, 152 int dmConvertFFTtoFreqDomain(DMFFTContext *ctx, DMFFTType *buffer,
153 DMFFTType *real, DMFFTType *imag) 153 DMFFTType *real, DMFFTType *imag)
154 { 154 {
155 int i; 155 if (ctx == NULL || buffer == NULL ||
156 if (ctx == NULL || buffer == NULL || real == NULL || imag == NULL) 156 real == NULL || imag == NULL)
157 return DMERR_NULLPTR; 157 return DMERR_NULLPTR;
158 158
159 // Convert from bitreversed form to normal frequency domain 159 // Convert from bitreversed form to normal frequency domain
160 const int *breversed = ctx->breversed; 160 const int *breversed = ctx->breversed;
161 161
162 for (i = 1; i < ctx->npoints; i++) 162 for (int i = 1; i < ctx->npoints; i++)
163 { 163 {
164 real[i] = buffer[breversed[i] ]; 164 real[i] = buffer[breversed[i] ];
165 imag[i] = buffer[breversed[i]+1]; 165 imag[i] = buffer[breversed[i]+1];
166 } 166 }
167 167
173 173
174 return DMERR_OK; 174 return DMERR_OK;
175 } 175 }
176 176
177 177
178 #define DM_PW(a, b, c) sqrt(a * a + b * b) * c
179
180
181 int dmConvertFFTtoFreqAndPower(DMFFTContext *ctx, DMFFTType *buffer, 178 int dmConvertFFTtoFreqAndPower(DMFFTContext *ctx, DMFFTType *buffer,
182 DMFFTType *real, DMFFTType *imag, DMFFTType *power, const DMFFTType scale) 179 DMFFTType *real, DMFFTType *imag, DMFFTType *power, const DMFFTType scale)
183 { 180 {
184 int i; 181 if (ctx == NULL || buffer == NULL ||
185 if (ctx == NULL || buffer == NULL || real == NULL || imag == NULL || power == NULL) 182 real == NULL || imag == NULL || power == NULL)
186 return DMERR_NULLPTR; 183 return DMERR_NULLPTR;
187 184
188 // Convert from bitreversed form to normal frequency domain 185 // Convert from bitreversed form to normal frequency domain
189 const int *breversed = ctx->breversed; 186 const int *breversed = ctx->breversed;
190 187
191 for (i = 1; i < ctx->npoints; i++) 188 for (int i = 1; i < ctx->npoints; i++)
192 { 189 {
193 DMFFTType fre = real[i] = buffer[breversed[i] ], 190 DMFFTType fre = real[i] = buffer[breversed[i] ],
194 fim = imag[i] = buffer[breversed[i]+1]; 191 fim = imag[i] = buffer[breversed[i]+1];
195 192
196 power[i] = DM_PW(fre, fim, scale); 193 power[i] = DM_PW(fre, fim, scale);
210 207
211 208
212 int dmConvertFFTtoPowerAndSum(DMFFTContext *ctx, DMFFTType *buffer, 209 int dmConvertFFTtoPowerAndSum(DMFFTContext *ctx, DMFFTType *buffer,
213 DMFFTType *power, const DMFFTType pscale, DMFFTType *sum, const DMFFTType sscale) 210 DMFFTType *power, const DMFFTType pscale, DMFFTType *sum, const DMFFTType sscale)
214 { 211 {
215 int i;
216 if (ctx == NULL || buffer == NULL || power == NULL || sum == NULL) 212 if (ctx == NULL || buffer == NULL || power == NULL || sum == NULL)
217 return DMERR_NULLPTR; 213 return DMERR_NULLPTR;
218 214
219 // Convert from bitreversed form to normal frequency domain 215 // Convert from bitreversed form to normal frequency domain
220 const int *breversed = ctx->breversed; 216 const int *breversed = ctx->breversed;
221 DMFFTType psum = 0; 217 DMFFTType psum = 0;
222 218
223 for (i = 1; i < ctx->npoints; i++) 219 for (int i = 1; i < ctx->npoints; i++)
224 { 220 {
225 DMFFTType fre = buffer[breversed[i] ], 221 DMFFTType fre = buffer[breversed[i] ],
226 fim = buffer[breversed[i]+1], 222 fim = buffer[breversed[i]+1],
227 fpw = sqrt(fre * fre + fim * fim); 223 fpw = sqrt(fre * fre + fim * fim);
228 224