libdspl-2.0
fft.c
1 /*
2 * Copyright (c) 2015-2018 Sergey Bakhurin
3 * Digital Signal Processing Library [http://dsplib.org]
4 *
5 * This file is part of libdspl-2.0.
6 *
7 * is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * DSPL is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with Foobar. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #include <stdlib.h>
22 #include <string.h>
23 #include "dspl.h"
24 
25 int fft_bit_reverse(complex_t* x, complex_t* y, int n, int p2);
26 
27 int fft_dit(fft_t *pfft, int n, complex_t* y);
28 
29 void fft_dit_krn(complex_t *x0, complex_t *x1, complex_t *w, int n,
30  complex_t *y0, complex_t *y1);
31 
32 int fft_p2(int n);
33 
34 
35 
36 
37 
38 /*******************************************************************************
39 COMPLEX vector IFFT
40 *******************************************************************************/
41 int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
42 {
43  int err, k;
44  double norm;
45 
46  if(!x || !pfft || !y)
47  return ERROR_PTR;
48  if(n<1)
49  return ERROR_SIZE;
50 
51 
52  err = fft_create(pfft, n);
53  if(err != RES_OK)
54  return err;
55 
56  memcpy(pfft->t1, x, n*sizeof(complex_t));
57  for(k = 0; k < n; k++)
58  IM(pfft->t1[k]) = -IM(pfft->t1[k]);
59 
60  err = fft_dit(pfft, n, y);
61 
62  if(err!=RES_OK)
63  return err;
64 
65  norm = 1.0 / (double)n;
66  for(k = 0; k < n; k++)
67  {
68  RE(y[k]) = RE(y[k])*norm;
69  IM(y[k]) = -IM(y[k])*norm;
70  }
71  return RES_OK;
72 }
73 
74 
75 
76 
77 
78 
79 
80 /*******************************************************************************
81 Real vector FFT
82 *******************************************************************************/
83 int DSPL_API fft(double *x, int n, fft_t* pfft, complex_t* y)
84 {
85  int err;
86 
87  if(!x || !pfft || !y)
88  return ERROR_PTR;
89  if(n<1)
90  return ERROR_SIZE;
91 
92 
93  err = fft_create(pfft, n);
94  if(err != RES_OK)
95  return err;
96 
97  re2cmplx(x, n, pfft->t1);
98 
99  return fft_dit(pfft, n, y);
100 }
101 
102 
103 
104 
105 /*******************************************************************************
106 COMPLEX vector FFT
107 *******************************************************************************/
108 int DSPL_API fft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
109 {
110  int err;
111 
112  if(!x || !pfft || !y)
113  return ERROR_PTR;
114  if(n<1)
115  return ERROR_SIZE;
116 
117 
118  err = fft_create(pfft, n);
119  if(err != RES_OK)
120  return err;
121 
122  memcpy(pfft->t1, x, n*sizeof(complex_t));
123 
124  return fft_dit(pfft, n, y);
125 }
126 
127 
128 
129 /*******************************************************************************
130 FFT bit reverse
131 *******************************************************************************/
132 int fft_bit_reverse(complex_t* x, complex_t* y, int n, int p2)
133 {
134  static unsigned char rb_table[256] =
135  {
136  0x00, 0x80, 0x40, 0xC0, 0x20, 0xA0, 0x60, 0xE0,
137  0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0,
138  0x08, 0x88, 0x48, 0xC8, 0x28, 0xA8, 0x68, 0xE8,
139  0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8,
140  0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4,
141  0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4,
142  0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC,
143  0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC,
144  0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2,
145  0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2,
146  0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA,
147  0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA,
148  0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6,
149  0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xF6,
150  0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE,
151  0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE,
152  0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1,
153  0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1,
154  0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9,
155  0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9,
156  0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5,
157  0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5,
158  0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0xED,
159  0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD,
160  0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3,
161  0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3,
162  0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB,
163  0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB,
164  0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7,
165  0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7,
166  0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF,
167  0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF
168  };
169 
170  if(!x || !y)
171  return ERROR_PTR;
172  if(n<1 || p2 < 1)
173  return ERROR_SIZE;
174 
175  unsigned int v, c;
176 
177  for(v = 0; v < (unsigned int)n; v++)
178  {
179  c = (unsigned int) ((rb_table[ v & 0xff] << 24) |
180  (rb_table[(v >> 8) & 0xff] << 16) |
181  (rb_table[(v >> 16) & 0xff] << 8) |
182  (rb_table[(v >> 24) & 0xff])) >> (32 - p2);
183 
184  RE(y[c]) = RE(x[v]);
185  IM(y[c]) = IM(x[v]);
186  }
187  return RES_OK;
188 }
189 
190 
191 
192 
193 
194 /*******************************************************************************
195 FFT create
196 *******************************************************************************/
197 int DSPL_API fft_create(fft_t *pfft, int n)
198 {
199  int p2, k,r,m,addr,s;
200  double phi;
201 
202  p2 = fft_p2(n);
203  if(p2 < 1)
204  return ERROR_FFT_SIZE;
205  if(n < pfft->n+1)
206  return RES_OK;
207 
208  pfft->n = n;
209  pfft->p2 = p2;
210 
211  pfft->w = pfft->w ? (complex_t*) realloc(pfft->w, n*sizeof(complex_t)):
212  (complex_t*) malloc( n*sizeof(complex_t));
213 
214 
215  pfft->t0 = pfft->t0 ? (complex_t*) realloc(pfft->t0, n*sizeof(complex_t)):
216  (complex_t*) malloc( n*sizeof(complex_t));
217 
218 
219  pfft->t1 = pfft->t1 ? (complex_t*) realloc(pfft->t1, n*sizeof(complex_t)):
220  (complex_t*) malloc( n*sizeof(complex_t));
221 
222  m = 0;
223  addr = 0;
224 
225  for(k = 0; k < p2; k++)
226  {
227  s = 1<<m;
228  for( r = 0; r < s; r++)
229  {
230  phi = -M_2PI *(double)r / (double)(2*s);
231  RE(pfft->w[addr+r]) = cos(phi);
232  IM(pfft->w[addr+r]) = sin(phi);
233  }
234  addr+=s;
235  m++;
236  }
237  return RES_OK;
238 }
239 
240 
241 
242 /*******************************************************************************
243 FFT decimation in time
244 *******************************************************************************/
245 int fft_dit(fft_t *pfft, int n, complex_t* y)
246 {
247  int k,s,m,waddr, dm,p2;
248  complex_t *t, *t0, *t1, *w;
249  int err;
250 
251  p2 = fft_p2(n);
252  if(p2<0)
253  return ERROR_FFT_SIZE;
254 
255  t0 = pfft->t0;
256  t1 = pfft->t1;
257  w = pfft->w;
258 
259  s = n>>1;
260  m = 1;
261  waddr = 0;
262 
263  err = fft_bit_reverse(t1, t0, n, p2);
264  if(err!= RES_OK)
265  return err;
266 
267  while(s)
268  {
269  dm = m<<1;
270  if(s > 1)
271  {
272  for(k = 0; k < n; k+=dm)
273  {
274  fft_dit_krn(t0+k, t0+k+m, w+waddr, m, t1+k, t1+k+m);
275  }
276  t = t1;
277  t1 = t0;
278  t0 = t;
279  waddr+=m;
280  m <<= 1;
281  }
282  else
283  {
284  fft_dit_krn(t0, t0+m, w+waddr, m, y, y+m);
285  }
286  s >>= 1;
287  }
288  return RES_OK;
289 }
290 
291 
292 
293 
294 
295 /*******************************************************************************
296 FFT decimation in time kernel
297 *******************************************************************************/
298 void fft_dit_krn(complex_t *x0, complex_t *x1, complex_t *w, int n,
299  complex_t *y0, complex_t *y1)
300 {
301  int k;
302  complex_t mul;
303  for(k = 0; k < n; k++)
304  {
305  RE(mul) = CMRE(x1[k], w[k]);
306  IM(mul) = CMIM(x1[k], w[k]);
307 
308  RE(y0[k]) = RE(x0[k]) + RE(mul);
309  IM(y0[k]) = IM(x0[k]) + IM(mul);
310 
311  RE(y1[k]) = RE(x0[k]) - RE(mul);
312  IM(y1[k]) = IM(x0[k]) - IM(mul);
313  }
314 }
315 
316 
317 
318 
319 
320 /*******************************************************************************
321 FFT free
322 *******************************************************************************/
323 void DSPL_API fft_free(fft_t *pfft)
324 {
325  if(!pfft)
326  return;
327  if(pfft->w)
328  free(pfft->w);
329  if(pfft->t0)
330  free(pfft->t0);
331  if(pfft->t1)
332  free(pfft->t1);
333 }
334 
335 
336 
337 
338 
339 /*******************************************************************************
340 FFT power of 2
341 *******************************************************************************/
342 int fft_p2(int n)
343 {
344  int p = (n>0) ? n : 0;
345  int p2 = 0;
346  while((p = p>>1))
347  p2++;
348  if((1<<p2)!=n)
349  return -1;
350  return p2;
351 }
352 
353 
354 
355 
356 
357 
358 
359 /*******************************************************************************
360 FFT shifting
361 *******************************************************************************/
362 int DSPL_API fft_shift(double* x, int n, double* y)
363 {
364  int n2, r;
365  int k;
366  double tmp;
367  double *buf;
368 
369  if(!x || !y)
370  return ERROR_PTR;
371 
372  if(n<1)
373  return ERROR_SIZE;
374 
375  r = n%2;
376  if(!r)
377  {
378  n2 = n>>1;
379  for(k = 0; k < n2; k++)
380  {
381  tmp = x[k];
382  y[k] = x[k+n2];
383  y[k+n2] = tmp;
384  }
385  }
386  else
387  {
388  n2 = (n-1) >> 1;
389  buf = (double*) malloc(n2*sizeof(double));
390  memcpy(buf, x, n2*sizeof(double));
391  memcpy(y, x+n2, (n2+1)*sizeof(double));
392  memcpy(y+n2+1, buf, n2*sizeof(double));
393  free(buf);
394  }
395  return RES_OK;
396 }
397 
398 
399 
400 /*******************************************************************************
401 FFT shifting for complex vector
402 *******************************************************************************/
403 int DSPL_API fft_shift_cmplx(complex_t* x, int n, complex_t* y)
404 {
405  int n2, r;
406  int k;
407  complex_t tmp;
408  complex_t *buf;
409 
410  if(!x || !y)
411  return ERROR_PTR;
412 
413  if(n<1)
414  return ERROR_SIZE;
415 
416  r = n%2;
417  if(!r)
418  {
419  n2 = n>>1;
420  for(k = 0; k < n2; k++)
421  {
422  RE(tmp) = RE(x[k]);
423  IM(tmp) = IM(x[k]);
424 
425  RE(y[k]) = RE(x[k+n2]);
426  IM(y[k]) = IM(x[k+n2]);
427 
428  RE(y[k+n2]) = RE(tmp);
429  IM(y[k+n2]) = IM(tmp);
430  }
431  }
432  else
433  {
434  n2 = (n-1) >> 1;
435  buf = (complex_t*) malloc(n2*sizeof(complex_t));
436  memcpy(buf, x, n2*sizeof(complex_t));
437  memcpy(y, x+n2, (n2+1)*sizeof(complex_t));
438  memcpy(y+n2+1, buf, n2*sizeof(complex_t));
439  free(buf);
440  }
441  return RES_OK;
442 }
443 
complex_t * w
Definition: dspl.h:46
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:123
Структура данных объекта быстрого преобразования Фурье
Definition: dspl.h:44
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:40
#define ERROR_PTR
Ошибка указателя. Данная ошибка означает, что один из обязательных указателей (память под который дол...
Definition: dspl.h:116
complex_t * t1
Definition: dspl.h:48
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:65
#define RE(x)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:64
int fft_shift(double *x, int n, double *y)
Перестановка спектральных отсчетов дискретного преобразования Фурье
Definition: fft.c:362
complex_t * t0
Definition: dspl.h:47
int re2cmplx(double *x, int n, complex_t *y)
Преобразование массива вещественных данных в массив комплексных данных.
Definition: complex.c:82
int n
Definition: dspl.h:49
void fft_free(fft_t *pfft)
Очистить структуру fft_t алгоритма БПФ
Definition: fft.c:323
int p2
Definition: dspl.h:50
#define ERROR_FFT_SIZE
Неверно задан размер БПФ.
Definition: dspl.h:93
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:80
int fft_create(fft_t *pfft, int n)
Заполнение структуры fft_t для алгоритма БПФ
Definition: fft.c:197